Airspace Congestion Smoothing by Stochastic Optimization Daniel Delahaye Amedeo Odoni C.E.N.A/M.I.T M.I.Ty August 15 1996

Abstract This paper addresses the general time-route assignment problem : One considers an air transportation network in a 2 dimensional space with asymmetric non-separable link cost and a eet of aircraft with their associated route and slot of departure. For each ight a set of alternative routes and a set of possible slots of departure are de ned. One must nd \optimal" route and slot allocation for each aircraft in a way that signi caly reduces the peak of workload in the most congested sectors and in the most congested airports, during one day of trac. A state of the art of the existing methods shows that this problem is usually partially treated and the whole problem remains unsolved due to the complexity induced. Genetic algorithm are then presented and adapted to the problem. New recombination operators are described and really improve the quality of the given solutions. Results on a large instance of the problem are then given and show how well genetic algorithm can manage this kind of combinatorial optimization problem.

Topics: Genetic Algorithm, Time of Departure and Route Optimization. Domain Area: Air Trac Control Status: Operational mock-up. Eort: 1 man/year Impact: Reduce the airspace congestion by changing the route and the time of departure of aircraft.

1 Introduction When joining two airports, an aircraft must follow routes and beacons ; these beacons are necessary for pilots to know their position during navigation and because of the small number of beacons on the ground they often represent crossing points of dierent airways. Crossing points may generate con icts between aircraft when their trajectories converge on it at the same time and induce a risk of collision. At the dawn of civil aviation, pilots solved con icts themselves because they always ew in good weather conditions (good visibility) with low speed aircraft. On the other hand, modern jet aircraft do not enable pilots to solve con icts because of their high speed and their ability to y with bad visibility. Therefore, pilots must be Centre d'Etudes de la Navigation Aerienne yMassachusett Institute of Technology

helped by an air trac controller on the ground who has a global view of the current trac distribution in the airspace and can give orders to the pilots to avoid collisions. As there are many aircraft simultaneously present in the sky, a single controller is not able to manage all of them. Then, airspace is partitioned into dierent sectors, each of them being assigned to a controller. Before taking o, when a ight is decided, the pilot has to le a Flight Plan (FPL) which describe the main steps of the ight (slot of departure, ight path, speeds, altitudes etc...), this FPL is necessary for the air trac control system to anticipate the ight monitoring. Usually, when there are few aircraft in the airspace the ight plan is well respected. As any human being, a controller has working limits, and when the number of aircraft increases, some parts of the airspace reach this limit and become congested. Two kinds of congestion can be identi ed according to the part of airspace involved : terminal congestion (around airports) and en route congestion (between airports). In the past, the rst way to reduce these congestions was to modify the structure of the airspace in a way that increases the capacity (increasing the number of runways, increasing the number of sectors by reducing their size). This has a limit due to the cost involved by new runways and the way to manage trac in too small sectors (a controller needs a minimum amount of airspace to be able to solve con icts). The other way to reduce congestion is to modify the ight plans in a way to adapt the demand to the available capacity. To reach this last goal three kinds of action can be applied :

Short term planning : Ground Holding, speed regulation and dynamic re-routing. Nowadays, when a delay

is expected to the arrival airport, a ground delay is applied in a way that will minimize the arrival airborne delay. This ground delay is safer (fewer aircraft waiting in the sky) and cheaper (according to the fuel consumption) than an airborne one. When the aircraft has taken o, this action can be completed by dynamic re-routing and speed regulation to adjust the actual demand to the present capacity. Long term planning : Modi cation of the slot of departure and the routes (static re-routing) of the whole

eet of aircraft. When we have a look at the congestion in the airspace one can remark that the peaks of congestion happen approximately at the same times during the day (one peak in the morning and one peak in the evening). The current time scheduling is usually xed by the airline companies without taking into account the ATC constraints. Then congestion is expected to be reduced by moving (in a limited domain) the time of departure of aircraft (in the past and in the future) and by changing the current ight paths (without too much extradistance). Slot and route fees. Another way to adjust the demand to the available capacity is to apply the principle : \Scarce resources are expensive". Then the departure slots and the arrival slots would have a higher price when the demand is strong and in a similar way the fastest route would have a higher price too. This action might help to spread the demand over time and over space in an equitable way.

This paper is focused on long term planning and shows how stochastic optimization is able to manage this kind of problem. In a rst part, a description of the previous related works is given and in a second one, a simpli ed model is developed and a mathematical formulation of our problem is given. In the third part a description of stochastic optimization is given and particularly of a genetic algorithm. The application of the algorithm is given on the fourth part and nally an improvement is given in the fth part.

2 Previous Related Works

2.1 Introduction

Trac assignment techniques have been developed in a way to reduce congestion in transportation networks by spreading the trac demand in time and in space. Historically, those techniques have been applied to road trac assignment due to the strong level of congestion encountered in this domain. The problem of users of a congested transportation network seeking to determine travel paths of minimal cost from origins to their respective destinations is a classical network problem. An equilibrium occurs when the number of trips between an origin and a destination equals the associated travel demand. Wardrop [93] stated the trac equilibrium conditions through two principles :

First Principle : The journey times of all routes actually used are equal, and less than those which would

be experienced by a single vehicle on any unused route. Second Principle : The average journey time is minimal.

From those two principles, Dafermos and Sparrow [26] coined the terms user-optimized and system-optimized transportation networks to distinguish between two distinct situations in which users act unilaterally, in their own self-interest, in selecting their routes, and in which users select routes according to what is optimal from the societal point of view, in that the total costs in the system are minimized. More detailed descriptions of those equilibria can be found in the following references : [38, 33, 86]. According to the separability1 of the link cost function, dierent algorithms have been developed and the associated references can be summarized in the following array :

symmetric non-separable link cost

User-optimized

System-optimized

[80, 88, 67]

[47, 43]

asymmetric non-separable link cost [34, 35, 85, 24, 48, 59]

[47, 43]

Improved trac assignment algorithms, take into account that demand can be in uenced by the induced cost but still with a static demand over time [41, 42, 25, 74]. By de nition all those techniques, are applied to static trac demand and are mainly used to optimize trac on a long time period and can only capture the macroscopic events. When a more precise matching between trac demand and capacity has to be found, microscopic events have to be taken into account, and dynamic trac assignment techniques have to be used ([78] gives a good description of those techniques). Those techniques try to nd an optimal route, or an optimal time of departure, or both, for each individual in a way to reach a dynamic user-equilibrium or a dynamic system-equilibrium. The complexity induced by the dynamic trac assignment is strong, especially when route and time of departure are simultaneously optimized. This problem is NP HARD [30, 11] and may have several optima [44]. Even when only the route decision is investigated, this problem does not have an exact solution when asymmetric non-separable link costs are used. Due to this complexity, the problem is always partially solved for simpli ed instances :

Ben Akiva, M [9] : \1" origin and \1"destination : time decisions only ; Hendrickson, C et al. [50] : \n" origins \1" destination, system-equilibrium, analytical approach with

simple link cost (time decisions only)) ; Mahmassani, H et al. [62, 22, 55] : \1" origin and \1" destination analytical approach with simple link cost (time decisons only) : Merchant, D.K et al. [69, 70] : network with separable link cost, route choices only (heuristic) ; Alfa, A.S [3] : heuristic for time-route trac assignment, the \convergence" to a xed point is not ensured ; Arnott, R et al. [7] : N parallel routes, time-route decisions, separable link cost function ; HO, J.K [51] : route choices, separable link cost, decomposition method ; Bell,M.G.H [8] : route choice, link with capacity, separable link cost : Carey, M [17] : route choices, separable link cost, system-equilibrium reached by convex programming.

1

separable link cost : ) cij = cij (fij ) 8 link (i; j ) ~ symmetric non-separable link cost : ) cij = cij (f~) (f~ ow on the whole network) and @[email protected] (f ) = @[email protected] (f~) ~ asymmetric non-separable link cost : ) cij = cij (f~) and @[email protected] (f ) 6= @[email protected] (f~) ij

kl

ij

kl

kl

ij

kl

ij

In the same time, speci c approaches have been developed to solve this route-time allocation problem. The main ones are the following : Space-time network ; Variational Inequality ; Optimal Control ; Simulation ; Dunamic Programming (Ground Holding Problem).

2.2 Space-Time Network This technique consists of expanding the space network into a space-time network by creating N space-time nodes for each space node (where N is the number of duration periods) and then assigning the dynamic demands using classical static trac assignment methods. The associated number of links depends on the speed of the vehicle traveling on the network which is the weak point of this model [96, 31, 49, 56, 54, 97].

2.3 Variational Inequality The variational inequality problem is a general problem formulation that encompasses a plethora of mathematical problems, including, among others, nonlinear equations, optimization problems, complementarity problems, and xed point problem. It is now known that the simultaneous route and departure-time choice equilibrium problem can be formulated as a variational inequality problem [75]. Bernstein, D et al. [15, 39, 13] successfully used this method to solve the simultaneous route and departure-time choice equilibrium problem in networks with separable link costs. This approach is still theoretical and can be applied only to very small networks.

2.4 Optimal Control The optimal control approach has been tried for solving small instances of the general route-time assignment problem. Usually, the state variables are the link ows at time t and the associated control variables are the ows entering links during the same period of time.

Wie, B.W et al. [95] : "n" origins and \1" destination, route choices only, system-equilibrium, separable

link costs ; Wie, B.W et al. [94] : multi origin-destination, route choices, user-equilibrium, separable link costs ; Papageorgiou, M [77] : route choices, user-equilibrium and system-equilibrium, separable link costs ; Ran, B et al. [81] : multi origin-destination, route choices, user-equilibrium, separable link costs ; Friesz, T.L et al. [40] : \n" origin and \1" destination, route choices, user-equilibrium, system equilibrium ; Leblanc, L.J et al. [60] : optimal control modeling route-time choices user-equilibrium (modeling only, no algorithm) ; Janson, B.N et al. [53] : route-time choices, separable link costs ;

Optimal control techniques can be used only for small networks.

2.5 Simulation Due to the complexity induced by the dynamic route-time assignment in a general network (asymmetric nonseparable link costs), simulation has been investigated in a way to reach a solution close to the user-equilibrium or system-equilibrium. The strong point of the simulation techniques is their ability to take into account the maximum features of the real problem and can be applied on general networks. Cascetta, E et al [20, 19, 21, 18], Cantarella, G.E et al. [16] and Vythoulkas, P.C [92] applied simulation techniques in the following way : A set of possible slots of departure and a set of possible routes is assigned to each user. Each user randomly selects a couple of decision variables and the trac is simulated for the whole period considered (usually a day). A user speci c utility function is associated with each selected individual choice after the simulation. The choices of the day k are then based on the utilities of the N previous days using a kind of learning process (limited memory from the day k ? N to the day k). It can be shown that this stochastic process is a stationary Markov chain that reaches a xed point. This xed point is very stable but its distance from the user-equilibrium point or the system-equilibrium point can not be quanti ed and depends of the initial conditions. Others works based on the same methods can be found in the following references : [10, 64, 12, 63, 65, 66, 89, 23]

2.6 Dynamic Programming (Ground Holding Problem) The ground holding Problem was rst addressed and formalized by Odoni, A.R in [76]. Andreata, G and Romanin-Jacur, G [5] developed a dynamic programming framework to compute optimal delays for n ights arriving in the same time period at a single destination with a deterministic capacity scenario (a description of airport capacity scenarios is given in [45]). Terrab, M and Odoni, A.R [87] solved the multiple time periods, single arrival airport case as a minimum cost ow problem in a capacitated network, and as an assignment problem. They also developed a dynamic programming model to cope with stochastic capacities. However, this model was impractical for typical ow instances. The problem \single airport and random capacities" was solved by Richetta, G and Odoni, A.R [83, 82]. Vranas, P et al. [91, 90] addressed the multi-airport case with a deterministic airport capacity scenario (if there are not connected ights the multi-airport problem can be decomposed into several single airport problem). A good summary of those previous approaches is given in [4]. Bertsimas, D.J and Stock, S [14] proposed 0-1 programming models which incorporate airport capacities, sector capacities and connected ights. The model only assigns delays to aircraft and does not address route decisions. Finally Maugis, L [68] extended the previous model to take sector constraints in a more smoothly way.

2.7 Conclusion From the trac assignment point of view, the problem we have to solve could be summarized in the following way :

a eet of aircraft is considered for which one must nd the optimal route and the optimal time of departure in a way to reach a \system-equilibrium" in a network with asymmetric non-separable link costs

All the previous approaches are not able to manage the whole problem due to the complexity induced. In the following, a model is proposed and a method is developed in a way to give \very good" solutions for realistic instances of the whole problem.

3 A Simpli ed Model 3.1 Introduction

Congestion in the airspace is due to aircraft which have close positions in a four dimensional space (one time dimension and 3 space dimensions). It is then relevant to investigate ways to separate those aircraft in this 4 dimensional space by changing their slot of departure (time separation) or by changing their route (spatial separation) or both. Those changes must be done in a way that takes into account the objectives of the airlines : the moving of the slot of departure must be done in a limited domain (otherwise, for instance, some aircraft will be forced to take o at 2:00AM to reduce the congestion of the day but at this time there will be no passengers to carry);

Workload Morning peak

Evening peak

Time (hours) 4

5

6

7

8

9

10

11 12

13 14

15 16

17 18

19

20 21

22 23

0

1

Figure 1: Congestion variations

the new slot of departure must take into account the connexions between ight (as a mater of fact some

aircraft have to wait the arrival of some previous ights to take o (hub phenomenon); the possible routes must not generate too large additional distances. So, for each ight, a new pair (slot of departure, route) will be chosen from two discrete and nite sets : a set of possible slots of departure (around the original slot of departure) ; a set of routes which do not increase too much the total path length and are approved by the airline company the ight belongs to. After taking o the aircraft will follow its ight path and will generate workload in the dierent sectors encountered and also a congestion at the arrival airport when arriving. According to the controllers themselves, the workload induced in a control sector is a function of the three main following criteria : the con ict workload that results from the dierent actions of the controller to solve con icts. the coordination workload corresponds to the information exchanges between a controller and the controller in charge of the bordering sector or between a controller and the pilots when an aircraft crosses a sector boundary; the monitoring aims at checking the dierent trajectories of the aircraft in a sector and induces a workload. From the airport point of view the congestion is de ned by the number of aircraft unable to land due to the runways use (so it depends on the number of aircraft waiting in the sky). This airport congestion is dependent on the available capacity which changes over the day according to the weather. A static capacity has been attributed to each airport for our experiments, but the model is able to manage dynamic airport capacity, as well. If one sector is considered, one can observe variations of the associated workload during the day according to the trac in it. So, a typical pro le of the sector control workload is very chaotic with dierent peaks correlated with the peaks of trac in it (usually one peak in the morning and one peak in the evening see gure 1). It must be noticed that the peaks do not appear at the same time in the dierent sectors due to the ow propagations. According to our de nition of workload, only the information related with the route and the slot of departure are relevant from the ight plan. As a mater of fact, this information enables us to compute the time when an aircraft enters (or exits) a sector and when it is over a crossing of two airways (point of con ict generation). When examining the physical air trac network, we notice that airways are superpositions of several ight routes which have the same projection on the ground but dierent altitudes according to their azimuth (semi circular rule2 ). So an airway can be modeled by a bidirectional link which combines several individual aircraft routes. Then, without loss of generality, the real 3-dimensional transportation network will be modeled by a classical 2-dimensional network on a horizontal plan. We can now de ne our goals more precisely in the following way : 2 Aircraft with heading between 0 and 180 have to y at odd altitudes (in thousands of feet) and even altitude for headings between 180 and 360

A

13

D

18

7 2

17

12

19

6

1

16

20 21

5 11 3

B

8 15

C

22

24

10

4 23

14

9

E Beacon events

1

5

11

15

22

24 Sector events

B

C

E

Figure 2: Simpli ed Flight Plan one considers an air trac transportation network in a 2 dimensional space and a eet of aircraft with their associated route and slot of departure. For each ight a set of alternative routes and a set of possible slots of departure are de ned. One must nd \optimal" route and slot allocation for each aircraft in a way that signi cally reduces the peak of workload in the most congested sectors and in the most congested airports, during one day of trac. This \optimal" bi-allocation (route-time) must take into account the objectives of the airlines.

3.2 Mathematical formulation A pair of decision variable (i ; ri ) is associated with each ight in which i is the advance or the delay from the original slot of departure and ri is the new route. With this notation (0; r0 ) will be considered as the most preferred choice from the user point of view. Those two decision variables (i ,ri ) will be chosen from two nite-discrete sets : for the slots and R for the routes. More precisely the structures of and R are the following : = ?m; ?m + 1; ::::; ?1; 0; 1; :::; p ? 1; p R = r0 ; r1 ; r2 ; :::; rmax For which m ; p are respectively the maximum advance and the maximum delay permitted for a ight. Those limits can be dierent for each ight. The routes are ordered according to cost induced to the associated ight. So r0 is the best one and rmax the worst. To be able to compute the control workload (according to our de nition), the route of an aircraft will be de ned by the two following lists :

list of the over own beacons in the air network with the associated passing time : Lb = (b1 ; t1 ); (b2 ; t2 ); :::; (bi ; ti ); :::

list of crossed sectors with the associated in-out times : Ls = (S1 ; Tin1; Tout1); (S2 ; Tin2; Tout2); :::; (Sk ; Tink ; Toutk ); ::: An example of a simpli ed ight plan is given in gure 2. In this case an aircraft is joining airports 1 and 24. It over ies 5 beacons (5,11,15,22) and crosses 3 sectors. Given this information, it is now possible to compute the three components of the control workload (monitoring, coordinations, con icts) in each sector of the airspace. To reach this goal, one rst begins to compute the exiting

ow on each link of the network and the crossing ow over each sector. With the de nition of our simpli ed ight plan, it is easy to compute the discrete event when an aircraft exits a link of the network or when an aircraft enters or exits a sector. The ows are then computed by using a moving window in which an average of the number of events is computed (the temporal size of the window is a parameter of the model). So the ow fij (t) exiting a link (i; j ) at time t is computed in the following way :

fij (t) = 2:D1+ 1

s=X t+D nX =N s=t?D n=1

t ijn

where D is the half extension of the temporal window (in slot) and ijn (t) is the kroneker symbol which is equal to 1 if the ight n exits the link (i; j ) during the slot t (N is the whole number of ight). Given the ows on links, it is now possible to give a mathematical formulation of the control workload induced in a sector. As it has been previously said, workload in a sector Sk at time t can be roughly expressed by the summation of three terms : WSt = WmoS (t) + WcoS (t) + WcfS (t); k

k

k

k

WmoS (t) Monitoring workload ; WcoS (t) Coordination workload ; WcfS (t) Con ict workload. k

k

k

3.2.1 Monitoring workload

The monitoring workload is directly related to the number of aircraft in a sector and can be expressed in the following way : WmoS (t) = :n S (t) weight coecient and n(t) is the average number of aircraft in the sector Sk at time t. k

k

n S (t) = 2:D1+ 1 : k

l=X t+D l=t?D

nS (l) k

nS (t) is the number of aircraft in the sector Sk during slot t. nS (t) = nS (t ? 1) + NinS (t) ? NoutS (t) NinS (t) is the number of aircraft entering sector Sk during the slot t and NoutS (t) is the number of aircraft exiting sector Sk during the slot t. k

k

k

k

k

k

k

3.2.2 Coordination workload

During a slot t the number of coordinations induced in a sector Sk is equal to the number of aircraft entering or exiting the sector during the same time period. So, it is related with the ows cut by the sector boundaries : WcoS (t) = (FinS (t) + FoutS (t)) where is a weight coecient and FinS (t) and FoutS (t) are the entering and exiting sector ows associated with Sk at time t. k

k

k

FinS (t) = 2:D1+ 1 k

l=X t+D l=t?D

k

k

NinS (l); FoutS (t) = 2:D1+ 1 k

k

l=X t+D l=t?D

NoutS (l); k

NinS (t) and NoutS (l) are the number of discrete events \an aircraft enters sector Sk ", \an aircraft exits sector Sk " during the slot l. k

k

δAm

δAp

A

δm

δp

B

B

δm

Previous Flights (Arriving Intervals)

B

δp

C

C

δm

C

δp

D

D

D

τE δEp

δEm

E

Initial Interval

δEm

δEp

Connected Flight (Departure Interval)

Updated Interval

τE

Figure 3: Example of connexion

3.2.3 Con ict workload

The con ict workload induced in a sector Sk depends on the crossing ows at each node in the sector. X X X [h(ijl )fij (t)flj (t)] WcfS (t) = 12 j 2NO i 2 NO l 2 NO i 6= j l 6= i; l 6= j where : is a weight coecient ; NOS is the set of nodes belonging to the sector Sk ; NO is the entire set of nodes in the network ; ijl is the crossing angle of the links (i; j ) and (l; j ) ; k

Sk

k

h(ijl ) is a weighted function to take into account the fact that a con ict is harder to solve when the angle of crossing is small.

3.2.4 Formulation of the objective function

With the previous de nition of the control workload , when a set of decision variables has been chosen, it is now straightforward to compute the dynamic workload pro le over time for each sector and to memorize the associated maxima. We de ne our objective in the following way : \ one must try to reduce congestion in the most overloaded sectors but not necessarily simultaneously in all the sectors" ; this will spread the congestion over several sectors. So, we have :

obj = min

kX =K k=1

max W (t) t2T S

k

3.3 Connected ights constraint When a ight has to be connected to some arriving ights, its slot of departure must be later than the time of arrival of the previous ight and separated from it by a minimum amount of time ( ). So, when a slot of departure is changed, one must rst check that this new scheduling respects the connecting constraint. To reach this goal, the slot moving set of the connected ight must take into account the slot moving sets of the previous

ights as it can be seen on gure 3. On this gure, the ight \E" is connected to 4 previous ights :"A,B,C,D". To always be able to assign a slot to a ight, the slot moving sets must satisfy the following properties :

Maximum boundary

? pE > k2fA;B;C;D max g Pk + E

This checking must be done by using the longest possible route for each previous ights A; B; C; D. Minimum Boundary ? mE > k2fA;B;C;D max g tka + E where tka is the actual arrival time, determined by the slot moving and routes assigned to the ight A; B; C; D. So, before assigning a moving slot to a connected ight, one must rst check that the previous associated ights have already been assigned a slot and a route from their own route and moving slots sets. This will ensure that the random point generated in the state domain will always satisfy the connexion constraint.

3.4 Problem complexity Before investigating an optimization method, the associated complexity of our problem must be studied. The model previously developed is discrete and induces a high combinatoric in a huge space. As a mater of fact, if Rn ; n are the route set and the slot moving set associated with ight n, the number of possible choice (Cn ) for this ight is : Cn = jRn j:jn j where jS j denotes the cardinality of the set S . Then the number of possible decision variable combinaisons is given by :

jStatej =

nY =N n=1

Cn

where N is the total number of ights. If the number of choices is the same for each ight (jRn j = C te = r and n j = C te = ), then the cardinality of the state domain becomes : jStatej = (r:)N For instance, for 20000 ights with 10 route choices and 10 possible slot movings, : 10020000. Moreover, those decision variables are not independent due to the connexion induced by the control workload, the airport congestions and the connexion constraint, so, decomposition methods cannot be applied. It must be noticed that the objective function : is not continuous (then it is not convex) : may have several equivalent optima (the objective is then said multimodal) ; This problem is then a strong NP HARD problem with non-separable state variables.

3.5 Conclusion The microscopic model previously developed takes into account the major features of the real problem with the associated constraints and is able toaccommodate speci c individual preferences. This model has been completed by a mascroscopic description of the sector workload which induces a robust computation of sector congestion. So, one must solve a multimodal combinatorial optimization problem in a huge space with non-separable decision variables. Currently, only stochastic optimization is well adapted to address this kind of problem.

4 Stochastic Optimization

4.1 Introduction

Stochastic optimization techniques usually addressed problem with strong complexity for which the objective function has no speci c feature such as convexity, continuity etc.... The main characteristic of those stochastic optimization techniques is to randomly move in the state domain in a way to improve the objective criterium. At the beginning, local search algorithms randomly moved in the state domain toward the closest optimum without being able to avoid local optima. As a mater of fact, this technique consist in rst selecting a random point in the state domain, second to nd a neighbor point and to compare the associated objectives : if the neighbor is better it begins the currant point else one must nd another neightbor point till the objective is improved. Afterward, more sophisticated and more powerful methods have been developed in a way that avoid local optima and converge toward global optima (Simulated Annealing and Genetic Algorithm) and which can identi ed several global optima (Genetic Algorithm). Those algorithm are mainly used in combinatorial optimization we now describe.

4.2 Combinatorial Optimization Combinatorial optimization problems represent a class of computational problems where one seeks the one \best" or \optimal" solution from among a large space of solutions [79]. Generally, the optimal solution is speci ed in terms of minimizing some given cost function subject to a set of constraints. Formally, a combinatorial optimization problem can be stated as a tuple < S; f > , where S , the state space, is the ( nite) set of possible solutions and f is a cost function of the form : f :S ??>R mapping each possible solution onto the reals. The optimal solution to any such combinatorial optimization problem is a solution, sopt 2 S , such that : f (sopt ) f (s); 8s 2 S sopt is called a globally optimal solution. In some cases, it can be more natural to formulate the optimal solution as that solution which maximizes, rather than minimizes, some quantity. Such cases can be readily handled, without loss of generality, within the formalism here by simply noting that maximization is equivalent to minimization with a reversal of sign in the cost function. Consider some solution, s 2 S , and the set of all solutions, Ss , which are in some sense \close" to s. The set Ss is called the neighborhood of s . We can de ne the neighborhood structure of S as a mapping N : S ! 2S which formalizes the notion of \closeness" for a particular optimization problem. If < S; f > is an optimization problem and N is a neighborhood structure, then s^ 2 S is called a local minimum, with respect to N provided that : f (^s) f (s); 8s 2 N The relationship between globally and locally optimal solutions can be visualized in the 2-dimensional case as in gure 4. One class of approximation algorithms used in solving optimization problems is known as local search. Local search algorithms operate by stepwise re nement based on a neighborhood structure. Given an initial solution, usually chosen at random, a local search algorithm evaluates the cost function on the initial solution. Next, subsequent candidate solutions are generated in the neighborhood of the currently optimal solution. If the cost of the newly generated candidate solution is less than the currently optimal solution, it becomes the currently optimal solution. This process is iterated until a satisfactory locally optimal solution is found, no improvement has been made for some number of previous iterations, or until some maximum number of iterations has been reached. By de nition, then, local search algorithms terminate in locally optimal solutions. Unless the neighborhood structure is de ned such that all locally optimal solutions are also globally optimal, then it is clear that local search cannot be guaranteed to result in a globally optimal solution. Furthermore, there is no easy answer to the question of how close the locally optimal solution found by any particular local search algorithm is to a globally optimal solution. In general, the quality of the solution found by local search is highly dependent on both the initial solution selected and the method of generating candidate solutions from the neighborhood of any particular solution. Those drawbacks can be avoided by using more sophisticated stochastic optimization techniques such as : genetic algorithms.

f(x)

Local optima

Global optimum x

Figure 4: Global vs. Local Minima

4.3 Classical Genetic Algorithm Genetic algorithms (GAs) [46, 52, 37, 72, 58, 36, 84] are problem solving systems based on principles of evolution and heredity. A GA maintains a population of individuals,P (t) = x1 ; x2 ; :::; xn at iteration t. Each individual represents a potential solution to the problem at hand and is implemented as some (possibly complex) data structure S . Each solution xi is evaluated to give some measure of tness. Then a new population at iteration t + 1 is formed by selecting the tter individuals. Some members of the new population undergo transformation by means of genetic operators to form new solutions. There are unary transformations mi (mutation type), which create new individuals by a small change of a single individual and higher order transformations ci (crossover type), which create new individuals by combining parts from several (two or more) individuals. For example, if parents are represented by a ve-dimensional vector (a1 ; a2 ; a3 ; a4 ; a5 ) and (b1 ; b2 ; b3; b4 ; b5 ), then a slicing crossover of chromosomes after the second gene produces ospring (a1 ; a2 ; b3 ; b4; b5 ) and (b1 ; b2 ; a3 ; a4 ; a5 ). The control parameters for genetic operators (probability of crossover and mutation) need to be carefully selected to provide better performance. The intuition behind the crossover operation is information exchange between dierent potential solutions. After some number of generations the program converges - the best individual hopefully represents the optimum solution. In [46], Goldberg presents a good summary of many recent GA applications in biology, computer science, engineering, operations research, physical sciences, and social sciences. Genetic algorithms use a vocabulary borrowed from natural genetics. We would talk about genes (or bits), individuals (chromosomes or gene strings), and population (of individuals). Populations evolve through generations. The execution steps of genetic algorithms are the following(see gure 5) : 1. Initialize population and evaluate tness : To initialize a population, we need rst to decide the number of genes for each individual and the total number of chromosomes (popsize) in the initial population. All the information contained in a state point must be summarized in the chromosome, this results from the coding which is one of the success keys of genetic algorithms. 2. Reproduction (Selection) : Reproduction is the selection of individuals with respect to the probability distribution based on the the tness values. Fitter individuals have better chances of getting selected for reproduction [72]. Dierent kinds of selection processes can be applied for a speci c problem; one of the most simple is the \Roulette Wheel Selection Process" : Each chromosome has a certain number of slots proportional to its tness value. The selection process is based on spinning the wheel many times; each time we select a single chromosome for a new population. Obviously, some chromosomes will be selected more than once. This is in accordance with genetic inheritance: the best chromosomes get more copies, the average stay even, and the worst die o. 3. Crossover : Now we are ready to apply the rst recombination operator, crossover, to selected individuals. The probability of crossover,pc , gives us the expected number pc:popsize of chromosomes which should undergo

population(k) selection P1

P2

crossover C1

C2

mutation

mutation

C1

C2

fitness evaluation

fitness evaluation

population(k+1)

Figure 5: Classical Genetic Algorithm the crossover operation. After having selected two parents (according to their tness), we generate a random number r between 0 and 1; if r < pc , then the parents undergo a crossover. We then generate a random number pos from the range of (1::m ? 1), where m is the total number of genes in a chromosome. The number pos indicates the position of the crossing point. The coupled chromosomes exchange genes at the crossing point as described earlier. 4. Mutation : Chromosomes which have not undergone a crossover may be randomly mutated on the gene by gene basis. To perform a mutation, a gene is randomly selected and modi ed by a problem dependent process in a way to produce new genes. The probability of mutation,pm , gives us the expected number of mutated gens pm :m:popsize. 5. Convergence : Following reproduction, crossover, and mutation, the new population is ready for its next generation. The rest of the evolutions are simply cyclic repetitions of the above steps until the system reaches a predetermined number of generations or converges (i.e., no improvement in the overall tness of the population). This classical scheme can be improved by using simulated annealing in the crossover operators.

4.4 Genetic Algorithm with Simulated Annealing Crossover This Simulated Annealing Crossover has been developed recently [28] and has been used in dierent applications with success [27, 29, 32]. It greatly improves the convergence speed of the genetic algorithm for all applications, so it is generic and not limited to a problem dependent improvement. This operator is based on another stochastic optimization technique called : Simulated Annealing.

4.4.1 Simulated Annealing

Simulated annealing is based on an analogy to the thermal process of metallurgical annealing from statistical mechanics; an examination of this process will help to better understand the simulated annealing process [57]. In metallurgical annealing, the goal is to bring a solid, in contact with heat bath, to a con guration with minimal energy. Distortions in the crystal lattice structure of the metal result in higher energy states. The central idea behind metallurigal annealing is that disruptions in the crystal lattice structure of the metal are

removed by thermal agitation at very high temperatures and that further disruptions can be prevented from forming if the metal is cooled slowly enough. Metropolis, et al. [71] succeeded in simulating the evolution of states a metal undergoes as it reaches thermal equilibrium, at constant temperature, using Monte Carlo methods. The technique, which has come to be known as the Metropolis Algorithm (MA), is as follows. Given a current state, generate a new state by randomly perturbing the current state (thus, simulating the random motion of particles). If the new state results in lower overall energy then retain the change as the current state; if the result is a higher energy state, then retain the new state with probability P , given by :

p = e? where E = Enew ? Eold is the energy change, T is the temperature of the heat bath and kb is the Boltzmann constant. For a system in thermal equilibrium at temperature T , the probability of nding the system in a state, s, with energy Es is given by E kb :T

PT (s) = Z (1T ) e?

Es kb :T

where Z (T ) , the so-called partition function, is given by

Z (T ) =

X

i

e?

Ei kb :T

where the summation is over all possible states of the system. PT (s) is known as the Boltzmann distribution. The similarities between the local search technique presented in the combinatorial optimization section and MA suggest the possibility of using MA to escape local minima. If < S; f > is an instance of an optimization problem and s1 ; s2 2 S are solutions with associated costs f (s1 ) and f (s2 ), then the acceptance criterion by which s2 is retained is determined by the acceptance probability : 1 if f (s2 ) f (s1 ) PT (s1 ! s2 ) = e ( 1 )? ( 2 ) otherwise f s

f s

T

where T 2 [0; 1] is a monotonically decreasing control parameter. In addition to the acceptance criterion and probability, full speci cation of the simulated annealing algorithm also requires a sequence fT0 = 1; T1 ; :::; Tk ; :::; Tfinal > 0g, called the annealing schedule and ,uk a further control parameter stipulating the number of solution candidates, or transitions, generated at T = Tk . So, given an optimization problem, < S; f >, then under reasonable restrictions on the neighborhood structure and a suciently large number of transitions at a xed value of the control parameter T , it has been shown [2] that the simulated annealing algorithm, using the acceptance probability, will nd a solution, s 2 S , with probability : () PT (s) = N 1(T ) e? 0 where N0 (T ) is given by : f s T

N)(T ) =

X

i2S

e?

()

f i T

PT (s) is known as the stationary distribution of the solution s , the equivalent of the Boltzmann distribution for physical systems. Denote the set of all globally optimal solutions to < S; f > by Sopt . It has been shown [61], for s 2 Sopt , that : lim P (s) = 1 (s) T !0 T

Sopt (S ) where (S ) (s) is de ned in the following way : Let A and A0 2 A be two sets. Then (A0 ) (a) = 1 i a 2 A0 , opt

opt

and 0 otherwise. The above result is extremely valuable, because it guarantees that the simulated annealing algorithm will asymptotically converge to Sopt , provided that the stationary distribution is reached for each value of the control parameter, T . However, asymptotic convergence is only possible given an in nite number of transitions [6]. Clearly, an in nite time convergence is not of practical interest; although it leads to the question of the eciency of approximating the asymptotic behavior. At least two important results in this regard have been

P1

P2

CLASSICAL CROSSOVER

P1

P2

C1

Pc

C2

TOURNAMENT

C1

C2

Figure 6: Simulated Annealing Crossover reported. Aarts and van Laarrhoven [1] have shown that the stationary distribution is approximated arbitrarily closely only under the condition that the number of transitions is at least quadratic in jS j , the cardinality of the solution space. Attacking the problem dierently, Mitra, et al. [73] have shown that for many problems, arbitrarily close approximation to the asymptotic behavior requires a number of transitions which is larger than jS j , leading to exponential-time execution. That is, enumerating all possible solutions requires less time than needed to perform the simulated annealing procedure until an arbitrarily close approximate solution is reached. We now return to consideration of the annealing schedule and uk , the number of transitions attempted at each value of T . By suitably choosing these two parameters (T and Uk ), polynomial-time execution is obtainable [1], at the expense of not being able to guarantee the degree of divergence (the distance) between the asymptotically optimal and the approximated solution. In general, the annealing schedule is chosen according to the following guidelines : 1. the initial value, T0 , is such that PT0 (s) = 1; 2. the decrement rule, Tk+1 = Tk ; k = 1; 2; :: where < 1, and typically 0:9 < < 0:99 ; 3. the nal value, or stopping criterion, Tfinal , is that value at which suciently many successive transitions have not resulted in a signi cant improvement in reducing the cost function The choice of initial value is motivated by noting that such a value results in virtually all attempted transitions being accepted, thus approximating P1 . It must be noticed that the larger the dierence between Tk and Tk+1 the more transitions, uk , that will be required to approximate the stationary distribution at Tk+1 [57, 2].

4.5 Simulated Annealing Crossover The way to introduce Simulated Annealing in the crossover operator is straightforward. Two parents (P1 ; P2 ) are selected in the population according to their tness and randomly crossed using the classical crossover operator. If the crossover happened, two children are generated (C1 ; C2 ). The rst child (C1 ) is then randomly associated with one parent and the second one (C2 ) with the remaining parent. For each created couple, the \winner" is then selected according to a simulated annealing process in the following way (see gure6) : if f (C ) > f (P ) ) C 0 = C ( )? ( ) if f (C ) < f (P ) ) C 0 = C with the probability P (C ! C 0 ) = e ( ) where T (gen) is the control parameter depending of the current generation gen. f C f P T gen

4.6 Conclusion The strong point of these stochastic optimization methods is their ability to give \very good" solution(s) to large instances of real problems. Usually the given solutions are very robust and very close to the optimum when it is known. They can manage very speci c constraints and do not need any properties concerning the objective function.

SECTOR CONGESTIONS

Traffic Simulator

Genetic Algorithm

FPL

CHROMOSOME

Flight Plans Generator

Figure 7: GA using simulated tness R1

n rk

1 r3

∆n

N r6

δnj

δ12 ∆1

RN

Rn

δN3 ∆N

∆n

: Set of slots for the flight n

Rn : Set of routes for the flight n

Figure 8: Structure of the Chromosome

5 Application to our problem 5.1 Introduction

As previously noted, Genetic Algorithms need a data coding, speci c recombination operators and a tness evaluation to be able to select best individuals. The way this speci c genetic algorithm works is the following. A set of ight plans is generated from each chromosome candidate and the whole associated day of trac is generated. Sector congestion and airport congestion are registered and the associated tness is computed (see gure 7). This is an appropriate method to use because the more realistic the simulation of the day of trac is the more relevant are the results. The problem speci c features of the Genetic Algorithm are now described.

5.2 Data Coding This step consists of converting each point of the state domain into a chromosome used by the genetic algorithm. In our problem, the state variables (which contain all the information needed to compute the sector workload) consist of the set of ight plans. The possible new path and new slot moving have been supposed to be chosen in two discrete- nite sets associated with each ight. In this case a straightforward coding has been used in the sense that each chromosome is built as a matrix which gather together the new slot moving (for the time of departure) and the new route number (for the ight path) (see gure 8). With this coding, a population of individuals can be created by randomly choosing a new slot moving number and a new route number from individual sets associated with each ight.

SLICING CROSSOVER

Figure 9: Crossover operator

5.3 Fitness Evaluation The previous description of the chromosome and the associated recombination operators can now be used to randomly generate some points in the state domain. To apply the selection operator, a tness must be associated with each chromosome in a way to evaluate the quality of each individual according to the optimization criterion. In our problem, one rst considers the default ight plans and computes the maximum workload in each sector in the airspace.

MaxW (ref ) =

kX =K k=1

max W (t; ref ) t2T S

k

Afterward, for each generated chromosome, maximum sector workloads are computed the same way :

MaxW (chrom) =

kX =K k=1

max W (t; chrom) t2T S

k

The tness is then de ned by the ratio of the two previous expressions. MaxW (ref ) fitness(chrom) = MaxW (chrom) So, when fitness(chrom) > 1, it means that the induced congestion is lower than the reference one.

5.4 Crossover Operator The crossover operator has been implemented as a slicing one on the rst attempt. So, to make a crossover, one rst select two parents and randomly chooses an intergene site from which the two induced ospring are exchanged (see gure 9).

5.5 Mutation Operator In our problem, the mutation is used to generate new slot moving or new route from the two discrete- nite sets associated with each ight. Then, a new couple (slot moving,route number) is generated from each individual associated set of slot moving and route for all the ights (see gure 9).

n

Rn

∆n

MUTATION

Figure 10: Mutation operator

6 Application to a simpli ed network 6.1 Introduction

A simpli ed network has been used to evaluate the performances of the algorithm. This example tries to reproduce the trac density encountered usually in airspace. The evaluations have been done with ight plans generated in a two dimensional network but the model could be applied to the real ight plans in three dimensional airspace, due to the fact that only ight plan events are used to compute the objective function. In this example, only sector congestion has been taken into account but as it has been mwntioned before, airport congestion could have been introduced by using airports with smaller capacitiey. The features of the network are the followings : Number of nodes : 100 Number of links : 544 Number of Origin-Destination pairs : 30 Number of Sectors : 15 Number of Airports : 24 Geographical size of the network : 1000NM x 1000NM A physical description is given in gure 11 with the sectoring and the associated sector number. Each airport is symbolized with a larger node. A hundred ights have been randomly generated for each Origin-Destination pair, so, a total of 3000 ight plans is used for this example (the speed of the aircraft has been xed at 500 kts). Those ights have been spread over 6 hours and the number of slots has been xed at 180. This induces a slot duration of 2 minutes. For each Origin-Destination, a set of possible dierent routes is randomly built in a way that never induces an extradistance greater than 30 percent of the original distance with a maximum of 10 routes. For each ight the maximum slot moving has been xed to 15 in the future and 15 in the past from the original slot of departure. Two kinds of experiments have been done to compare the gain given by the delay allocation only (the departure slot can only be moved in the future (not in the past) and the original route is used) and the time-route allocation (dierent routes can be used and the moving of the slot of departure can be done in the past or in the future). The rst is a kind of Ground Holding Policy(GHP) and the second one is a Full Planning Policy(FPP). The original plannings induce a control workload in the airspace with a maximum for each sector at a given slot. The workload weight coecients used for these experiment are : = 1:0; = 1:0; = 1:0. The half extension of the temporal window used to compute the ows from the discrete events has been xed at D = 2 and induces a window duration of 10 minutes. The maximum control workload registered in each sector over the six hour duration is given in gure 12 (with this data set). On this gure, sector indices are indicated on the horizontal axis and the associated maximum sector workload on the vertical axis. One can notice that sector 9 is the most overloaded. For this problem, with the de nition of the tness, the genetic algorithm is expected to reduce those peaks of congestion for the most congested sectors.

6

12

8 15

3

13 14

4 9

10 11

2 7

5 1

Figure 11: Structure of the test network Maximum sector workloads over the day 2000

1500

1000

500

1

2

3

4

5

6

7

8

9

10 11 12 13 14 15

Figure 12: Max workload distributions over sectors

Fitness evolution (ALEAT100) 2.5

Bestfit

2

1.5

1

5

10

15

20

5

10

15

20

25

30

35

40

45

50

25 30 generations

35

40

45

50

2.5

Avgfit

2

1.5

1

Figure 13: Best and average tness evolution over generations

6.2 Results The results presented in this section are obtained from runs on a SPARC5 workstation. The large size of the chromosomes (6000 bytes) limits the maximum size of the population used by the genetic algorithm at 90 due to memory limitations (a population of 50 individuals has been used for all the experiments). It must be noticed that this small population size will reduce the possibilities of genetic algorithms, and better results could be reached on a machine with more memory (even if the following results are already very good). The others parameters (probability of crossover, probability of mutation, number of generations) have been xed by experimentation at Pc = 0:6, Pm = 0:2 and NB gen=100. The technical features of the genetic algorithm were the following : Selection : Stochastic remainder without replacement ; Crossover : Crossover with a Simulated Annealing Tournament ; Scaling : Power Law Scaling ; Elitist Strategy Applied ; No scharing. Given all these parameters, each run lasts about 15 minutes on a SPARC5 workstation. This running time is not critical because the problem solved has not any real time constraints (the improved time-route planning could be applied for several months). The associated tness evolution is given in gure 13. On this gure, the two kinds of policies (GHP and FPP) can be identi ed by the dierent lines used (plain for FPP and dashed for GHP). As expected, the FPP policy gives better results than the GHP policy with respectively a best tness of 1.79(FPP) and 1.45(GHP). These results are relevant due to the fact that the FPP policy has more degrees of freedom. It means that original \congestion" in the sense given by the optimization criterion, could be respectively divided by 1.79 and 1.45 with the new ight plans. Even with the small population size used, the results given by the genetic algorithm are very encouraging. The associated max workload distributions over sectors is given on gure 14. On this gure, it can be noticed that the max workload on the most overloaded sector (number 9) has been 1841 divided by 1.673 ( 1100 ) for the GHP policy and by 2.77 ( 1841 ) for the FPP policy. A more detailed description of 663 the sector control workload evolutions is given in gures 15,16,17 (on these diagrams workload has been divided by 1000 for presentation). As expected, the workload is spread around the peak as in a smoothing process. The previous curves make us notice that the algorithm only aects the most congested sectors without modifying the less loaded ones. This is due to the de nition of the criterion which takes into account only the most loaded sectors. So, when the recombination operator is applied, if they aect aircraft involved in the underloaded sectors (or uncongested airports) the induced tness will not change a lot and this makes the crossover and mutation irrelivant. Then, those operators have to be changed in a way that aect only aircraft involved in the most loaded sectors (or in the most congested airports) during the peak periods. This is the subject of the following section.

Maximum sector workloads over the day 2000

1500

1000

500

1

2

3

4

5

6

7

8

9

10 11 12 13 14 15

Figure 14: Max workload over sectors for the two associated policies (-. GHP ... FPP)

W1

2 1

W2

0 0 2

W3

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

1 0 0 2

W4

40

1 0 0 2

1 0 0 2

W5

20

1 0 0

Figure 15: Sector control workload evolution over the day (-. GHP ... FPP) W6

2 1

W7

0 0 2

W8

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

1 0 0 2

W9

40

1 0 0 2

1 0 0 2

W10

20

1 0 0

Figure 16: Sector control workload evolution over the day (-. GHP ... FPP)

W11

2 1

W12

0 0 2

W13

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

1 0 0 2

W14

40

1 0 0 2

1 0 0 2

W15

20

1 0 0

Figure 17: Sector control workload evolution over the day (-. GHP ... FPP)

7 Improvement of the recombination operators 7.1 Introduction

To be able to recognize the aircraft involved in the biggest sector congestion peak or in the biggest airport congestion peak, new information must be added to the chromosome which indicates for each gene, the maximum level of sector congestion encountered during a ight and the level of congestion at the destination airport when the aircraft land. To reach this goal, the biggest peak of congestion is determined over all the sectors and airports, and a number of congestion levels is de ned respectively for the airports and the sectors. With these two parameters, and for each slot, a level of congestion can be associated with each sector and each airport. So, after each evaluation, the associated list of sectors of each individual simpli ed ight plan is read in a way to identify the maximum level of congestion encountered during the ight and at the destination airport. Those congestion levels will be used to improve the recombination operators.

7.2 New Crossover Operator The successive steps of this new crossover operator are the following :

two parents are rst selected according to their tness ; the summation of the congestion levels (sector and airport) is computed for each ight in bothp parents. For p p p

a ight n, total congestion level in the parent p will be noted Wn =WSn + WAn (with WSn : maximum level of sector congestion encountered by the ight n in the parent p, WApn : airport congestion level encountered by the ight n when arriving at destination (in the parent p)) ; ; an order relationship is then constructed with the total congestion level in the following way :

{ ight planing n in parent 1 is said to be \much better" than ight planing n in parent 2 if Wn < :Wn ; where 2 [0:7; 0:95]; 1

2

Fitness evolution (ALEAT100) 2.5

Bestfit

2

1.5

1

5

10

15

20

5

10

15

20

25

30

35

40

45

50

25 30 generations

35

40

45

50

2.5

Avgfit

2

1.5

1

Figure 18: New Best and average tness evalution over generations

{ ight planing n in parent 2 is said to be \much better" than ight planing n in parent 1 if Wn < :Wn ; { ight planing n in parent 1 and in parent 2 are said to be \equivalent" if none of the previous relations 2

1

matches; if a ight planning \is much better" in the rst parent than in the second then it is copied in the second ; if a ight planning \is much better" in the second parent than in the rst then it is copied in the second ; if the two ight plannings \are equivalent" they are randomly exchanged with a constant probability (0.5) ; This is done for all the ights in the chromosome.

7.3 New Mutation Operator As already noted, this new operator must only aect the ights involved in the highest peaks of congestion. This new operator works in the following way :

two threshold congestion levels are randomly chosen : one for the sector congestion (ThS ) and one for the

airport congestion (ThA) ; then for each ight n in the chromosome the following are applied : if (WAn > ThA) or(WSn > ThS ) then the associated ight plan is modi ed by randomly choosing either a new route or a new slot moving (but not both) in the associated data set (Rn ; n ) ; else the ight planing is unchanged;

7.4 New Results The problem described in section 6 has been resolved using this new genetic algorithm on the same machine (SPARC5) with the same parameters (population size, probability of crossing etc...). The same policies have been compared (GHP and FPP) one the same random trac sample. These new operators greatly improve the tness evolution as it can be seen in gure 18. The associated new best tnesses are 2.32 for the FPP and 1.68 for the GHP, as compared with the results given by the classical recombination operators (FPP : 1.79 and GHP : 1.45). The associated max workload distributions over sectors are given on gure 19. The improvement on the most overloaded sector (number 9) is now 2.73 with the GHP ( 1841 ) and 6.71 with 672 ). The associated time evolution of the sector workload is given in gures 20,21,22. the FPP( 1841 274

Maximum sector workloads over the day 2000

1500

1000

500

1

2

3

4

5

6

7

8

9

10 11 12 13 14 15

Figure 19: New max workload over sectors for the two associated policies (-. GHP ... FPP)

W1

2 1

W2

0 0 2

W3

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

1 0 0 2

W4

40

1 0 0 2

1 0 0 2

W5

20

1 0 0

Figure 20: Sector control workload evolution over the day (-. GHP ... FPP) W6

2 1

W7

0 0 2

W8

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

1 0 0 2

W9

40

1 0 0 2

1 0 0 2

W10

20

1 0 0

Figure 21: Sector control workload evolution over the day (-. GHP ... FPP)

W11

2 1

W12

0 0 2

W13

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

1 0 0 2

W14

40

1 0 0 2

1 0 0 2

W15

20

1 0 0

Figure 22: Sector control workload evolution over the day (-. GHP ... FPP)

8 Conclusion This dicult problem has been handled well by the proposed genetic algorithm even with a small population (50 individuals) and a small number of generations (100). Even better results could be expected with larger populations. A mathematical model has been developed to formulate this bi-allocation (route and slot of departure) problem in a way that could be adapted to the real ight plans in three dimensional airspace. So, the next step, is to try this optimization technique on a real trac sample. The strength of this model is its ability to manage the constraint of the airline companies in a microscopic way by using individual sets of decision variables associated with each ight. Using those individual decision variables, macroscopic models of the sector control workload and of airport congestion have been developed by computing the associated ows with a time moving window on the strategic points (crossings, sector boundaries, airports). This reduces the sensitivity of the \optimum" given by the genetic algorithm with the proposed solution. This model can also take into account the connexion between

ights (hub phenomenon) by changing the decisions variable sets of the following ight. This model can also use dynamic capacities for both airports and sectors. A fast simulation process has been used to compute the associated sector workload and the airport congestion to improve the performance of the genetic algorithms. A realistic example has been treated and two policies have been tested : GHP (delay decision only) and FPP (route and slot of departure decision). As expected the second one gives much better results due to the larger state domain induced by the two extra degrees of freedom of the FPP policy. Two kinds of genetic algorithms have been tested : the rst one has a classical crossover(slicing) and mutation( random gene) operators and the second one focuses its recombination actions on the most congested sectors by assigning two levels of congestion for each ight(one for the most congested sector encountered during the ight and one for the congestion of the destination airport when the ight is arriving). This new information greatly improves the quality of the given solutions. It must be noticed, that only aircraft involved in the most congested parts of the airspace (sector or airport) receive a ight plan modi cation. A FPL change cost could be added to classify aircraft and identi ed the one which must be modi ed rst.

References [1] A.H.L Aarts and P.J.M Van Laarhoven. Statistical cooling : A general approach to combinatorial optimization. Philips Journal of Research, 40:193{226, 1985. [2] E.H.L Aarts and J Korst. Simulated Annealing and Boltzmann Machines : A Stochastic Approach to Combinatorial Optimization and Neural Computing. John Willey and Sons, New York, 1989. [3] A.S Alfa. Departure rate and route assignment of commuter trac during peak period. Transportation Research, 23B(5):337{344, 1989. [4] G Andreatta, A.R Odoni, and O Richetta. Models for the ground holding problem. In L Bianco and A.R Odoni, editors, Large Scale Computation and Information Processing in Air Trac Control, Transportation Analysis, pages 125{168. Springer-Verlag, 1993. [5] G Andreatta and G Romanin-Jacur. Aircraft ow management under congestion. Transportation Science, 21(4), 1987. [6] S Anily and A Federgruen. Simulated annealing methods with general acceptance probabilities. Advances in Applied Probability, 24:657{667, 1987. [7] R Arnott, A DePalma, and R Lindsey. Departure time and route choice for the morning commute. Transportation Research, 24B(3):209{228, 1990. [8] M.G.H Bell. Stochastic user equilibrium assignment in networks with queues. Transportation Research, 29B(2):125{137, 1995. [9] M Ben-Akiva. Dynamic network equilibrium research. Transportation Research A, 19A(5):429{431, 1985. [10] M Ben-Akiva, A DePalma, and P Kanaroglou. Dynamic model of peak period trac congestion with elastic arrival rates. Transportation Science, 20(2):164{181, 1986. [11] M Ben-Akiva, A DePalma, and I Kaysi. Dynamic network models and driver information systems. Transportation Research, 25A(5):251{266, 1991. [12] M Ben-Akiva, H.N Koutsopoulos, and A Mukundan. A dynamic trac model system for atms/atis operation. IVHS Journal, 2(1):1{9, 1994. [13] D Bernstein, T.L Friesz, R.L Tobin, and B.W Wie. A variational control formulation of the simultaneous route and departure-time choice equilibrium problem. In Transportation and Trac Theory, 1993. [14] D.J Bertsimas and S Stock. The air trac ow management problem with en-route capacities. Technical report, A.P Sloan School of Management. M.I.T, 1994. [15] Wie B.W., R.L Tobin, T.L Friesz, and D Bernstein. A discrete time, nested cost operator approach to the dynamic network user equilibrium problem. Transportation Science, 29(1):79{92, 1995. [16] G.E Cantarella and E Cascetta. Dynamic processes and equilibrium in transportation networks : Toward a unifying theory. Transportation Science, 29(4):305{328, 1995. [17] M Carey. Optimal time-varying ows on congested networks. Operations Research, 35(1):58{69, 1987. [18] E Cascetta. Static and dynamic models of stochastic assignment to transportation networks. In Flow Control of Congested Networks, 1987. [19] E Cascetta. A stochastic process approach to the analysis of temporal dynamics in transportation networks. Transportation Research, 23B(1):1{17, 1989. [20] E Cascetta and G.E Cantarella. A day-to-day and within-day dynamic stochastic assignment model. Transportation Research, 25A(5):277{291, 1991. [21] E Cascetta and G.E Cantarella. Modeling dynamics in transportation networks : State of the art and future developments. In Simulation Practice and Theory 1. SIMPRA, 1993.

[22] G.L Chang and H.S Mahmassani. Travel time prediction and departure time ajustment behavior dynamics in a congested trac system. Transportation Research, 22B(3):217{232, 1988. [23] G.L Chang and H.S Mahmassani. Travel time prediction and departure time ajustment behavior dynamics in a congested trac system. Transportation Research, 22B(3):217{232, 1988. [24] S Dafermos. Trac equilibrium an variational inequalities. Transportation Science, 14(1):42{54, 1980. [25] S Dafermos. Relaxation algorithms for the general asymmetric trac equilibrium problem. Transportation Science, 16(2):231{239, 1982. [26] S Dafermos and F.T Sparrow. The trac assignment problem for a general network. Journal of Research of the National Bureau of Standards, 73B:91{118, 1969. [27] D Delahaye, J.M Alliot, M Schoenauer, and J.L Farges. Genetic algorithms for air trac assignment. In Proceedings of the European Conference on Arti cial Intelligence. ECAI, 1994. [28] D Delahaye, J.M Alliot, M Schoenauer, and J.L Farges. Genetic algorithms for partitioning airspace. In Proceedings of the Tenth IEEE Conference on Arti cial Intelligence for Application. CAIA, 1994. [29] D Delahaye, J.M Alliot, M Schoenauer, and J.L Farges. Genetic algorithms automatic regrouping of air trac sectors. In Proceedings of the Fourth International Conference on Evolutionary Programming. Natural Selection inc., 1995. [30] A DePalma, P Hansen, and M Labbe. Commuters'paths with penalties for early or late arrival time. Transportation Science, 24(4):276{286, 1990. [31] O Drissi-Kaitouni and A Hameda-Benchekroun. A dynamic trac assignment model and a solution algorithm. Transportation Science, 26(2):119{128, 1992. [32] N Durand and Alliot. Automatic aircraft con ict resolution using genetic algorithms. In Proceedings of the 11th Annual ACM Conference on Applied Computing, Philadelphia, 1996. (ACM/SAC. [33] J.E Fernandez and T.L Friesz. Equilibrium predictions in transportation markets : the state of the ar. Transportation Research B, 17B(2):155{172, 1983. [34] C Fisk and S Nguyen. Solution algorithms for network equilibrium models with asymmetric user costs. Transportation Science, 16(3):361{381, 1982. [35] M Florian and H Spiess. The convergence of diagonalization algorithms for asymmetric network equilibrium problems. Transportation Research B, 16B(6):477{483, 1982. [36] D.B Fogel. Evolutionary Computation. Toward a new Philosophy of Machine Intelligence. IEEE press, 1994. [37] L.J Fogel, A.J Owens, and M.J Walsh. Arti cial Intelligence Through Simulated Evolution. Wiley and sons. NY, 1966. [38] T.L Friesz. Transportation network equilibrium, design and aggregation : Key developements and research opportunities. Transportation Research, 19A(5):413{427, 1985. [39] T.L Friesz, D Bernstein, T.E Smith, and B.W Wie. A variational inequality formulation of the dynamic network user equilibrium problem. Operations Research, 41(1):179{191, 1993. [40] T.L Friesz, J Luque, R.L Tobin, and B.W Wie. Dynamic network trac assignment considered as a continuous time optimal control problem. Operation Research, 37(6):893{901, 1989. [41] N.H Gartner. Optimal trac assignment with elastic demands : A review. part i : Analysis framework. Transportation Science, 14(2):174{191, 1980. [42] N.H Gartner. Optimal trac assignment with elastic demands : A review. part ii : Algorithmic approaches. Transportation Science, 14(2):191{208, 1980. [43] N.H Gartner, S.B Gershwin, J.D.C Little, and P Ross. Pilot study of computer based urban trac management. Transportation Research B, 14B(1):203{217, 1980.

[44] M.O Ghali and M.J Smith. A model for the dynamic system optimum trac assignment. Transportation Research, 29B(3):155{170, 1995. [45] E.P Gilbo. Airport capacity : Representation, estimation, optimization. IEEE Transaction on Control Systems Technology, 1(3):144{154, 1993. [46] D.E Goldberg. Genetic Algorithms in Search, Optimization and Machine Learning. Reading MA Addison Wesley, 1989. [47] D.W Hearn. A dual ascent algorithm for trac assignment problems. Transportation Research B, 24B(6):423{ 430, 1990. [48] D.W Hearn, S Lawphongpanich, and S Nguyen. Pilot study of computer based urban trac management. Transportation Research B, 18B(4):357{365, 1984. [49] M.P Helme. Reducing air trac delay in a space-time network. IEEE, 1992. [50] C Hendrickson and G Kogur. Schedule delay and departure time decisions in a deterministic model. Transportation Science, 15(1):62{77, 1981. [51] J.K Ho. Solving the dynamic trac assignment problem on a hypercube multicomputer. Transportation Research, 24B(6):443{451, 1990. [52] J.H Holland. Adaptation in Natural and Arti cial Systems. University of Michigan press, 1975. [53] B.N Janson and J Robles. Dynamic trac assignment with arrival time cost. In Transportation and Trac Theory, 1993. [54] R Jayakrishnan, K.W.K Tsai, and A Chen. A dynamic trac assignment model with trac- ow relationship. Transportation Research, 3C(1):51{72, 1995. [55] R.C Jou and H.S Mahmassani. Comparability and transferability of commuter behavior chracteristics between cities : Departure time and route switching decisions. In Transportation Research Record (73 Annual Meeting, 1993. [56] D.E Kaufman, J Nonis, and R.L Smith. A mixed integer linear programming formulation of the dynamic trac assignment problem. IEEE, 1992. [57] S Kirkpatrick, C.D Gelatt, and M.P Vecchi. Optimization by simulated annealing. Science, 220(4598):671{ 680, 1983. [58] J.R Koza. Genetic Programming. MIT press, 1992. [59] S Lawphongpanich and D.W Hearn. Simplicial decomposition of the assignment problem. Transportation Research B, 18B(2):123{133, 1984. [60] L.J Leblanc, B Ran, and D.E Boyce. Dynamic travel choice models for urban transportation networks. IEEE, 1992. [61] M Lundy and A Mees. Convergence of an annealing algorithm. Mathematical Programming, 34:111{124, 1986. [62] H Mahmassani and R Herman. Dynamic user equilibrium departure time and route choice on idealized trac arterials. Transportation Science, 18(4):362{384, 1984. [63] H. Mahmassani and S Peeta. System optimal dynamic assignment for electronic route guidance in a congested trac network. In Proceeding of the Second International Capri Seminar on Urban Trac Networks, 1992. [64] H.S Mahmassani and G.L Chang. Experiments with departure time choice dynamics of urban commuters. Transportation Research, 20B(4):297{320, 1986. [65] H.S Mahmassani and R Jayakrishnan. System performance and user response under real-time information in a congested trac corridor. Transportation Research, 25A(5):293{307, 1991.

[66] F.L Mannering, S.A Abu-eisheh, and A.T Arnadottir. Dynamic trac equilibrium with discrete/continuous econometric models. Transportation Science, 24(2):105{116, 1990. [67] P Marcotte. Network optimization with continuous control science. Transportation Science, 17(21):181{197, 1983. [68] L Maugis. Mathematical programming for the air trac ow management problem with en-route capacities. IFOR, 1996. [69] D.K Merchant and G.L Nemhauser. A model and an algorithm for the dynamic trac assignment problems. Transportation Science, 12(3):183{199, 1978. [70] D.K Merchant and G.L Nemhauser. Optimality conditions for a dynamic trac assignment model. Transportation Science, 12(3):200{207, 1978. [71] N Metropolis, M Rosenbluth, M Rosenbluth, A Teller, and E Teller. Equation of state calculations for fast computing machines. Journal of Chemical Physics, 6:1087{1092, 1953. [72] Z Michalewicz. Genetic algorithms + Data Structures = Evolution Programs. Springer-verlag, 1992. [73] D Mitra, F Romeo, and A.L Sangiovanni-Vincentelli. Convergence and nite-time behavior od simulated annealing. Advances in Applied Probability, 18:747{771, 1986. [74] A Nagurney. An equilibrium scheme for the trac assignment problem with elastic demands. Transportation Research B, 22B(1):73{79, 1988. [75] A Nagurney. NETWORK ECONOMICS A Variational Inequality Approach, volume 1 of Advances in Computational Economics. Kluwer Academic, 1993. [76] A.R Odoni. The ow management problem in air trac control. In A.R Odoni et al, editor, Flow Control of Cogested Networks, volume F38 of ASI Series, pages 269{288. NATO, 1987. [77] M Papageorgiou. Dynamic modeling, assignment, and route guidance in trac networks. Transportation Research, 24B(6):471{495, 1990. [78] M. Papageorgiou. Concise encyclopedia of trac and transportation systems. Pergamon Press, 1991. [79] C.H Papdimitriou and K Steiglitz. Combinatorial Optimization : Algorithms and Complexity. Prentice-Hall, New York, 1982. [80] W.B Powell and Y She. The convergence of equilibrium algorithms with predetermined step size. Transportation Science, 16(1):45{55, 1982. [81] B Ran, D.E Boyce, and L.J LeBlanc. A new class of instantaneous dynamic user-optimal trac assignment models. Operations Research, 41(1):192{202, 1993. [82] O Richetta and A.R Odoni. Solving optimally the static ground holding policy problem in air trac control. Transportation Science, 27(3):228{238, 1993. [83] O Richetta and A.R Odoni. Dynamic solution to the ground-holding problem in air trac control. Transportation Research, 28A(3):167{185, 1994. [84] H.P Schwefel. Evolution and Optimum Seeking. Wiley, New York, 1995. [85] M.J Smith. An algorithm for solving asymmetric equilibrium problem with a continuous cost- ow function. Transportation Research B, 17B(5):365{371, 1983. [86] M.J Smith. The existence and calculation of trac equilibria. Transportation Research B, 17B(4):291{303, 1983. [87] M Terrab and A.R Odoni. Strategic ow management for air trac control. Operations Research, 41(1):138{ 152, 1993.

[88] D Van Vliet. The frank-wolfe algorithm for equilibrium trac assignment viewed as a variational inequality. Transportation Research B, 21B(1):87{89, 1987. [89] D.V Vliet and P.D.C Dow. Capacity-restrainned road assignment. Trac Engineering and Control, pages 296{305, 1979. [90] P Vranas, D Bertsimas, and A.R Odoni. The multi-airport ground-holding problem in air trac control. Operation Research, 42(2):249{261, 1994. [91] P.B.M Vranas, D Bertsimas, and A.R Odoni. Dynamic ground-holding policies for a network of airports. Transportation Science, 28(4):275{291, 1994. [92] P.C Vythoulkas. A dynamic stochastic assignment model for the analysis of general networks. Transportation Research, 24B(6):453{469, 1990. [93] J.G Wardrop. Some theoritical aspects of road trac research. In Proceedings of the Institute of Civil Engineers. Part II, pages 325{378, 1952. [94] B.W Wie, T.L Friesz, and R.L Tobin. Dynamic user optimal trac assignment on congested multidestination networks. Transportation Research, 24B(6):431{442, 1990. [95] B.W Wie, R.L Tobin, and T.L Friesz. The augmented lagrangian method for solving dynamic network trac assignment models in discrete time. Transportation Science, 28(3):204{220, 1994. [96] D.J Zawack and G.L Thompson. A dynamic space-time network ow model for city trac congestion. Transportation Science, 21(3):153{162, 1987. [97] S Zenios. Network based models for air trac control. European Journal of Operational Research, 50:166{178, 1991.

Abstract This paper addresses the general time-route assignment problem : One considers an air transportation network in a 2 dimensional space with asymmetric non-separable link cost and a eet of aircraft with their associated route and slot of departure. For each ight a set of alternative routes and a set of possible slots of departure are de ned. One must nd \optimal" route and slot allocation for each aircraft in a way that signi caly reduces the peak of workload in the most congested sectors and in the most congested airports, during one day of trac. A state of the art of the existing methods shows that this problem is usually partially treated and the whole problem remains unsolved due to the complexity induced. Genetic algorithm are then presented and adapted to the problem. New recombination operators are described and really improve the quality of the given solutions. Results on a large instance of the problem are then given and show how well genetic algorithm can manage this kind of combinatorial optimization problem.

Topics: Genetic Algorithm, Time of Departure and Route Optimization. Domain Area: Air Trac Control Status: Operational mock-up. Eort: 1 man/year Impact: Reduce the airspace congestion by changing the route and the time of departure of aircraft.

1 Introduction When joining two airports, an aircraft must follow routes and beacons ; these beacons are necessary for pilots to know their position during navigation and because of the small number of beacons on the ground they often represent crossing points of dierent airways. Crossing points may generate con icts between aircraft when their trajectories converge on it at the same time and induce a risk of collision. At the dawn of civil aviation, pilots solved con icts themselves because they always ew in good weather conditions (good visibility) with low speed aircraft. On the other hand, modern jet aircraft do not enable pilots to solve con icts because of their high speed and their ability to y with bad visibility. Therefore, pilots must be Centre d'Etudes de la Navigation Aerienne yMassachusett Institute of Technology

helped by an air trac controller on the ground who has a global view of the current trac distribution in the airspace and can give orders to the pilots to avoid collisions. As there are many aircraft simultaneously present in the sky, a single controller is not able to manage all of them. Then, airspace is partitioned into dierent sectors, each of them being assigned to a controller. Before taking o, when a ight is decided, the pilot has to le a Flight Plan (FPL) which describe the main steps of the ight (slot of departure, ight path, speeds, altitudes etc...), this FPL is necessary for the air trac control system to anticipate the ight monitoring. Usually, when there are few aircraft in the airspace the ight plan is well respected. As any human being, a controller has working limits, and when the number of aircraft increases, some parts of the airspace reach this limit and become congested. Two kinds of congestion can be identi ed according to the part of airspace involved : terminal congestion (around airports) and en route congestion (between airports). In the past, the rst way to reduce these congestions was to modify the structure of the airspace in a way that increases the capacity (increasing the number of runways, increasing the number of sectors by reducing their size). This has a limit due to the cost involved by new runways and the way to manage trac in too small sectors (a controller needs a minimum amount of airspace to be able to solve con icts). The other way to reduce congestion is to modify the ight plans in a way to adapt the demand to the available capacity. To reach this last goal three kinds of action can be applied :

Short term planning : Ground Holding, speed regulation and dynamic re-routing. Nowadays, when a delay

is expected to the arrival airport, a ground delay is applied in a way that will minimize the arrival airborne delay. This ground delay is safer (fewer aircraft waiting in the sky) and cheaper (according to the fuel consumption) than an airborne one. When the aircraft has taken o, this action can be completed by dynamic re-routing and speed regulation to adjust the actual demand to the present capacity. Long term planning : Modi cation of the slot of departure and the routes (static re-routing) of the whole

eet of aircraft. When we have a look at the congestion in the airspace one can remark that the peaks of congestion happen approximately at the same times during the day (one peak in the morning and one peak in the evening). The current time scheduling is usually xed by the airline companies without taking into account the ATC constraints. Then congestion is expected to be reduced by moving (in a limited domain) the time of departure of aircraft (in the past and in the future) and by changing the current ight paths (without too much extradistance). Slot and route fees. Another way to adjust the demand to the available capacity is to apply the principle : \Scarce resources are expensive". Then the departure slots and the arrival slots would have a higher price when the demand is strong and in a similar way the fastest route would have a higher price too. This action might help to spread the demand over time and over space in an equitable way.

This paper is focused on long term planning and shows how stochastic optimization is able to manage this kind of problem. In a rst part, a description of the previous related works is given and in a second one, a simpli ed model is developed and a mathematical formulation of our problem is given. In the third part a description of stochastic optimization is given and particularly of a genetic algorithm. The application of the algorithm is given on the fourth part and nally an improvement is given in the fth part.

2 Previous Related Works

2.1 Introduction

Trac assignment techniques have been developed in a way to reduce congestion in transportation networks by spreading the trac demand in time and in space. Historically, those techniques have been applied to road trac assignment due to the strong level of congestion encountered in this domain. The problem of users of a congested transportation network seeking to determine travel paths of minimal cost from origins to their respective destinations is a classical network problem. An equilibrium occurs when the number of trips between an origin and a destination equals the associated travel demand. Wardrop [93] stated the trac equilibrium conditions through two principles :

First Principle : The journey times of all routes actually used are equal, and less than those which would

be experienced by a single vehicle on any unused route. Second Principle : The average journey time is minimal.

From those two principles, Dafermos and Sparrow [26] coined the terms user-optimized and system-optimized transportation networks to distinguish between two distinct situations in which users act unilaterally, in their own self-interest, in selecting their routes, and in which users select routes according to what is optimal from the societal point of view, in that the total costs in the system are minimized. More detailed descriptions of those equilibria can be found in the following references : [38, 33, 86]. According to the separability1 of the link cost function, dierent algorithms have been developed and the associated references can be summarized in the following array :

symmetric non-separable link cost

User-optimized

System-optimized

[80, 88, 67]

[47, 43]

asymmetric non-separable link cost [34, 35, 85, 24, 48, 59]

[47, 43]

Improved trac assignment algorithms, take into account that demand can be in uenced by the induced cost but still with a static demand over time [41, 42, 25, 74]. By de nition all those techniques, are applied to static trac demand and are mainly used to optimize trac on a long time period and can only capture the macroscopic events. When a more precise matching between trac demand and capacity has to be found, microscopic events have to be taken into account, and dynamic trac assignment techniques have to be used ([78] gives a good description of those techniques). Those techniques try to nd an optimal route, or an optimal time of departure, or both, for each individual in a way to reach a dynamic user-equilibrium or a dynamic system-equilibrium. The complexity induced by the dynamic trac assignment is strong, especially when route and time of departure are simultaneously optimized. This problem is NP HARD [30, 11] and may have several optima [44]. Even when only the route decision is investigated, this problem does not have an exact solution when asymmetric non-separable link costs are used. Due to this complexity, the problem is always partially solved for simpli ed instances :

Ben Akiva, M [9] : \1" origin and \1"destination : time decisions only ; Hendrickson, C et al. [50] : \n" origins \1" destination, system-equilibrium, analytical approach with

simple link cost (time decisions only)) ; Mahmassani, H et al. [62, 22, 55] : \1" origin and \1" destination analytical approach with simple link cost (time decisons only) : Merchant, D.K et al. [69, 70] : network with separable link cost, route choices only (heuristic) ; Alfa, A.S [3] : heuristic for time-route trac assignment, the \convergence" to a xed point is not ensured ; Arnott, R et al. [7] : N parallel routes, time-route decisions, separable link cost function ; HO, J.K [51] : route choices, separable link cost, decomposition method ; Bell,M.G.H [8] : route choice, link with capacity, separable link cost : Carey, M [17] : route choices, separable link cost, system-equilibrium reached by convex programming.

1

separable link cost : ) cij = cij (fij ) 8 link (i; j ) ~ symmetric non-separable link cost : ) cij = cij (f~) (f~ ow on the whole network) and @[email protected] (f ) = @[email protected] (f~) ~ asymmetric non-separable link cost : ) cij = cij (f~) and @[email protected] (f ) 6= @[email protected] (f~) ij

kl

ij

kl

kl

ij

kl

ij

In the same time, speci c approaches have been developed to solve this route-time allocation problem. The main ones are the following : Space-time network ; Variational Inequality ; Optimal Control ; Simulation ; Dunamic Programming (Ground Holding Problem).

2.2 Space-Time Network This technique consists of expanding the space network into a space-time network by creating N space-time nodes for each space node (where N is the number of duration periods) and then assigning the dynamic demands using classical static trac assignment methods. The associated number of links depends on the speed of the vehicle traveling on the network which is the weak point of this model [96, 31, 49, 56, 54, 97].

2.3 Variational Inequality The variational inequality problem is a general problem formulation that encompasses a plethora of mathematical problems, including, among others, nonlinear equations, optimization problems, complementarity problems, and xed point problem. It is now known that the simultaneous route and departure-time choice equilibrium problem can be formulated as a variational inequality problem [75]. Bernstein, D et al. [15, 39, 13] successfully used this method to solve the simultaneous route and departure-time choice equilibrium problem in networks with separable link costs. This approach is still theoretical and can be applied only to very small networks.

2.4 Optimal Control The optimal control approach has been tried for solving small instances of the general route-time assignment problem. Usually, the state variables are the link ows at time t and the associated control variables are the ows entering links during the same period of time.

Wie, B.W et al. [95] : "n" origins and \1" destination, route choices only, system-equilibrium, separable

link costs ; Wie, B.W et al. [94] : multi origin-destination, route choices, user-equilibrium, separable link costs ; Papageorgiou, M [77] : route choices, user-equilibrium and system-equilibrium, separable link costs ; Ran, B et al. [81] : multi origin-destination, route choices, user-equilibrium, separable link costs ; Friesz, T.L et al. [40] : \n" origin and \1" destination, route choices, user-equilibrium, system equilibrium ; Leblanc, L.J et al. [60] : optimal control modeling route-time choices user-equilibrium (modeling only, no algorithm) ; Janson, B.N et al. [53] : route-time choices, separable link costs ;

Optimal control techniques can be used only for small networks.

2.5 Simulation Due to the complexity induced by the dynamic route-time assignment in a general network (asymmetric nonseparable link costs), simulation has been investigated in a way to reach a solution close to the user-equilibrium or system-equilibrium. The strong point of the simulation techniques is their ability to take into account the maximum features of the real problem and can be applied on general networks. Cascetta, E et al [20, 19, 21, 18], Cantarella, G.E et al. [16] and Vythoulkas, P.C [92] applied simulation techniques in the following way : A set of possible slots of departure and a set of possible routes is assigned to each user. Each user randomly selects a couple of decision variables and the trac is simulated for the whole period considered (usually a day). A user speci c utility function is associated with each selected individual choice after the simulation. The choices of the day k are then based on the utilities of the N previous days using a kind of learning process (limited memory from the day k ? N to the day k). It can be shown that this stochastic process is a stationary Markov chain that reaches a xed point. This xed point is very stable but its distance from the user-equilibrium point or the system-equilibrium point can not be quanti ed and depends of the initial conditions. Others works based on the same methods can be found in the following references : [10, 64, 12, 63, 65, 66, 89, 23]

2.6 Dynamic Programming (Ground Holding Problem) The ground holding Problem was rst addressed and formalized by Odoni, A.R in [76]. Andreata, G and Romanin-Jacur, G [5] developed a dynamic programming framework to compute optimal delays for n ights arriving in the same time period at a single destination with a deterministic capacity scenario (a description of airport capacity scenarios is given in [45]). Terrab, M and Odoni, A.R [87] solved the multiple time periods, single arrival airport case as a minimum cost ow problem in a capacitated network, and as an assignment problem. They also developed a dynamic programming model to cope with stochastic capacities. However, this model was impractical for typical ow instances. The problem \single airport and random capacities" was solved by Richetta, G and Odoni, A.R [83, 82]. Vranas, P et al. [91, 90] addressed the multi-airport case with a deterministic airport capacity scenario (if there are not connected ights the multi-airport problem can be decomposed into several single airport problem). A good summary of those previous approaches is given in [4]. Bertsimas, D.J and Stock, S [14] proposed 0-1 programming models which incorporate airport capacities, sector capacities and connected ights. The model only assigns delays to aircraft and does not address route decisions. Finally Maugis, L [68] extended the previous model to take sector constraints in a more smoothly way.

2.7 Conclusion From the trac assignment point of view, the problem we have to solve could be summarized in the following way :

a eet of aircraft is considered for which one must nd the optimal route and the optimal time of departure in a way to reach a \system-equilibrium" in a network with asymmetric non-separable link costs

All the previous approaches are not able to manage the whole problem due to the complexity induced. In the following, a model is proposed and a method is developed in a way to give \very good" solutions for realistic instances of the whole problem.

3 A Simpli ed Model 3.1 Introduction

Congestion in the airspace is due to aircraft which have close positions in a four dimensional space (one time dimension and 3 space dimensions). It is then relevant to investigate ways to separate those aircraft in this 4 dimensional space by changing their slot of departure (time separation) or by changing their route (spatial separation) or both. Those changes must be done in a way that takes into account the objectives of the airlines : the moving of the slot of departure must be done in a limited domain (otherwise, for instance, some aircraft will be forced to take o at 2:00AM to reduce the congestion of the day but at this time there will be no passengers to carry);

Workload Morning peak

Evening peak

Time (hours) 4

5

6

7

8

9

10

11 12

13 14

15 16

17 18

19

20 21

22 23

0

1

Figure 1: Congestion variations

the new slot of departure must take into account the connexions between ight (as a mater of fact some

aircraft have to wait the arrival of some previous ights to take o (hub phenomenon); the possible routes must not generate too large additional distances. So, for each ight, a new pair (slot of departure, route) will be chosen from two discrete and nite sets : a set of possible slots of departure (around the original slot of departure) ; a set of routes which do not increase too much the total path length and are approved by the airline company the ight belongs to. After taking o the aircraft will follow its ight path and will generate workload in the dierent sectors encountered and also a congestion at the arrival airport when arriving. According to the controllers themselves, the workload induced in a control sector is a function of the three main following criteria : the con ict workload that results from the dierent actions of the controller to solve con icts. the coordination workload corresponds to the information exchanges between a controller and the controller in charge of the bordering sector or between a controller and the pilots when an aircraft crosses a sector boundary; the monitoring aims at checking the dierent trajectories of the aircraft in a sector and induces a workload. From the airport point of view the congestion is de ned by the number of aircraft unable to land due to the runways use (so it depends on the number of aircraft waiting in the sky). This airport congestion is dependent on the available capacity which changes over the day according to the weather. A static capacity has been attributed to each airport for our experiments, but the model is able to manage dynamic airport capacity, as well. If one sector is considered, one can observe variations of the associated workload during the day according to the trac in it. So, a typical pro le of the sector control workload is very chaotic with dierent peaks correlated with the peaks of trac in it (usually one peak in the morning and one peak in the evening see gure 1). It must be noticed that the peaks do not appear at the same time in the dierent sectors due to the ow propagations. According to our de nition of workload, only the information related with the route and the slot of departure are relevant from the ight plan. As a mater of fact, this information enables us to compute the time when an aircraft enters (or exits) a sector and when it is over a crossing of two airways (point of con ict generation). When examining the physical air trac network, we notice that airways are superpositions of several ight routes which have the same projection on the ground but dierent altitudes according to their azimuth (semi circular rule2 ). So an airway can be modeled by a bidirectional link which combines several individual aircraft routes. Then, without loss of generality, the real 3-dimensional transportation network will be modeled by a classical 2-dimensional network on a horizontal plan. We can now de ne our goals more precisely in the following way : 2 Aircraft with heading between 0 and 180 have to y at odd altitudes (in thousands of feet) and even altitude for headings between 180 and 360

A

13

D

18

7 2

17

12

19

6

1

16

20 21

5 11 3

B

8 15

C

22

24

10

4 23

14

9

E Beacon events

1

5

11

15

22

24 Sector events

B

C

E

Figure 2: Simpli ed Flight Plan one considers an air trac transportation network in a 2 dimensional space and a eet of aircraft with their associated route and slot of departure. For each ight a set of alternative routes and a set of possible slots of departure are de ned. One must nd \optimal" route and slot allocation for each aircraft in a way that signi cally reduces the peak of workload in the most congested sectors and in the most congested airports, during one day of trac. This \optimal" bi-allocation (route-time) must take into account the objectives of the airlines.

3.2 Mathematical formulation A pair of decision variable (i ; ri ) is associated with each ight in which i is the advance or the delay from the original slot of departure and ri is the new route. With this notation (0; r0 ) will be considered as the most preferred choice from the user point of view. Those two decision variables (i ,ri ) will be chosen from two nite-discrete sets : for the slots and R for the routes. More precisely the structures of and R are the following : = ?m; ?m + 1; ::::; ?1; 0; 1; :::; p ? 1; p R = r0 ; r1 ; r2 ; :::; rmax For which m ; p are respectively the maximum advance and the maximum delay permitted for a ight. Those limits can be dierent for each ight. The routes are ordered according to cost induced to the associated ight. So r0 is the best one and rmax the worst. To be able to compute the control workload (according to our de nition), the route of an aircraft will be de ned by the two following lists :

list of the over own beacons in the air network with the associated passing time : Lb = (b1 ; t1 ); (b2 ; t2 ); :::; (bi ; ti ); :::

list of crossed sectors with the associated in-out times : Ls = (S1 ; Tin1; Tout1); (S2 ; Tin2; Tout2); :::; (Sk ; Tink ; Toutk ); ::: An example of a simpli ed ight plan is given in gure 2. In this case an aircraft is joining airports 1 and 24. It over ies 5 beacons (5,11,15,22) and crosses 3 sectors. Given this information, it is now possible to compute the three components of the control workload (monitoring, coordinations, con icts) in each sector of the airspace. To reach this goal, one rst begins to compute the exiting

ow on each link of the network and the crossing ow over each sector. With the de nition of our simpli ed ight plan, it is easy to compute the discrete event when an aircraft exits a link of the network or when an aircraft enters or exits a sector. The ows are then computed by using a moving window in which an average of the number of events is computed (the temporal size of the window is a parameter of the model). So the ow fij (t) exiting a link (i; j ) at time t is computed in the following way :

fij (t) = 2:D1+ 1

s=X t+D nX =N s=t?D n=1

t ijn

where D is the half extension of the temporal window (in slot) and ijn (t) is the kroneker symbol which is equal to 1 if the ight n exits the link (i; j ) during the slot t (N is the whole number of ight). Given the ows on links, it is now possible to give a mathematical formulation of the control workload induced in a sector. As it has been previously said, workload in a sector Sk at time t can be roughly expressed by the summation of three terms : WSt = WmoS (t) + WcoS (t) + WcfS (t); k

k

k

k

WmoS (t) Monitoring workload ; WcoS (t) Coordination workload ; WcfS (t) Con ict workload. k

k

k

3.2.1 Monitoring workload

The monitoring workload is directly related to the number of aircraft in a sector and can be expressed in the following way : WmoS (t) = :n S (t) weight coecient and n(t) is the average number of aircraft in the sector Sk at time t. k

k

n S (t) = 2:D1+ 1 : k

l=X t+D l=t?D

nS (l) k

nS (t) is the number of aircraft in the sector Sk during slot t. nS (t) = nS (t ? 1) + NinS (t) ? NoutS (t) NinS (t) is the number of aircraft entering sector Sk during the slot t and NoutS (t) is the number of aircraft exiting sector Sk during the slot t. k

k

k

k

k

k

k

3.2.2 Coordination workload

During a slot t the number of coordinations induced in a sector Sk is equal to the number of aircraft entering or exiting the sector during the same time period. So, it is related with the ows cut by the sector boundaries : WcoS (t) = (FinS (t) + FoutS (t)) where is a weight coecient and FinS (t) and FoutS (t) are the entering and exiting sector ows associated with Sk at time t. k

k

k

FinS (t) = 2:D1+ 1 k

l=X t+D l=t?D

k

k

NinS (l); FoutS (t) = 2:D1+ 1 k

k

l=X t+D l=t?D

NoutS (l); k

NinS (t) and NoutS (l) are the number of discrete events \an aircraft enters sector Sk ", \an aircraft exits sector Sk " during the slot l. k

k

δAm

δAp

A

δm

δp

B

B

δm

Previous Flights (Arriving Intervals)

B

δp

C

C

δm

C

δp

D

D

D

τE δEp

δEm

E

Initial Interval

δEm

δEp

Connected Flight (Departure Interval)

Updated Interval

τE

Figure 3: Example of connexion

3.2.3 Con ict workload

The con ict workload induced in a sector Sk depends on the crossing ows at each node in the sector. X X X [h(ijl )fij (t)flj (t)] WcfS (t) = 12 j 2NO i 2 NO l 2 NO i 6= j l 6= i; l 6= j where : is a weight coecient ; NOS is the set of nodes belonging to the sector Sk ; NO is the entire set of nodes in the network ; ijl is the crossing angle of the links (i; j ) and (l; j ) ; k

Sk

k

h(ijl ) is a weighted function to take into account the fact that a con ict is harder to solve when the angle of crossing is small.

3.2.4 Formulation of the objective function

With the previous de nition of the control workload , when a set of decision variables has been chosen, it is now straightforward to compute the dynamic workload pro le over time for each sector and to memorize the associated maxima. We de ne our objective in the following way : \ one must try to reduce congestion in the most overloaded sectors but not necessarily simultaneously in all the sectors" ; this will spread the congestion over several sectors. So, we have :

obj = min

kX =K k=1

max W (t) t2T S

k

3.3 Connected ights constraint When a ight has to be connected to some arriving ights, its slot of departure must be later than the time of arrival of the previous ight and separated from it by a minimum amount of time ( ). So, when a slot of departure is changed, one must rst check that this new scheduling respects the connecting constraint. To reach this goal, the slot moving set of the connected ight must take into account the slot moving sets of the previous

ights as it can be seen on gure 3. On this gure, the ight \E" is connected to 4 previous ights :"A,B,C,D". To always be able to assign a slot to a ight, the slot moving sets must satisfy the following properties :

Maximum boundary

? pE > k2fA;B;C;D max g Pk + E

This checking must be done by using the longest possible route for each previous ights A; B; C; D. Minimum Boundary ? mE > k2fA;B;C;D max g tka + E where tka is the actual arrival time, determined by the slot moving and routes assigned to the ight A; B; C; D. So, before assigning a moving slot to a connected ight, one must rst check that the previous associated ights have already been assigned a slot and a route from their own route and moving slots sets. This will ensure that the random point generated in the state domain will always satisfy the connexion constraint.

3.4 Problem complexity Before investigating an optimization method, the associated complexity of our problem must be studied. The model previously developed is discrete and induces a high combinatoric in a huge space. As a mater of fact, if Rn ; n are the route set and the slot moving set associated with ight n, the number of possible choice (Cn ) for this ight is : Cn = jRn j:jn j where jS j denotes the cardinality of the set S . Then the number of possible decision variable combinaisons is given by :

jStatej =

nY =N n=1

Cn

where N is the total number of ights. If the number of choices is the same for each ight (jRn j = C te = r and n j = C te = ), then the cardinality of the state domain becomes : jStatej = (r:)N For instance, for 20000 ights with 10 route choices and 10 possible slot movings, : 10020000. Moreover, those decision variables are not independent due to the connexion induced by the control workload, the airport congestions and the connexion constraint, so, decomposition methods cannot be applied. It must be noticed that the objective function : is not continuous (then it is not convex) : may have several equivalent optima (the objective is then said multimodal) ; This problem is then a strong NP HARD problem with non-separable state variables.

3.5 Conclusion The microscopic model previously developed takes into account the major features of the real problem with the associated constraints and is able toaccommodate speci c individual preferences. This model has been completed by a mascroscopic description of the sector workload which induces a robust computation of sector congestion. So, one must solve a multimodal combinatorial optimization problem in a huge space with non-separable decision variables. Currently, only stochastic optimization is well adapted to address this kind of problem.

4 Stochastic Optimization

4.1 Introduction

Stochastic optimization techniques usually addressed problem with strong complexity for which the objective function has no speci c feature such as convexity, continuity etc.... The main characteristic of those stochastic optimization techniques is to randomly move in the state domain in a way to improve the objective criterium. At the beginning, local search algorithms randomly moved in the state domain toward the closest optimum without being able to avoid local optima. As a mater of fact, this technique consist in rst selecting a random point in the state domain, second to nd a neighbor point and to compare the associated objectives : if the neighbor is better it begins the currant point else one must nd another neightbor point till the objective is improved. Afterward, more sophisticated and more powerful methods have been developed in a way that avoid local optima and converge toward global optima (Simulated Annealing and Genetic Algorithm) and which can identi ed several global optima (Genetic Algorithm). Those algorithm are mainly used in combinatorial optimization we now describe.

4.2 Combinatorial Optimization Combinatorial optimization problems represent a class of computational problems where one seeks the one \best" or \optimal" solution from among a large space of solutions [79]. Generally, the optimal solution is speci ed in terms of minimizing some given cost function subject to a set of constraints. Formally, a combinatorial optimization problem can be stated as a tuple < S; f > , where S , the state space, is the ( nite) set of possible solutions and f is a cost function of the form : f :S ??>R mapping each possible solution onto the reals. The optimal solution to any such combinatorial optimization problem is a solution, sopt 2 S , such that : f (sopt ) f (s); 8s 2 S sopt is called a globally optimal solution. In some cases, it can be more natural to formulate the optimal solution as that solution which maximizes, rather than minimizes, some quantity. Such cases can be readily handled, without loss of generality, within the formalism here by simply noting that maximization is equivalent to minimization with a reversal of sign in the cost function. Consider some solution, s 2 S , and the set of all solutions, Ss , which are in some sense \close" to s. The set Ss is called the neighborhood of s . We can de ne the neighborhood structure of S as a mapping N : S ! 2S which formalizes the notion of \closeness" for a particular optimization problem. If < S; f > is an optimization problem and N is a neighborhood structure, then s^ 2 S is called a local minimum, with respect to N provided that : f (^s) f (s); 8s 2 N The relationship between globally and locally optimal solutions can be visualized in the 2-dimensional case as in gure 4. One class of approximation algorithms used in solving optimization problems is known as local search. Local search algorithms operate by stepwise re nement based on a neighborhood structure. Given an initial solution, usually chosen at random, a local search algorithm evaluates the cost function on the initial solution. Next, subsequent candidate solutions are generated in the neighborhood of the currently optimal solution. If the cost of the newly generated candidate solution is less than the currently optimal solution, it becomes the currently optimal solution. This process is iterated until a satisfactory locally optimal solution is found, no improvement has been made for some number of previous iterations, or until some maximum number of iterations has been reached. By de nition, then, local search algorithms terminate in locally optimal solutions. Unless the neighborhood structure is de ned such that all locally optimal solutions are also globally optimal, then it is clear that local search cannot be guaranteed to result in a globally optimal solution. Furthermore, there is no easy answer to the question of how close the locally optimal solution found by any particular local search algorithm is to a globally optimal solution. In general, the quality of the solution found by local search is highly dependent on both the initial solution selected and the method of generating candidate solutions from the neighborhood of any particular solution. Those drawbacks can be avoided by using more sophisticated stochastic optimization techniques such as : genetic algorithms.

f(x)

Local optima

Global optimum x

Figure 4: Global vs. Local Minima

4.3 Classical Genetic Algorithm Genetic algorithms (GAs) [46, 52, 37, 72, 58, 36, 84] are problem solving systems based on principles of evolution and heredity. A GA maintains a population of individuals,P (t) = x1 ; x2 ; :::; xn at iteration t. Each individual represents a potential solution to the problem at hand and is implemented as some (possibly complex) data structure S . Each solution xi is evaluated to give some measure of tness. Then a new population at iteration t + 1 is formed by selecting the tter individuals. Some members of the new population undergo transformation by means of genetic operators to form new solutions. There are unary transformations mi (mutation type), which create new individuals by a small change of a single individual and higher order transformations ci (crossover type), which create new individuals by combining parts from several (two or more) individuals. For example, if parents are represented by a ve-dimensional vector (a1 ; a2 ; a3 ; a4 ; a5 ) and (b1 ; b2 ; b3; b4 ; b5 ), then a slicing crossover of chromosomes after the second gene produces ospring (a1 ; a2 ; b3 ; b4; b5 ) and (b1 ; b2 ; a3 ; a4 ; a5 ). The control parameters for genetic operators (probability of crossover and mutation) need to be carefully selected to provide better performance. The intuition behind the crossover operation is information exchange between dierent potential solutions. After some number of generations the program converges - the best individual hopefully represents the optimum solution. In [46], Goldberg presents a good summary of many recent GA applications in biology, computer science, engineering, operations research, physical sciences, and social sciences. Genetic algorithms use a vocabulary borrowed from natural genetics. We would talk about genes (or bits), individuals (chromosomes or gene strings), and population (of individuals). Populations evolve through generations. The execution steps of genetic algorithms are the following(see gure 5) : 1. Initialize population and evaluate tness : To initialize a population, we need rst to decide the number of genes for each individual and the total number of chromosomes (popsize) in the initial population. All the information contained in a state point must be summarized in the chromosome, this results from the coding which is one of the success keys of genetic algorithms. 2. Reproduction (Selection) : Reproduction is the selection of individuals with respect to the probability distribution based on the the tness values. Fitter individuals have better chances of getting selected for reproduction [72]. Dierent kinds of selection processes can be applied for a speci c problem; one of the most simple is the \Roulette Wheel Selection Process" : Each chromosome has a certain number of slots proportional to its tness value. The selection process is based on spinning the wheel many times; each time we select a single chromosome for a new population. Obviously, some chromosomes will be selected more than once. This is in accordance with genetic inheritance: the best chromosomes get more copies, the average stay even, and the worst die o. 3. Crossover : Now we are ready to apply the rst recombination operator, crossover, to selected individuals. The probability of crossover,pc , gives us the expected number pc:popsize of chromosomes which should undergo

population(k) selection P1

P2

crossover C1

C2

mutation

mutation

C1

C2

fitness evaluation

fitness evaluation

population(k+1)

Figure 5: Classical Genetic Algorithm the crossover operation. After having selected two parents (according to their tness), we generate a random number r between 0 and 1; if r < pc , then the parents undergo a crossover. We then generate a random number pos from the range of (1::m ? 1), where m is the total number of genes in a chromosome. The number pos indicates the position of the crossing point. The coupled chromosomes exchange genes at the crossing point as described earlier. 4. Mutation : Chromosomes which have not undergone a crossover may be randomly mutated on the gene by gene basis. To perform a mutation, a gene is randomly selected and modi ed by a problem dependent process in a way to produce new genes. The probability of mutation,pm , gives us the expected number of mutated gens pm :m:popsize. 5. Convergence : Following reproduction, crossover, and mutation, the new population is ready for its next generation. The rest of the evolutions are simply cyclic repetitions of the above steps until the system reaches a predetermined number of generations or converges (i.e., no improvement in the overall tness of the population). This classical scheme can be improved by using simulated annealing in the crossover operators.

4.4 Genetic Algorithm with Simulated Annealing Crossover This Simulated Annealing Crossover has been developed recently [28] and has been used in dierent applications with success [27, 29, 32]. It greatly improves the convergence speed of the genetic algorithm for all applications, so it is generic and not limited to a problem dependent improvement. This operator is based on another stochastic optimization technique called : Simulated Annealing.

4.4.1 Simulated Annealing

Simulated annealing is based on an analogy to the thermal process of metallurgical annealing from statistical mechanics; an examination of this process will help to better understand the simulated annealing process [57]. In metallurgical annealing, the goal is to bring a solid, in contact with heat bath, to a con guration with minimal energy. Distortions in the crystal lattice structure of the metal result in higher energy states. The central idea behind metallurigal annealing is that disruptions in the crystal lattice structure of the metal are

removed by thermal agitation at very high temperatures and that further disruptions can be prevented from forming if the metal is cooled slowly enough. Metropolis, et al. [71] succeeded in simulating the evolution of states a metal undergoes as it reaches thermal equilibrium, at constant temperature, using Monte Carlo methods. The technique, which has come to be known as the Metropolis Algorithm (MA), is as follows. Given a current state, generate a new state by randomly perturbing the current state (thus, simulating the random motion of particles). If the new state results in lower overall energy then retain the change as the current state; if the result is a higher energy state, then retain the new state with probability P , given by :

p = e? where E = Enew ? Eold is the energy change, T is the temperature of the heat bath and kb is the Boltzmann constant. For a system in thermal equilibrium at temperature T , the probability of nding the system in a state, s, with energy Es is given by E kb :T

PT (s) = Z (1T ) e?

Es kb :T

where Z (T ) , the so-called partition function, is given by

Z (T ) =

X

i

e?

Ei kb :T

where the summation is over all possible states of the system. PT (s) is known as the Boltzmann distribution. The similarities between the local search technique presented in the combinatorial optimization section and MA suggest the possibility of using MA to escape local minima. If < S; f > is an instance of an optimization problem and s1 ; s2 2 S are solutions with associated costs f (s1 ) and f (s2 ), then the acceptance criterion by which s2 is retained is determined by the acceptance probability : 1 if f (s2 ) f (s1 ) PT (s1 ! s2 ) = e ( 1 )? ( 2 ) otherwise f s

f s

T

where T 2 [0; 1] is a monotonically decreasing control parameter. In addition to the acceptance criterion and probability, full speci cation of the simulated annealing algorithm also requires a sequence fT0 = 1; T1 ; :::; Tk ; :::; Tfinal > 0g, called the annealing schedule and ,uk a further control parameter stipulating the number of solution candidates, or transitions, generated at T = Tk . So, given an optimization problem, < S; f >, then under reasonable restrictions on the neighborhood structure and a suciently large number of transitions at a xed value of the control parameter T , it has been shown [2] that the simulated annealing algorithm, using the acceptance probability, will nd a solution, s 2 S , with probability : () PT (s) = N 1(T ) e? 0 where N0 (T ) is given by : f s T

N)(T ) =

X

i2S

e?

()

f i T

PT (s) is known as the stationary distribution of the solution s , the equivalent of the Boltzmann distribution for physical systems. Denote the set of all globally optimal solutions to < S; f > by Sopt . It has been shown [61], for s 2 Sopt , that : lim P (s) = 1 (s) T !0 T

Sopt (S ) where (S ) (s) is de ned in the following way : Let A and A0 2 A be two sets. Then (A0 ) (a) = 1 i a 2 A0 , opt

opt

and 0 otherwise. The above result is extremely valuable, because it guarantees that the simulated annealing algorithm will asymptotically converge to Sopt , provided that the stationary distribution is reached for each value of the control parameter, T . However, asymptotic convergence is only possible given an in nite number of transitions [6]. Clearly, an in nite time convergence is not of practical interest; although it leads to the question of the eciency of approximating the asymptotic behavior. At least two important results in this regard have been

P1

P2

CLASSICAL CROSSOVER

P1

P2

C1

Pc

C2

TOURNAMENT

C1

C2

Figure 6: Simulated Annealing Crossover reported. Aarts and van Laarrhoven [1] have shown that the stationary distribution is approximated arbitrarily closely only under the condition that the number of transitions is at least quadratic in jS j , the cardinality of the solution space. Attacking the problem dierently, Mitra, et al. [73] have shown that for many problems, arbitrarily close approximation to the asymptotic behavior requires a number of transitions which is larger than jS j , leading to exponential-time execution. That is, enumerating all possible solutions requires less time than needed to perform the simulated annealing procedure until an arbitrarily close approximate solution is reached. We now return to consideration of the annealing schedule and uk , the number of transitions attempted at each value of T . By suitably choosing these two parameters (T and Uk ), polynomial-time execution is obtainable [1], at the expense of not being able to guarantee the degree of divergence (the distance) between the asymptotically optimal and the approximated solution. In general, the annealing schedule is chosen according to the following guidelines : 1. the initial value, T0 , is such that PT0 (s) = 1; 2. the decrement rule, Tk+1 = Tk ; k = 1; 2; :: where < 1, and typically 0:9 < < 0:99 ; 3. the nal value, or stopping criterion, Tfinal , is that value at which suciently many successive transitions have not resulted in a signi cant improvement in reducing the cost function The choice of initial value is motivated by noting that such a value results in virtually all attempted transitions being accepted, thus approximating P1 . It must be noticed that the larger the dierence between Tk and Tk+1 the more transitions, uk , that will be required to approximate the stationary distribution at Tk+1 [57, 2].

4.5 Simulated Annealing Crossover The way to introduce Simulated Annealing in the crossover operator is straightforward. Two parents (P1 ; P2 ) are selected in the population according to their tness and randomly crossed using the classical crossover operator. If the crossover happened, two children are generated (C1 ; C2 ). The rst child (C1 ) is then randomly associated with one parent and the second one (C2 ) with the remaining parent. For each created couple, the \winner" is then selected according to a simulated annealing process in the following way (see gure6) : if f (C ) > f (P ) ) C 0 = C ( )? ( ) if f (C ) < f (P ) ) C 0 = C with the probability P (C ! C 0 ) = e ( ) where T (gen) is the control parameter depending of the current generation gen. f C f P T gen

4.6 Conclusion The strong point of these stochastic optimization methods is their ability to give \very good" solution(s) to large instances of real problems. Usually the given solutions are very robust and very close to the optimum when it is known. They can manage very speci c constraints and do not need any properties concerning the objective function.

SECTOR CONGESTIONS

Traffic Simulator

Genetic Algorithm

FPL

CHROMOSOME

Flight Plans Generator

Figure 7: GA using simulated tness R1

n rk

1 r3

∆n

N r6

δnj

δ12 ∆1

RN

Rn

δN3 ∆N

∆n

: Set of slots for the flight n

Rn : Set of routes for the flight n

Figure 8: Structure of the Chromosome

5 Application to our problem 5.1 Introduction

As previously noted, Genetic Algorithms need a data coding, speci c recombination operators and a tness evaluation to be able to select best individuals. The way this speci c genetic algorithm works is the following. A set of ight plans is generated from each chromosome candidate and the whole associated day of trac is generated. Sector congestion and airport congestion are registered and the associated tness is computed (see gure 7). This is an appropriate method to use because the more realistic the simulation of the day of trac is the more relevant are the results. The problem speci c features of the Genetic Algorithm are now described.

5.2 Data Coding This step consists of converting each point of the state domain into a chromosome used by the genetic algorithm. In our problem, the state variables (which contain all the information needed to compute the sector workload) consist of the set of ight plans. The possible new path and new slot moving have been supposed to be chosen in two discrete- nite sets associated with each ight. In this case a straightforward coding has been used in the sense that each chromosome is built as a matrix which gather together the new slot moving (for the time of departure) and the new route number (for the ight path) (see gure 8). With this coding, a population of individuals can be created by randomly choosing a new slot moving number and a new route number from individual sets associated with each ight.

SLICING CROSSOVER

Figure 9: Crossover operator

5.3 Fitness Evaluation The previous description of the chromosome and the associated recombination operators can now be used to randomly generate some points in the state domain. To apply the selection operator, a tness must be associated with each chromosome in a way to evaluate the quality of each individual according to the optimization criterion. In our problem, one rst considers the default ight plans and computes the maximum workload in each sector in the airspace.

MaxW (ref ) =

kX =K k=1

max W (t; ref ) t2T S

k

Afterward, for each generated chromosome, maximum sector workloads are computed the same way :

MaxW (chrom) =

kX =K k=1

max W (t; chrom) t2T S

k

The tness is then de ned by the ratio of the two previous expressions. MaxW (ref ) fitness(chrom) = MaxW (chrom) So, when fitness(chrom) > 1, it means that the induced congestion is lower than the reference one.

5.4 Crossover Operator The crossover operator has been implemented as a slicing one on the rst attempt. So, to make a crossover, one rst select two parents and randomly chooses an intergene site from which the two induced ospring are exchanged (see gure 9).

5.5 Mutation Operator In our problem, the mutation is used to generate new slot moving or new route from the two discrete- nite sets associated with each ight. Then, a new couple (slot moving,route number) is generated from each individual associated set of slot moving and route for all the ights (see gure 9).

n

Rn

∆n

MUTATION

Figure 10: Mutation operator

6 Application to a simpli ed network 6.1 Introduction

A simpli ed network has been used to evaluate the performances of the algorithm. This example tries to reproduce the trac density encountered usually in airspace. The evaluations have been done with ight plans generated in a two dimensional network but the model could be applied to the real ight plans in three dimensional airspace, due to the fact that only ight plan events are used to compute the objective function. In this example, only sector congestion has been taken into account but as it has been mwntioned before, airport congestion could have been introduced by using airports with smaller capacitiey. The features of the network are the followings : Number of nodes : 100 Number of links : 544 Number of Origin-Destination pairs : 30 Number of Sectors : 15 Number of Airports : 24 Geographical size of the network : 1000NM x 1000NM A physical description is given in gure 11 with the sectoring and the associated sector number. Each airport is symbolized with a larger node. A hundred ights have been randomly generated for each Origin-Destination pair, so, a total of 3000 ight plans is used for this example (the speed of the aircraft has been xed at 500 kts). Those ights have been spread over 6 hours and the number of slots has been xed at 180. This induces a slot duration of 2 minutes. For each Origin-Destination, a set of possible dierent routes is randomly built in a way that never induces an extradistance greater than 30 percent of the original distance with a maximum of 10 routes. For each ight the maximum slot moving has been xed to 15 in the future and 15 in the past from the original slot of departure. Two kinds of experiments have been done to compare the gain given by the delay allocation only (the departure slot can only be moved in the future (not in the past) and the original route is used) and the time-route allocation (dierent routes can be used and the moving of the slot of departure can be done in the past or in the future). The rst is a kind of Ground Holding Policy(GHP) and the second one is a Full Planning Policy(FPP). The original plannings induce a control workload in the airspace with a maximum for each sector at a given slot. The workload weight coecients used for these experiment are : = 1:0; = 1:0; = 1:0. The half extension of the temporal window used to compute the ows from the discrete events has been xed at D = 2 and induces a window duration of 10 minutes. The maximum control workload registered in each sector over the six hour duration is given in gure 12 (with this data set). On this gure, sector indices are indicated on the horizontal axis and the associated maximum sector workload on the vertical axis. One can notice that sector 9 is the most overloaded. For this problem, with the de nition of the tness, the genetic algorithm is expected to reduce those peaks of congestion for the most congested sectors.

6

12

8 15

3

13 14

4 9

10 11

2 7

5 1

Figure 11: Structure of the test network Maximum sector workloads over the day 2000

1500

1000

500

1

2

3

4

5

6

7

8

9

10 11 12 13 14 15

Figure 12: Max workload distributions over sectors

Fitness evolution (ALEAT100) 2.5

Bestfit

2

1.5

1

5

10

15

20

5

10

15

20

25

30

35

40

45

50

25 30 generations

35

40

45

50

2.5

Avgfit

2

1.5

1

Figure 13: Best and average tness evolution over generations

6.2 Results The results presented in this section are obtained from runs on a SPARC5 workstation. The large size of the chromosomes (6000 bytes) limits the maximum size of the population used by the genetic algorithm at 90 due to memory limitations (a population of 50 individuals has been used for all the experiments). It must be noticed that this small population size will reduce the possibilities of genetic algorithms, and better results could be reached on a machine with more memory (even if the following results are already very good). The others parameters (probability of crossover, probability of mutation, number of generations) have been xed by experimentation at Pc = 0:6, Pm = 0:2 and NB gen=100. The technical features of the genetic algorithm were the following : Selection : Stochastic remainder without replacement ; Crossover : Crossover with a Simulated Annealing Tournament ; Scaling : Power Law Scaling ; Elitist Strategy Applied ; No scharing. Given all these parameters, each run lasts about 15 minutes on a SPARC5 workstation. This running time is not critical because the problem solved has not any real time constraints (the improved time-route planning could be applied for several months). The associated tness evolution is given in gure 13. On this gure, the two kinds of policies (GHP and FPP) can be identi ed by the dierent lines used (plain for FPP and dashed for GHP). As expected, the FPP policy gives better results than the GHP policy with respectively a best tness of 1.79(FPP) and 1.45(GHP). These results are relevant due to the fact that the FPP policy has more degrees of freedom. It means that original \congestion" in the sense given by the optimization criterion, could be respectively divided by 1.79 and 1.45 with the new ight plans. Even with the small population size used, the results given by the genetic algorithm are very encouraging. The associated max workload distributions over sectors is given on gure 14. On this gure, it can be noticed that the max workload on the most overloaded sector (number 9) has been 1841 divided by 1.673 ( 1100 ) for the GHP policy and by 2.77 ( 1841 ) for the FPP policy. A more detailed description of 663 the sector control workload evolutions is given in gures 15,16,17 (on these diagrams workload has been divided by 1000 for presentation). As expected, the workload is spread around the peak as in a smoothing process. The previous curves make us notice that the algorithm only aects the most congested sectors without modifying the less loaded ones. This is due to the de nition of the criterion which takes into account only the most loaded sectors. So, when the recombination operator is applied, if they aect aircraft involved in the underloaded sectors (or uncongested airports) the induced tness will not change a lot and this makes the crossover and mutation irrelivant. Then, those operators have to be changed in a way that aect only aircraft involved in the most loaded sectors (or in the most congested airports) during the peak periods. This is the subject of the following section.

Maximum sector workloads over the day 2000

1500

1000

500

1

2

3

4

5

6

7

8

9

10 11 12 13 14 15

Figure 14: Max workload over sectors for the two associated policies (-. GHP ... FPP)

W1

2 1

W2

0 0 2

W3

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

1 0 0 2

W4

40

1 0 0 2

1 0 0 2

W5

20

1 0 0

Figure 15: Sector control workload evolution over the day (-. GHP ... FPP) W6

2 1

W7

0 0 2

W8

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

1 0 0 2

W9

40

1 0 0 2

1 0 0 2

W10

20

1 0 0

Figure 16: Sector control workload evolution over the day (-. GHP ... FPP)

W11

2 1

W12

0 0 2

W13

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

1 0 0 2

W14

40

1 0 0 2

1 0 0 2

W15

20

1 0 0

Figure 17: Sector control workload evolution over the day (-. GHP ... FPP)

7 Improvement of the recombination operators 7.1 Introduction

To be able to recognize the aircraft involved in the biggest sector congestion peak or in the biggest airport congestion peak, new information must be added to the chromosome which indicates for each gene, the maximum level of sector congestion encountered during a ight and the level of congestion at the destination airport when the aircraft land. To reach this goal, the biggest peak of congestion is determined over all the sectors and airports, and a number of congestion levels is de ned respectively for the airports and the sectors. With these two parameters, and for each slot, a level of congestion can be associated with each sector and each airport. So, after each evaluation, the associated list of sectors of each individual simpli ed ight plan is read in a way to identify the maximum level of congestion encountered during the ight and at the destination airport. Those congestion levels will be used to improve the recombination operators.

7.2 New Crossover Operator The successive steps of this new crossover operator are the following :

two parents are rst selected according to their tness ; the summation of the congestion levels (sector and airport) is computed for each ight in bothp parents. For p p p

a ight n, total congestion level in the parent p will be noted Wn =WSn + WAn (with WSn : maximum level of sector congestion encountered by the ight n in the parent p, WApn : airport congestion level encountered by the ight n when arriving at destination (in the parent p)) ; ; an order relationship is then constructed with the total congestion level in the following way :

{ ight planing n in parent 1 is said to be \much better" than ight planing n in parent 2 if Wn < :Wn ; where 2 [0:7; 0:95]; 1

2

Fitness evolution (ALEAT100) 2.5

Bestfit

2

1.5

1

5

10

15

20

5

10

15

20

25

30

35

40

45

50

25 30 generations

35

40

45

50

2.5

Avgfit

2

1.5

1

Figure 18: New Best and average tness evalution over generations

{ ight planing n in parent 2 is said to be \much better" than ight planing n in parent 1 if Wn < :Wn ; { ight planing n in parent 1 and in parent 2 are said to be \equivalent" if none of the previous relations 2

1

matches; if a ight planning \is much better" in the rst parent than in the second then it is copied in the second ; if a ight planning \is much better" in the second parent than in the rst then it is copied in the second ; if the two ight plannings \are equivalent" they are randomly exchanged with a constant probability (0.5) ; This is done for all the ights in the chromosome.

7.3 New Mutation Operator As already noted, this new operator must only aect the ights involved in the highest peaks of congestion. This new operator works in the following way :

two threshold congestion levels are randomly chosen : one for the sector congestion (ThS ) and one for the

airport congestion (ThA) ; then for each ight n in the chromosome the following are applied : if (WAn > ThA) or(WSn > ThS ) then the associated ight plan is modi ed by randomly choosing either a new route or a new slot moving (but not both) in the associated data set (Rn ; n ) ; else the ight planing is unchanged;

7.4 New Results The problem described in section 6 has been resolved using this new genetic algorithm on the same machine (SPARC5) with the same parameters (population size, probability of crossing etc...). The same policies have been compared (GHP and FPP) one the same random trac sample. These new operators greatly improve the tness evolution as it can be seen in gure 18. The associated new best tnesses are 2.32 for the FPP and 1.68 for the GHP, as compared with the results given by the classical recombination operators (FPP : 1.79 and GHP : 1.45). The associated max workload distributions over sectors are given on gure 19. The improvement on the most overloaded sector (number 9) is now 2.73 with the GHP ( 1841 ) and 6.71 with 672 ). The associated time evolution of the sector workload is given in gures 20,21,22. the FPP( 1841 274

Maximum sector workloads over the day 2000

1500

1000

500

1

2

3

4

5

6

7

8

9

10 11 12 13 14 15

Figure 19: New max workload over sectors for the two associated policies (-. GHP ... FPP)

W1

2 1

W2

0 0 2

W3

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

1 0 0 2

W4

40

1 0 0 2

1 0 0 2

W5

20

1 0 0

Figure 20: Sector control workload evolution over the day (-. GHP ... FPP) W6

2 1

W7

0 0 2

W8

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

1 0 0 2

W9

40

1 0 0 2

1 0 0 2

W10

20

1 0 0

Figure 21: Sector control workload evolution over the day (-. GHP ... FPP)

W11

2 1

W12

0 0 2

W13

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

20

40

60

80

100

120

140

160

180

1 0 0 2

W14

40

1 0 0 2

1 0 0 2

W15

20

1 0 0

Figure 22: Sector control workload evolution over the day (-. GHP ... FPP)

8 Conclusion This dicult problem has been handled well by the proposed genetic algorithm even with a small population (50 individuals) and a small number of generations (100). Even better results could be expected with larger populations. A mathematical model has been developed to formulate this bi-allocation (route and slot of departure) problem in a way that could be adapted to the real ight plans in three dimensional airspace. So, the next step, is to try this optimization technique on a real trac sample. The strength of this model is its ability to manage the constraint of the airline companies in a microscopic way by using individual sets of decision variables associated with each ight. Using those individual decision variables, macroscopic models of the sector control workload and of airport congestion have been developed by computing the associated ows with a time moving window on the strategic points (crossings, sector boundaries, airports). This reduces the sensitivity of the \optimum" given by the genetic algorithm with the proposed solution. This model can also take into account the connexion between

ights (hub phenomenon) by changing the decisions variable sets of the following ight. This model can also use dynamic capacities for both airports and sectors. A fast simulation process has been used to compute the associated sector workload and the airport congestion to improve the performance of the genetic algorithms. A realistic example has been treated and two policies have been tested : GHP (delay decision only) and FPP (route and slot of departure decision). As expected the second one gives much better results due to the larger state domain induced by the two extra degrees of freedom of the FPP policy. Two kinds of genetic algorithms have been tested : the rst one has a classical crossover(slicing) and mutation( random gene) operators and the second one focuses its recombination actions on the most congested sectors by assigning two levels of congestion for each ight(one for the most congested sector encountered during the ight and one for the congestion of the destination airport when the ight is arriving). This new information greatly improves the quality of the given solutions. It must be noticed, that only aircraft involved in the most congested parts of the airspace (sector or airport) receive a ight plan modi cation. A FPL change cost could be added to classify aircraft and identi ed the one which must be modi ed rst.

References [1] A.H.L Aarts and P.J.M Van Laarhoven. Statistical cooling : A general approach to combinatorial optimization. Philips Journal of Research, 40:193{226, 1985. [2] E.H.L Aarts and J Korst. Simulated Annealing and Boltzmann Machines : A Stochastic Approach to Combinatorial Optimization and Neural Computing. John Willey and Sons, New York, 1989. [3] A.S Alfa. Departure rate and route assignment of commuter trac during peak period. Transportation Research, 23B(5):337{344, 1989. [4] G Andreatta, A.R Odoni, and O Richetta. Models for the ground holding problem. In L Bianco and A.R Odoni, editors, Large Scale Computation and Information Processing in Air Trac Control, Transportation Analysis, pages 125{168. Springer-Verlag, 1993. [5] G Andreatta and G Romanin-Jacur. Aircraft ow management under congestion. Transportation Science, 21(4), 1987. [6] S Anily and A Federgruen. Simulated annealing methods with general acceptance probabilities. Advances in Applied Probability, 24:657{667, 1987. [7] R Arnott, A DePalma, and R Lindsey. Departure time and route choice for the morning commute. Transportation Research, 24B(3):209{228, 1990. [8] M.G.H Bell. Stochastic user equilibrium assignment in networks with queues. Transportation Research, 29B(2):125{137, 1995. [9] M Ben-Akiva. Dynamic network equilibrium research. Transportation Research A, 19A(5):429{431, 1985. [10] M Ben-Akiva, A DePalma, and P Kanaroglou. Dynamic model of peak period trac congestion with elastic arrival rates. Transportation Science, 20(2):164{181, 1986. [11] M Ben-Akiva, A DePalma, and I Kaysi. Dynamic network models and driver information systems. Transportation Research, 25A(5):251{266, 1991. [12] M Ben-Akiva, H.N Koutsopoulos, and A Mukundan. A dynamic trac model system for atms/atis operation. IVHS Journal, 2(1):1{9, 1994. [13] D Bernstein, T.L Friesz, R.L Tobin, and B.W Wie. A variational control formulation of the simultaneous route and departure-time choice equilibrium problem. In Transportation and Trac Theory, 1993. [14] D.J Bertsimas and S Stock. The air trac ow management problem with en-route capacities. Technical report, A.P Sloan School of Management. M.I.T, 1994. [15] Wie B.W., R.L Tobin, T.L Friesz, and D Bernstein. A discrete time, nested cost operator approach to the dynamic network user equilibrium problem. Transportation Science, 29(1):79{92, 1995. [16] G.E Cantarella and E Cascetta. Dynamic processes and equilibrium in transportation networks : Toward a unifying theory. Transportation Science, 29(4):305{328, 1995. [17] M Carey. Optimal time-varying ows on congested networks. Operations Research, 35(1):58{69, 1987. [18] E Cascetta. Static and dynamic models of stochastic assignment to transportation networks. In Flow Control of Congested Networks, 1987. [19] E Cascetta. A stochastic process approach to the analysis of temporal dynamics in transportation networks. Transportation Research, 23B(1):1{17, 1989. [20] E Cascetta and G.E Cantarella. A day-to-day and within-day dynamic stochastic assignment model. Transportation Research, 25A(5):277{291, 1991. [21] E Cascetta and G.E Cantarella. Modeling dynamics in transportation networks : State of the art and future developments. In Simulation Practice and Theory 1. SIMPRA, 1993.

[22] G.L Chang and H.S Mahmassani. Travel time prediction and departure time ajustment behavior dynamics in a congested trac system. Transportation Research, 22B(3):217{232, 1988. [23] G.L Chang and H.S Mahmassani. Travel time prediction and departure time ajustment behavior dynamics in a congested trac system. Transportation Research, 22B(3):217{232, 1988. [24] S Dafermos. Trac equilibrium an variational inequalities. Transportation Science, 14(1):42{54, 1980. [25] S Dafermos. Relaxation algorithms for the general asymmetric trac equilibrium problem. Transportation Science, 16(2):231{239, 1982. [26] S Dafermos and F.T Sparrow. The trac assignment problem for a general network. Journal of Research of the National Bureau of Standards, 73B:91{118, 1969. [27] D Delahaye, J.M Alliot, M Schoenauer, and J.L Farges. Genetic algorithms for air trac assignment. In Proceedings of the European Conference on Arti cial Intelligence. ECAI, 1994. [28] D Delahaye, J.M Alliot, M Schoenauer, and J.L Farges. Genetic algorithms for partitioning airspace. In Proceedings of the Tenth IEEE Conference on Arti cial Intelligence for Application. CAIA, 1994. [29] D Delahaye, J.M Alliot, M Schoenauer, and J.L Farges. Genetic algorithms automatic regrouping of air trac sectors. In Proceedings of the Fourth International Conference on Evolutionary Programming. Natural Selection inc., 1995. [30] A DePalma, P Hansen, and M Labbe. Commuters'paths with penalties for early or late arrival time. Transportation Science, 24(4):276{286, 1990. [31] O Drissi-Kaitouni and A Hameda-Benchekroun. A dynamic trac assignment model and a solution algorithm. Transportation Science, 26(2):119{128, 1992. [32] N Durand and Alliot. Automatic aircraft con ict resolution using genetic algorithms. In Proceedings of the 11th Annual ACM Conference on Applied Computing, Philadelphia, 1996. (ACM/SAC. [33] J.E Fernandez and T.L Friesz. Equilibrium predictions in transportation markets : the state of the ar. Transportation Research B, 17B(2):155{172, 1983. [34] C Fisk and S Nguyen. Solution algorithms for network equilibrium models with asymmetric user costs. Transportation Science, 16(3):361{381, 1982. [35] M Florian and H Spiess. The convergence of diagonalization algorithms for asymmetric network equilibrium problems. Transportation Research B, 16B(6):477{483, 1982. [36] D.B Fogel. Evolutionary Computation. Toward a new Philosophy of Machine Intelligence. IEEE press, 1994. [37] L.J Fogel, A.J Owens, and M.J Walsh. Arti cial Intelligence Through Simulated Evolution. Wiley and sons. NY, 1966. [38] T.L Friesz. Transportation network equilibrium, design and aggregation : Key developements and research opportunities. Transportation Research, 19A(5):413{427, 1985. [39] T.L Friesz, D Bernstein, T.E Smith, and B.W Wie. A variational inequality formulation of the dynamic network user equilibrium problem. Operations Research, 41(1):179{191, 1993. [40] T.L Friesz, J Luque, R.L Tobin, and B.W Wie. Dynamic network trac assignment considered as a continuous time optimal control problem. Operation Research, 37(6):893{901, 1989. [41] N.H Gartner. Optimal trac assignment with elastic demands : A review. part i : Analysis framework. Transportation Science, 14(2):174{191, 1980. [42] N.H Gartner. Optimal trac assignment with elastic demands : A review. part ii : Algorithmic approaches. Transportation Science, 14(2):191{208, 1980. [43] N.H Gartner, S.B Gershwin, J.D.C Little, and P Ross. Pilot study of computer based urban trac management. Transportation Research B, 14B(1):203{217, 1980.

[44] M.O Ghali and M.J Smith. A model for the dynamic system optimum trac assignment. Transportation Research, 29B(3):155{170, 1995. [45] E.P Gilbo. Airport capacity : Representation, estimation, optimization. IEEE Transaction on Control Systems Technology, 1(3):144{154, 1993. [46] D.E Goldberg. Genetic Algorithms in Search, Optimization and Machine Learning. Reading MA Addison Wesley, 1989. [47] D.W Hearn. A dual ascent algorithm for trac assignment problems. Transportation Research B, 24B(6):423{ 430, 1990. [48] D.W Hearn, S Lawphongpanich, and S Nguyen. Pilot study of computer based urban trac management. Transportation Research B, 18B(4):357{365, 1984. [49] M.P Helme. Reducing air trac delay in a space-time network. IEEE, 1992. [50] C Hendrickson and G Kogur. Schedule delay and departure time decisions in a deterministic model. Transportation Science, 15(1):62{77, 1981. [51] J.K Ho. Solving the dynamic trac assignment problem on a hypercube multicomputer. Transportation Research, 24B(6):443{451, 1990. [52] J.H Holland. Adaptation in Natural and Arti cial Systems. University of Michigan press, 1975. [53] B.N Janson and J Robles. Dynamic trac assignment with arrival time cost. In Transportation and Trac Theory, 1993. [54] R Jayakrishnan, K.W.K Tsai, and A Chen. A dynamic trac assignment model with trac- ow relationship. Transportation Research, 3C(1):51{72, 1995. [55] R.C Jou and H.S Mahmassani. Comparability and transferability of commuter behavior chracteristics between cities : Departure time and route switching decisions. In Transportation Research Record (73 Annual Meeting, 1993. [56] D.E Kaufman, J Nonis, and R.L Smith. A mixed integer linear programming formulation of the dynamic trac assignment problem. IEEE, 1992. [57] S Kirkpatrick, C.D Gelatt, and M.P Vecchi. Optimization by simulated annealing. Science, 220(4598):671{ 680, 1983. [58] J.R Koza. Genetic Programming. MIT press, 1992. [59] S Lawphongpanich and D.W Hearn. Simplicial decomposition of the assignment problem. Transportation Research B, 18B(2):123{133, 1984. [60] L.J Leblanc, B Ran, and D.E Boyce. Dynamic travel choice models for urban transportation networks. IEEE, 1992. [61] M Lundy and A Mees. Convergence of an annealing algorithm. Mathematical Programming, 34:111{124, 1986. [62] H Mahmassani and R Herman. Dynamic user equilibrium departure time and route choice on idealized trac arterials. Transportation Science, 18(4):362{384, 1984. [63] H. Mahmassani and S Peeta. System optimal dynamic assignment for electronic route guidance in a congested trac network. In Proceeding of the Second International Capri Seminar on Urban Trac Networks, 1992. [64] H.S Mahmassani and G.L Chang. Experiments with departure time choice dynamics of urban commuters. Transportation Research, 20B(4):297{320, 1986. [65] H.S Mahmassani and R Jayakrishnan. System performance and user response under real-time information in a congested trac corridor. Transportation Research, 25A(5):293{307, 1991.

[66] F.L Mannering, S.A Abu-eisheh, and A.T Arnadottir. Dynamic trac equilibrium with discrete/continuous econometric models. Transportation Science, 24(2):105{116, 1990. [67] P Marcotte. Network optimization with continuous control science. Transportation Science, 17(21):181{197, 1983. [68] L Maugis. Mathematical programming for the air trac ow management problem with en-route capacities. IFOR, 1996. [69] D.K Merchant and G.L Nemhauser. A model and an algorithm for the dynamic trac assignment problems. Transportation Science, 12(3):183{199, 1978. [70] D.K Merchant and G.L Nemhauser. Optimality conditions for a dynamic trac assignment model. Transportation Science, 12(3):200{207, 1978. [71] N Metropolis, M Rosenbluth, M Rosenbluth, A Teller, and E Teller. Equation of state calculations for fast computing machines. Journal of Chemical Physics, 6:1087{1092, 1953. [72] Z Michalewicz. Genetic algorithms + Data Structures = Evolution Programs. Springer-verlag, 1992. [73] D Mitra, F Romeo, and A.L Sangiovanni-Vincentelli. Convergence and nite-time behavior od simulated annealing. Advances in Applied Probability, 18:747{771, 1986. [74] A Nagurney. An equilibrium scheme for the trac assignment problem with elastic demands. Transportation Research B, 22B(1):73{79, 1988. [75] A Nagurney. NETWORK ECONOMICS A Variational Inequality Approach, volume 1 of Advances in Computational Economics. Kluwer Academic, 1993. [76] A.R Odoni. The ow management problem in air trac control. In A.R Odoni et al, editor, Flow Control of Cogested Networks, volume F38 of ASI Series, pages 269{288. NATO, 1987. [77] M Papageorgiou. Dynamic modeling, assignment, and route guidance in trac networks. Transportation Research, 24B(6):471{495, 1990. [78] M. Papageorgiou. Concise encyclopedia of trac and transportation systems. Pergamon Press, 1991. [79] C.H Papdimitriou and K Steiglitz. Combinatorial Optimization : Algorithms and Complexity. Prentice-Hall, New York, 1982. [80] W.B Powell and Y She. The convergence of equilibrium algorithms with predetermined step size. Transportation Science, 16(1):45{55, 1982. [81] B Ran, D.E Boyce, and L.J LeBlanc. A new class of instantaneous dynamic user-optimal trac assignment models. Operations Research, 41(1):192{202, 1993. [82] O Richetta and A.R Odoni. Solving optimally the static ground holding policy problem in air trac control. Transportation Science, 27(3):228{238, 1993. [83] O Richetta and A.R Odoni. Dynamic solution to the ground-holding problem in air trac control. Transportation Research, 28A(3):167{185, 1994. [84] H.P Schwefel. Evolution and Optimum Seeking. Wiley, New York, 1995. [85] M.J Smith. An algorithm for solving asymmetric equilibrium problem with a continuous cost- ow function. Transportation Research B, 17B(5):365{371, 1983. [86] M.J Smith. The existence and calculation of trac equilibria. Transportation Research B, 17B(4):291{303, 1983. [87] M Terrab and A.R Odoni. Strategic ow management for air trac control. Operations Research, 41(1):138{ 152, 1993.

[88] D Van Vliet. The frank-wolfe algorithm for equilibrium trac assignment viewed as a variational inequality. Transportation Research B, 21B(1):87{89, 1987. [89] D.V Vliet and P.D.C Dow. Capacity-restrainned road assignment. Trac Engineering and Control, pages 296{305, 1979. [90] P Vranas, D Bertsimas, and A.R Odoni. The multi-airport ground-holding problem in air trac control. Operation Research, 42(2):249{261, 1994. [91] P.B.M Vranas, D Bertsimas, and A.R Odoni. Dynamic ground-holding policies for a network of airports. Transportation Science, 28(4):275{291, 1994. [92] P.C Vythoulkas. A dynamic stochastic assignment model for the analysis of general networks. Transportation Research, 24B(6):453{469, 1990. [93] J.G Wardrop. Some theoritical aspects of road trac research. In Proceedings of the Institute of Civil Engineers. Part II, pages 325{378, 1952. [94] B.W Wie, T.L Friesz, and R.L Tobin. Dynamic user optimal trac assignment on congested multidestination networks. Transportation Research, 24B(6):431{442, 1990. [95] B.W Wie, R.L Tobin, and T.L Friesz. The augmented lagrangian method for solving dynamic network trac assignment models in discrete time. Transportation Science, 28(3):204{220, 1994. [96] D.J Zawack and G.L Thompson. A dynamic space-time network ow model for city trac congestion. Transportation Science, 21(3):153{162, 1987. [97] S Zenios. Network based models for air trac control. European Journal of Operational Research, 50:166{178, 1991.