- Email: [email protected]
Network based models for air-traffic control
166
European Journal of OperationalResearch50 (1991) 166-178 North-Holland
Theory and Methodology
Network based models for air-traffic control Stavros A. Zenios Department of Decision Sciences, The Wharton School, University of Pennsylvania, Philadelphia, PA 19104, USA
Abstract: Several aspects of air-traffic control fit into the framework of network optimization models with nonlinear and nonseparable objectives. We develop two network based models for air-traffic control over the U.S. airspace. The first model deals with the dual objectives of monitoring congestion over high altitude jet routes while minimizing transportation cost. The second deals with the optimal flow of departing traffic through major airports to avoid delays at the destination. Of particular interest is the first model, due to the high increase of traffic volume over the last few years and the resulting congestion and increased number of accidents. A prototype model is built using data from a control sector in the Midwest. By integrating recent advances in numerical optimization, network flow theory and supercomputing we were able to solve the resulting problem very efficiently, thus demonstrating its feasibility. Computational experiments with real data on a CRAY X-MP/24 supercomputer are discussed.
Keywords: Traffic management, network optimization, supercomputing
1. Introduction
On an average day in 1984 approximately 17 000 aircraft traveled the high altitude jet routes over the continental United States. Takeoff and landing of one such plane are only two operations among millions that the Air Traffic Control (ATC) system handles each year. Some 14 000 controllers monitor the flow of traffic on these routes working in three basic facilities: 398 airport towers, 18 TRACON (Terminal Radar Approach Control) and 20 en route control centers. The complexities of the ATC are analyzed in a series of articles in a special issue of IEEE Spectrum (1986).
ReceivedNovember1987; revised May 1989
Air-traffic controllers strive to prevent collisions and to keep the traffic moving efficiently. They are responsible for keeping each plane in their jurisdiction separated by at least 1 000 feet vertically and by at least 5 miles horizontally. Remarkably, despite the potentials for failure, most airline flights are uneventful. However, traffic volume has been increasing steadily and with it have increased the risk of accidents and potential inefficiences. The number of reported near-collisions in mid-air increased from 584 in 1984 to 758 in 1985. It is expected that traffic volume will increase by 62% over the next decade. In response, the Federal Aviation Administration (FAA) plans to spend $11 billion on a program to update the air-traffic system - called the National Airspace System Plan. One of the primary objectives of this plan is to develop suitable procedures for assessing, and ultimately controlling, the efficiency and
0377-2217/91/$03.50 © 1991 - ElsevierSciencePublishersB.V.(North-Holland)
S.A. Zenios / Network based modelsfor air-traffic control
safety of the system. In this paper we propose two network based models for the U.S. air-traffic control. The first model deals with the dual objectives of monitoring congestion over the high altitude jet routes while minimizing transportation costs. The second model deals with the optimal flow of departing flights through major airports to avoid delays at the destination. The use of mathematical programming tools for studying problems in civil aviation originated with Ferguson and Dantzig (1957). More recently Bielli et al. (1982) developed a linear multicommodity network model for traffic flow near the Rome airport. Booth and Harvey (1983) proposed an integer programming formulation for minimizing traffic delays. Friedman (1986) discussed a pacing algorithm for determining the optimal velocity of aircrafts for meeting a prescribed schedule. The aforementioned models brought insight to the flight control problem and were well received by practitioners - refer to the proceedings of the Princeton Conference (1983). However, although the plausibility of these models is evident, they are far from being implementable. First, they do not consider in an integrated fashion the risk associated with a particular economic-efficient policy. Second, no attempt was made to evaluate these models using real data from large scale instances or to develop computationally tractable algorithms for their analysis. The work described here deviates from earlier research in several factors: (1) it integrates in a unified framework the congestion (as a risk surrogate) and cost components of air-traffic control; (2) emphasizes computationally efficient solution of these models by integrating advances in numerical optimization, network flow theory and supercomputing; and (3) evaluates the applicability of these models using real data from a representative sector of the Indianapolis control center over a six-hour time interval. The paper is organized as follows: Section 2 gives details of the air-traffic control models describing the network structure. The development of a prototype for the Indianapolis control sector is discussed in Section 3, followed by a presentation of the solution methodology (i.e. numerical optimization and supercomputing) in Section 4. Computational results are presented in Section 5.
167
The final section discusses possible uses of the model, identifies directions for further research and summarizes the conclusions.
2. Air-traffic control modeling Airspace, for the purpose of air-traffic control, is divided in a number of ways. At heights above 18000 feet all airspace is considered part of the ATC system. This controlled airspace is subdivided horizontally and vertically into some 500 sectors, each one handled by an air-traffic controller and an assistant. As a plane passes from one sector to another it is handled off to a different controller. At the lower altitudes the airspace over busy terminals is designated a terminal control area. All aircraft entering this space - shaped like an inverted cone - must be in contact with a controller. Below the sector airspace and around the terminal control areas general aviation aircraft operates under visual flight rules. In addition to monitoring traffic through individual sectors by the controllers, a computer network calculates the nationwide flow of traffic into some 40 key airports. If an airport is becoming clogged and more airplanes can arrive than can be accomodated, a schedule of ground delays is computed to keep planes from taking off for the congested airports. The overall flow of traffic is therefore determined by possible delays introduced by the centralized system or by alternative routing decisions made by the controllers. We describe now mathematical modeling approaches to these two components of the system. The first model looks at the problem from the controller's point of view, determining optimal altitudes and possible delays to avoid congestion. The second model studies the global problem of smooth traffic flow, prescribing optimal delays at key airports. 2.1. Integrated cost / congestion planning model
Figure 1 depicts a multi-period model for a single sector over all possible altitudes. Nodes indicate the origin or destination of flights whose paths travel through the target sector during time intervals of interest (i.e. the time horizon). Arcs indicate interconnecting paths between origin and destination through the target sector. We do not require the origin or destination airports of a
168
S.A. Zenios / Network based models for air-traffic control
Fig. 1. Integrated cost/congestion network model
flight to be in the jurisdiction of the target sector; the flight paths however should cross the sector for a flight to be included in this model. Note that an aircraft can follow a number of alternative altitudes in traveling through a sector on a predetermined path. As we extend the model for multiple time periods we introduce additional nodes in the obvious manner: one node for each airport during every time period. Flow on arcs interconnecting the multiple time periods represent delays - either at the origin or destination airports. A collector node for every destination airport is introduced to monitor terminal facilities. The flow on the arcs entering a collector node during a particular time period indicates total arrivals and is limited due to airway capacities. The use of multipliers on these arcs (generalized network model) changes units from planes to passengers and can be used to study the utilization of airport facilities. The underlying optimization model includes the dual objectives of optimizing congestion and cost. Cost is measured in terms of fuel consumption, either en route or due to delays at takeoff or landing. Other costs, like crew salaries and a surrogate for customer dissatisfaction, can also be included. Quantifying the risk component of the system is a critical and highly controversial issue. For example, Odoni and Endoh (1983) consider a probabilistic measure of risk at the micro-event
level and McDonald (1985) measures risk as a function of controllers' workload. For the purpose of our analysis we are optimizing a relative risk measure: the key idea is to compare some risk norm under different scenaria or system states. It is recognized that this measure of risk does not translate into a probabilistic notion of collisions. However, it provides a measure of a preferred quantity (less congestion is preferable than more congestion). As such it can assist in identifying an efficient state of a control sector with low congestion. Congestion, as specified through appropriate occupancy coefficients, is employed as a risk surrogate. The coefficients are defined as follows: Lt Ot(i,j),(k,I)
= (Time spent by aircrafts on route (i, j ) on same altitude L as aircrafts on route (k, 1) during period t) / ( T i m e spent by aircrafts on route (i, j ) on altitude L during period t). Occupancy coefficients indicate the fraction of time that two planes are in the vicinity of each other thus contributing to the congestion of the sector. See Section 3 for a sample calculation of the coefficients.
S.A. Zenios / Network based models for air-traffic control
The system-wide relative risk is then defined as
R(x) = ~, ~., ~ t
L,
1
Lt Lt~ 2 gt(i,jXk,l)~XijXkl ) '
L (i,j)(k,l)
where x ,L, and x[5' indicate the number of flights on paths (i, j ) and (k, l) at altitude L and time period t. We introduce some additional notation in order to formalize the model: i, j,
= Origin airport node at time period t. = Destination airport node at time period
J
= Collector node (sink) for destination airport j, for all t. = N u m b e r of aircrafts staying at origin node i during time period t to t + 1. This is a convenient notation for x66+. = N u m b e r of aircrafts delayed at destination node j during period t to t + 1. This is a convenient notation for xj,j,+. = N u m b e r of aircrafts from airport i to airport j during period t at flight altitude L. This is a convenient notation for x ttj (~(t" = N u m b e r of flights at airport j granted permission to land at time period t. This is a convenient notation for x j, j. = Fuel consumption of flights delayed on the ground at origin airport i. = Fuel consumption of flights delayed airborne at destination airport j. = Fuel consumption of flights from i to j traveling at altitude L. = N u m b e r of flights originating at airport i during period t. = Bound on the number of flights that airport i can accomodate in takeoff queue. = Bound on the number of flights that airport j can accomodate in landing queue. = Bound on the number of passengers that airport j can accomodate. = Average number of passengers per plane landing at airport j during period t.
t.
x~
x~
x~ '
x~j
CG, CAj L C&, s6 UT,
ULj
dJ
m',.,
The subscript t of a parameter or variable is dropped if there is no ambiguity in doing so. For example we write CG~ instead of CGi, since this
169
parameter is time-independent; likewise for CA j, U T i, U L j and dj. The time indices t are used to indicate a time slice of fixed duration - in our empirical investigation it is set equal to approximately ½ hour. During this interval we consider all the flights through the sector both for congestion and fuel consumption via the variables xij. Lt Of course this assumes that all the planes through the target sector will take the same amount of time to travel through. This is a simplification that assumes a circular shape for the sector and identical speeds for all the aircraft. While the first assumption is not unreasonable, the second one is in general not true. Hence we use as a time interval the m a x i m u m time tmax it takes for an aircraft in the model to go through the sector. Any aircraft that travels through the sector is less than tmaX is still assumed to be in the sector for the duration tmax. Hence the risk measure is an upper bound on the true congestion. A limitation with this choice of time interval is that airplanes are assumed to arrive at their destination at the same time (since they spent identical times to travel through the sector). As a result the impact of the schedule on the destination is not known and the model is of little value there. The purpose of the model is primarily to optimize flow of traffic and congestion at the target control sector. The impact of this model on congestion at the origin and destination airport can be limited by allowing only a small number of short delays in our prototype up to three 10-minute intervals. Hence the developed schedule will be close to a prototype macro-schedule developed to optimize congestion at airports (see the model of Section 2.2). The more flexibility is given to (OPT1) the less congested this sector will be at the expense of long delays and congestion at the origin and destination airports. With this discussion on notation and congestion function we may now formalize the model: (OPT1) minimize W1
{ ~ ~ CGix[ + ~2 ~ " t
it
t
CAx;- + Y'. ~
j,
+ w2R(x) subject to t--1
E E x L t "-~ X~ = Si-]- Xi L j,
L
~
t (i,j)
Ci9x[)t}
S.A. Zenios / Network based raodels for air-traffic control
170
for all origin airports i and periods t,
(1)
~f~xLt ~,.., ,j +xjt - 1 = x jtj + x j t L
it
for all destination airports j and periods t, (2) m jtj x j jt <~dj
for all destination airports j and periods t,
(3)
0 ~
(4)
0 ~
(5)
for all destination airports j,
W 1 "{" W2 = 1 .
(6)
The first two terms of the objective represent fuel consumption either due to delays on the ground before take off or delayed landing at the destination airport. The third term accounts for transportation costs for flights over all possible altitudes. The last term in the objective reflects the risk measure as discussed earlier. Constraints (1) are flow continuity constraints at the origin airports. All flights scheduled for departure from an airport are allowed to do so with some (possible) delay as indicated by the variable xj. Constraints (2) impose similar conditions at the destination airports. Constraints (3) account for all passengers destined for airport J and it can handle limits on airport facilities. If all multipliers are set equal to one (m~s = 1), then the model accounts for aircrafts and not passengers. This set of constraints is useful in monitoring passenger service facilities (e.g. limits imposed by waiting room capacity, baggage claim areas, customs inspection queues and so on). Constraints (4) and (5) impose a limit on the number of aircrafts a particular airport can handle in the takeoff or landing queues respectively. The parameters wt and w2 are goal programruing weights used in the context of the multicriteria optimization model developed. They are not meant to be set a priori. Instead they are used as part of a sensitivity analysis to identify an efficient frontier for the operation of the system. This issue is discussed further in a later section.
2.2. Air-traffic flow control Model (OPT1) considers flights with paths through a control sector of interest. It can be used as a real time operational model to determine the number of flights allocated to the different al-
titudes of the target sector given a schedule of flights over a short time period (for example a few hours). It also estimates (short) delays for different flights to optimize the operation of the target sector. We study now a global model that considers the flow of traffic through major airports at a nationwide level. The objective of this model is to minimize delays of departing flight while avoiding congestion at the destination airports. This target is a achieved by imposing bounds on the throughput capacity of both origin and destination airports. The optimization is achieved on a system-wide basis without any consideration of congestion at intermediate sectors. This problem is currently solved by simulation as discussed in the special issue of Spectrum (1986). As a result of this model control sectors may become highly congested when flights from different origin/destination pairs are going through. Model (OPT1) would be used to optimize the control sectors, creating of course more congestion at the airports than the congestion dictated by the model we introduce next. Integrating the two models successfully still remains an open question. A simple network model studying four airports is depicted in Figure 2. Every airport (represented by triangular nodes) has a set of inbound and a set of outgoing arcs indicating arriving and departing flights. Every outgoing arc corresponds to a particular time period - flow on these arcs indicates departing flights during every time interval. Incoming arcs represent the arrival time of flights. A node is used to represent departing flights for every time period and a second node represents arriving flights. A third set of nodes indicate the destination (or origin) airports for all flights during a particular time period. Figure 3 summarizes the attributes of the nodes modeling a single airport. Some additional notation is needed: I = Set of nodes indicating airports. Jt i = Node indicating destination of flights departing airport i in period t. K/ = Node indicating origin of flights arriving at airport i in period t. D/ - Set of destination airports for flights departing airport i in period t. O/ = Set of origin airports for flights arriving at airport i in period t. Ckt = Cost of flight from node k to node L s~j = Flights scheduled from airport i to airport j during the time horizon.
S.A. Zenios / Network based models for air-traffic control
171
ORD STL
ATL
ORD
ORD
STL PHL
ORD
ORD
PHL
ORD
/
~.
ORD
PHL
PHI.
Fig. 2. Network model for traffic flow control
Departures from airport i
J
Destination airports of flights departing
at period t, t+l, etc. J
airport i at period r This is the set Dlt.
o
n
t
Origin au'pons of flights arriving at airport i at period r This is the set 0 t .
,l~ }0~+1 period t, t+l, etc. F i g . 3. N o d e a t t r i b u t e s o f f l o w c o n t r o l m o d e l
172
S.A. Zenios / Network based models for air-traffic control
dj~
= Flights scheduled for airport i from airport j during the time horizon. Uot,i UO t~ = airport capacities for departing and arriving flights respectively during time period t for airport i. With this notation we may now formalize the model: (OPT2) i i E rZ - , c iiJ t ¢i + E E cr,,f/¢,, JiJ t
minimize
i~l
t
i¢l
t
subject to
Sig = ~_, x~j for all origin-destination ' pairs (i, j ) , j ~ J ] ,
(7)
dji = Y'~x~i for all origin-destination t pairs(j,i), j~K/,
(8)
X x,j 1~o/
for all destined airports i in period t, and originating airport j ~ K t, f/), = ][~ xj'
(9)
for all airports i in period t,
(10)
I~D/
Xk/=Xlm
for all k = J / ,
l ~ D t k,
m~O/,
(11) i 0 <~fK,i, <~ Uot,
(12)
i 0 ~
(13)
The objective function represents the penalties imposed for possible delays. A linear relation is employed which is consistent with technology considerations (i.e. fuel consumption, salaries, maintenance and so on). Constraints (7) and (8) account for flights leaving and arriving at every airport. Constraints (9), (10) and (11) ensure that every flight reaches its destination - and arrives from the anticipated origin. These are essentially multicommodity network constraints. All planes are considered identical for the purpose of determining optimal delays; thereafter every plane is scheduled for is destination. Finally, constraints (12) and (13) impose bounds on the maximum number of flights that an airport can accomodate, either for takeoff or landing.
3. A prototype model for the Indianapolis control sector
A modeling system was built for the Indianapolis control sector - based on the multiperiod model (OPT1) of Section 2. It serves as a prototype to demonstrate the plausibility of the proposed methodology and study computational issues related to the efficient analysis of these models. One of the key issues in the success of modeling real world systems is the ability to collect the required data in a timely fashion. The comprehensive National Airspace System Plan will provide a centralized database of aircraft scheduled for, or actually flying, the high altitude jet routes. It can serve as the core database for the model. For our prototype the following data were required: airport information, flight information and fuel burn data. Table 1 provides more details. Some data, like the airport coordinates, were available through sources used in the past by FAA, while other data had to be collected for the network model. The following sources were employed: • International Official Airline Guide (IOAG) tape providing information about the airports. • Flight Progress Strip data collected by the control center providing flight information. • Fuel Burn Model (FUB) developed by the FAA providing data on fuel burn rates for different types of aircrafts at different altitudes.
Table 1 Data requirements for network models Airport information (1) Airport Identification Code (2) Geographical coordinates Flight information (1) Flight Identification Code (2) Origin airport (3) Destination airport (4) Cruise altitude upon entrance into target sector (5) C~is e altitude o n exit from targe t sector (6) Time of entrance !nt0 target sector (7) Time of exit from target sector (8) Flight hemi-code specifying legitimate cruise altitudes Fuel burn data (1) Aircraft type (2) Fuel burn rate per hour for every legitimate cruise altitude (3) Fuel burn rate per nautical mile for every legitimate cruise altitude
S.A. Zenios / Networkbasedmodelsfor air-traffic control Flight progress data were collected for a high traffic period on January 9, 1985• A total of 185 aircraft crossed the sector over a six-hour period. The duration of flights through the sector ranged from 4 to 23 minutes. Five distinct cruise altitudes above 29 000 feet were selected by the planes. The model was built based on the following considerations: (1) Allow for delays at the origin airport - up to three 10-minute intervals - and similar delays at the destination airports. This time grid can be made finer by considering a larger number of progressively smaller delay intervals (for example six five-minute intervals). The added accuracy will be balanced by the larger network that has to be solved. (2) Allow for every plane to follow one of two alternative cruise altitudes besides the one currently preferred• Choice was restricted to altitudes one level above and one level below the primary altitude• Again, this restriction can be relaxed at the expense of generating larger network models• The aircrafts could be instructed to follow any one of the five legitimate cruise altitudes specified for a particular flight. This model includes the dual objectives of assessing risk and cost as discussed in Section 2. Transportation and delay costs are linear functions of fuel burn data as obtained from the FUB model. The risk analysis is based on the occupancy coefficients. Example. If a single plane travels at altitude 35 000 feet, and the plane is in the target sector during the whole interval of 30 min its contribution to the risk function is 0. If a second aircraft travels at 35 000 feet for 30 rain the corresponding risk coefficient is 1. If the second plane stays in the target sector for 15 min the risk coefficient is 0.5. The resulting nonlinear network problem has 1295 nodes, 2873 arcs and approximating 15 000 non-zero occupancy coefficients for interacting flights. Using a finer time discretization and allowing the planes to follow five alternative altitudes would result in a model with 2405 nodes and 8510 arcs. A global model of all 500 sectors of the U.S. airspace gives a network with millions of arcs and nodes and is beyond the computing capabilities of existing algorithms. A modeling system - called F A A N E T - was built to integrate the available databases with network generators and the optimi-
173
ooT GENERATE NETWORK • ALTERNATIVE CRUISE ALTITUDES • POSSIBLE DELAYS • T R A N S P O R T A T I O N COSTS
IDENTIFY
RISK
1
• I N T E R A C T I N G ARCS • FLIGHTS AT S A M E A L T I T U D E • FLIGHTS IN T A R G E T SECTOR S I M U L T A N E O U S L Y
PREPROCESSOR
l
• W R I T E N E T W O R K IN T H E F O R M A T OF THE OPTIMIZATION CODE • WRITE RISK COEFFICIENTS
1
OPTIMIZATION
CODE
Fig. 4• The FAANET modeling system zation software. The overall modeling process is shown in Figure 4.
4. Solution methodologies 4.1. Numerical optimization The size of the models studied in this paper lies beyond what is considered computationally feasible with off-the-shelf optimization packages on general purpose computers• However, it is recognized that models (OPT1) and (OPT2) have the special structure of network optimization problems with nonlinear and nonseparable objectives: (NLGN) minimize
F(x)
subject to
A x = b,
(14)
1 ~ x ~ u,
(15)
where F ( x ) is a convex (possibly nonseparable) continuously differentiable real-value function, x is the vector of decision variable, A is a generalized node-arc constraint matrix, l and u are bounds on the decision variables and b is a vector of supplies and demands. The basis of these models is a forest of l-trees (i.e. trees with one extra arc creating a loop)• This special graph structure can
174
S.A. Zenios / Network based models for air-traffic control
be exploited in designing efficient software in two ways: (1) Sparse matrix techniques reduce the data storage requirements and therefore allow us to solve very large problems. (2) The inverse of the basis B appears in products of the form B - l y and B-l~r where y is the pivot column and ~r is the first order multiplier estimates. These products can be computed efficiently using graph data structure as described in Mulvey and Zenios (1986) or Zenios (1986) and their references. Motivated by the success of several researchers in solving linear generalized networks with hundreds of thousands of nodes and arcs [see for example Glover et al. (1978), Mulvey and Zenios (1986), Brown and McBride (1986)] we developed a nonlinear programming system for the same class of problems. The software - called GENOS implements the primal truncated Newton algorithm of Dembo and Steihaug (1983) using the ideas of active set partitioning of Murtaugh and Saunders (1978). In a manner similar to most nonlinear programming procedures each truncated Newton iteration consists of two stages: a search direction routine and a step length routine. In its simplest form - presented here for completeness the algorithm proceeds as follows: Step O. Initialization. Set k = 0 and determine a starting point x k. Step 1. Compute a search direction using a quadratic approximation to the objective function. Solve inexactly Newton's equation -
Although interested readers should consult the references, some explanations are in order. The constraint matrix is partitioned into basic (B), superbasic (S) and nonbasic (N) submatrices respectively: A = [ B I S I N ] . Decision variables are partitioned in a similar fashion x = [x s I xslxN]. The partitioning is based on first order estimates of Lagrange multipliers. In Step 1 a descent direction (Ps) is computed for the superbasic variables and the basic variables are adjusted along PB to maintain feasibility. Nonbasic variables remain fixed. Matrix Z is a projection matrix that spans the null space of A,
Z=
-B-1S) I 0
•
This matrix is computed efficiently using techniques from network flow theory to form the product B-aS. I and 0 indicate the identity and null matrices respectively and they are not computed or stored explicitly The forcing sequence parameter r k determines the accuracy in solving Newton's equations in Step 1, in a fashion that enhances the efficienc3>of our implementation while preserving the convergence properties of the algorithm. The truncated Newton algorithm and its implementation is described in detail by Ahlfeld et al. (1987) and Zenios (1986) where interested readers are referred.
4.2. Supercomputing [ Z t H ( x k ) Z ] p s = - Z t g ( x k) + r k, where H(xk), g(x k) denotes the Hessian matrix and gradient vector of the objective function respectively. Compute p B -- ( - B - ~S) Ps and set p N =0.
Step 2. Perform one-dimensional search along pk = [Pal Psi PN]: minimize
F( x t + a p t ) .
Let a* denote the solution. Step 3. Take a descent step: x t + l = X k + o t * p k.
Step 4. Check for convergence, set k = k + 1 and go to Step 1.
The development of specialized algorithms that capitalize on the problems intrinsic structure is only one approach to solving the air-traffic control models. Going a step further we examined the possible use of supercomputing in solving the developed models. The software system GENOS was installed on a CRAY X-MP/24 supercomputer at Boeing Computer Services. The CRAY X-MP is a vector computer with multiple independent functional units capable of achieving peak computational performance of 210 MFLOPS (Million Floating Point Operations per Second). In practice, performance of approximately 180 MFLOPS is not uncommon. It became evident from the first attempts to use GENOS on the CRAY that the graph structure of the underlying network prob-
S.A. Zenios / Networkbasedmodelsfor air-trafficcontrol lem inhibited the use of vector functional units. The program achieve an average performance of 18 MFLOPS. A timing analysis of the program for a typical test problem revealed that the following parts of the program were the most computationally intensive: (1) Conjugate Gradient routines (CG). (2) F u n c t i o n / G r a d i e n t / H e s s i a n evaluation (CALFGH). (3) Matrix partitioning routines. (4) Basis inversion routines. (5) Superbasics data management routines. (6) General purpose sorting routine. The C G routines implement the linear conjugate gradient algorithm for solving Newton's equation in Step 1. For large problems as such as 90% of the program is spent in this step. Attention was focused on the efficient coding of this part of the algorithm for the vector capabilities of the CRAY. The advantage of C G for large scale computing is that it does not require explicit computation of the inverse reduced Hessian ZtH(xk)Z. Instead it forms a series of products of the form Ztv, H ( x k) V and Zo, where v is a vector defined iteratively during the algorithm. The sparsity of matrices Z and H(x k) inhibited the efficient vectorization of the matrix-vector products in CG. The implementation had to be changed as described below in order to achieve the required efficiency:
175
(1) Matrices Z and H(x k) were accessed by column and row respectively and not on an element by element basis as in the original implementation. Therefore the required matrix-vector products are formed as a series of inner vector products, and not as a series of scalar multiplications and additions. (2) The vector products of item (1) are sparse and would not execute in CRAY vector mode. We first have to gather the elements of the vector v corresponding to the non-zero components of Z and H(xk). The vector products are now dense; they are computed in vector mode by the add and multiply functional units of the CRAY and the results are scattered according to the sparsity structure of the matrices. (3) Sufficient care is also taken to arrange data in memory so that the required quantities would be available during every clock cycle to avoid delaying the vector function units in item (2) above. The performance of the modified GENOS code is illustrated in Table 2. In particular the last two columns of the table indicate that the vectorized program is about three times faster than the original version. Additional work was required to structure the program in taking advantage of the multiple independent CPUs of the CRAY. An additional improvement by a factor of 1.8 was observed on a CRAY X - M P / 4 8 with two CPUs dedicated to the program. The restructured program would achieve performance in excess of 160
Table 2 Performanceof the network optimization softwareGENOS Problem (nodel/arcs)
NLPNETG Solution time in CPU seconds VAX 11/750 IBM3081-D
CRAY X-MP Original
CRAY X-MP Vectorized
PTN150 (150/196)
23.86
1.98
0.328
0.009
PTN660 (666/906)
297.93
22.85
2.435
1.227
SMBANK (64/117)
21.50
2.64
0.358
0.129
BIGBANK (1116/2230)
9100.00
376.74
244.107
46.194
GROUPlac (200/500)
1652.00
218.64
18.668
4.983
176
S.A. Zenios / Network based models for air-traffic control
MFLOPS. The modified program on the CRAY achieved a thirty-fold improvement in performance over the IBM 3081 version and more than 300-fold improvement over the VAX 11/750. Further details on the use of vector computers in solving large scale network problems are given in Zenios and Mulvey (1986).
5. Computational results A series of computational experiments were conducted in order to investigate the efficiency of the developed software, demonstrate the feasibility of the proposed methodology and study the developed prototype models. The Indianapolis prototype - and a companion miniature test problem - were first solved with the general purpose optimizer MINOS of Murtaugh and Saunders (1977) on a VAX 11/750 running UNIX 4.3BSD. The performance of MINOS in solving these problems is summarized below: Test problem (28 nodes, 65 arcs): 15.22 CPU seconds, Indianapolis prototype (1295 nodes, 2873 arcs):
ficiency with which these models are solved. To study the ability of the models to give acceptable responses in managing the traffic over high altitude jet routes we developed a histogram (Figure 5) of the number of flights at 35000ft under current system conditions, and following the use of the model. We note that 35000ft is currently the most congested altitude. The model indicates how flights should be diverted to 31000ft and 39000ft in order to achieve uniform congestion over all altitudes. The diversion of flights is determined taking into consideration the cost incurred by any decision in terms of fuel consumption in the context of the multi-criteria objective function employed. For example, the solution of Figure 5 was obtained by assigning equal weights to the cost and congestion components of the objective (i.e. w1 = w2 = 0.5). To study potential tradeoffs between fuel efficient operation of the system and the even distribution of traffic on different altitudes we developed the efficient frontier of Figure 6. This is obtained by varying the weights wl, wE such that their sum is equal to one. Point A of the figure indicates the current state of the system. The
-- 31 CPU hours.
Indianapolis sector 39000 fx
The same problems were solved with the network optimizer under identical computing conditions as follows: Test problem:
1.24 CPU seconds,
Indianapolis prototype:
--- 15 CPU minutes.
Both algorithms solve the problem to the same level of accuracy and the objective function values are identical to six decimal points. G E N O S was then used to solve the problems on a CRAY X - M P / 2 4 at Boeing Computer Services. The results are summarized below: Test Problem: 0.17 CPU second, Indianapolis prototype: 3.50 CPU seconds.
,~ D
i' TT ,T 0
8
TT 12
TTTT; 16
20
. ~ Before OptJm,zauon t ~ A f l e r Optimization
12
9~i 6~3[0 0
I TT1 TIT T 16
3 1 0 0 0 ft
This set of runs demonstrates the feasibility of the proposed methodology. Large scale instances of air-traffic control models can be solved within a few seconds with network optimization algorithms on a supercomputer. This makes possible the use of the developed models as operational planning tools under real time conditions. The quality of the results obtained from the model is certainly more important that the el-
Before Optimization After Optimization
[TTTT 4
Before Optimization After Optimization
8
T TTTTTT 12
16
20
Fig. 5. Histogram of flight congestion at 31, 35 and 39 thousand feet
S.A. Zenios / Network based models for air-traffic control
177
RISK 20
A
c;
0 19
I
I
I
I
20
21
22
23
I~
24
(x 10:) C O S T (lb)
Fig. 6. Cost-congestion efficient frontier for the Indianapolis sector
efficient frontier is not meant to serve as a tool for direct tradeoffs between congestion and cost. Such choices are often difficult and highly controversial and should be made through an informed political process. The model provides essential input to such a process. It is also useful in comparing the current state of the system with possible modifications or extensions.
6. Conclusions and future research
This paper discussed an optimization planning model for the U.S. air-traffic system. Use of this model was demonstrated with real data from a representative sector. By developing specialized algorithms we were able to solve the resulting problems efficiently, thus demonstrating the feasibility of the approach for strategic planning. In addition, taking advantage of recent technological developments in supercomputer design, we were able to solve very large problems in a matter of seconds. This research has established that network models can be used as the basis for assessing some aspects of the U.S. airspace. The technology - in terms of computer systems and algorithms is available and the required data can be collected. Further research is needed, however, in the mathematical representation of the air-traffic system and in integrating the global solution of model (OPT2) with the more myopic recommendations of (OPT1). This study considers a deterministic model of risk. In general, a stochastic program
which considers systemwide control aspects would be more appropriate. While congestion was used as a risk surrogate, other criteria - like expected number of aircraft conflicts, or expected number of controller interventions - should be evaluated as alternatives. This is an area of potential interface between the optimization model described here and the probabilistic models of Odoni and Endoh (1983). While substantial progress has been made towards building and verifying the model, additional work needs to be done in model validation. Alternative strategies generated by the model have to be compared with existing methodologies to establish the correspondence of the model to the perceived reality. This is another area of potential interface between the optimization model and the general probabilistic framework presented in Powell, Mulvey and Babu (1984). Additional work is needed in developing efficient solution algorithms for model (OPT2) that capitalize on the multicommodity structure of the problem and the dynamic nature of the air-traffic system. The two models considered in this paper, namely (OPT1) and (OPT2), should finally be integrated in a modeling system in order to complete a practical tool.
Acknowledgement
The work described in this report was partially supported by NSF grant DCR84-01098, NASA
178
S.A. Zenias / Network based models for air-traffic control
g r a n t N A G - I - 5 2 0 , a n d a g r a n t f r o m the R e s e a r c h F o u n d a t i o n o f the U n i v e r s i t y o f P e n n s y l v a n i a .
References Our Burdened Skies, Special Issue of 1EEE Spectrum, November, 1986. Safety Issues in Air Traffic Systems Planning Design (1983), Princeton University, Princeton, NJ. Ahlfeld, D.P., Dembo, R., Mulvey, J.M., and Zenios, S.A. (1987), "Nonlinear programming on generalized networks", ACM Transactions on Mathematical Software 13 (4), 350367. BieUi, M., Calicchio, G., Nicoletti, B., and Ricciardelli, S. (1982), "The air-traffic flow control problem as an application of network theory", Computers and Operations Research 9 (4), 265-278. Booth, G.R., and Harvey, C.W. (1983), "Minimizing air-traffic delays by mathematical programming", in: Safety Issues in Air Traffic System Planning Design, Princeton University, Princeton, NJ. Brown, G.G., and McBride, R.D. (1984), "Solving generalized networks", Management Science 30 (12). Dembo, R.S., and Steihaug, T. (1983), "Truncated-Newton algorithms for large-scale unconstrained optimization", Mathematical Programming 26, 190-212. Ferguson, A.R., and Dantzig, G.B. (1957), "The allocation of aircraft to routes", Management Science, 3, 45-73.
Friedman, M.F. (1986), "Optimal timing of speed control in air-traffic control", Paper presented at the TIMS/ORSA National Meeting, Miami, FL. Glover, F., Hultz, J., Klingman, D., and Stutz, J. (1978), "Generalized networks: A fundamental computer-based planning tool", Management Science 24 (12), 1209-1220. McDonald, B. (1985), Private communication, 1985. Mulvey, J.M., and Zenios, S. (1985), "Solving large scale generalized networks", Journal of Information and Optimization Science 6, 94-11. Murtaugh, B., and Saunders, M. (1978), "Large-scale linearly constrained optimization", Mathematical Programming 14, 41-72. Odoni, A., and Endoh, S. (1983), "A general model for predicting the frequency of air conflicts", in: Safety Issues in Air Traffic Systems Planning Design, Princeton University, Princeton, NJ. Powell, W.B., Mulvey, J.M., and Babu, S. (1984), "Development and validation of models of air-traffic control performance", Final report, DOT DTFA 03-82-c-400046, Princeton University. Zenios, S.A. (1986), "Sequential and Parallel Algorithms for Convex Generalized Network Problems and Related Applications", Doctoral dissertation, Civil Engineering Department, Princeton University, Princeton, NJ. Zenios, S.A., and Mulvey, J.M. (1986), "Nonlinear network programming on vector supercomputers: A study on the CRAY X-MP", Operations Research 34 (5) 667-682.
European Journal of OperationalResearch50 (1991) 166-178 North-Holland
Theory and Methodology
Network based models for air-traffic control Stavros A. Zenios Department of Decision Sciences, The Wharton School, University of Pennsylvania, Philadelphia, PA 19104, USA
Abstract: Several aspects of air-traffic control fit into the framework of network optimization models with nonlinear and nonseparable objectives. We develop two network based models for air-traffic control over the U.S. airspace. The first model deals with the dual objectives of monitoring congestion over high altitude jet routes while minimizing transportation cost. The second deals with the optimal flow of departing traffic through major airports to avoid delays at the destination. Of particular interest is the first model, due to the high increase of traffic volume over the last few years and the resulting congestion and increased number of accidents. A prototype model is built using data from a control sector in the Midwest. By integrating recent advances in numerical optimization, network flow theory and supercomputing we were able to solve the resulting problem very efficiently, thus demonstrating its feasibility. Computational experiments with real data on a CRAY X-MP/24 supercomputer are discussed.
Keywords: Traffic management, network optimization, supercomputing
1. Introduction
On an average day in 1984 approximately 17 000 aircraft traveled the high altitude jet routes over the continental United States. Takeoff and landing of one such plane are only two operations among millions that the Air Traffic Control (ATC) system handles each year. Some 14 000 controllers monitor the flow of traffic on these routes working in three basic facilities: 398 airport towers, 18 TRACON (Terminal Radar Approach Control) and 20 en route control centers. The complexities of the ATC are analyzed in a series of articles in a special issue of IEEE Spectrum (1986).
ReceivedNovember1987; revised May 1989
Air-traffic controllers strive to prevent collisions and to keep the traffic moving efficiently. They are responsible for keeping each plane in their jurisdiction separated by at least 1 000 feet vertically and by at least 5 miles horizontally. Remarkably, despite the potentials for failure, most airline flights are uneventful. However, traffic volume has been increasing steadily and with it have increased the risk of accidents and potential inefficiences. The number of reported near-collisions in mid-air increased from 584 in 1984 to 758 in 1985. It is expected that traffic volume will increase by 62% over the next decade. In response, the Federal Aviation Administration (FAA) plans to spend $11 billion on a program to update the air-traffic system - called the National Airspace System Plan. One of the primary objectives of this plan is to develop suitable procedures for assessing, and ultimately controlling, the efficiency and
0377-2217/91/$03.50 © 1991 - ElsevierSciencePublishersB.V.(North-Holland)
S.A. Zenios / Network based modelsfor air-traffic control
safety of the system. In this paper we propose two network based models for the U.S. air-traffic control. The first model deals with the dual objectives of monitoring congestion over the high altitude jet routes while minimizing transportation costs. The second model deals with the optimal flow of departing flights through major airports to avoid delays at the destination. The use of mathematical programming tools for studying problems in civil aviation originated with Ferguson and Dantzig (1957). More recently Bielli et al. (1982) developed a linear multicommodity network model for traffic flow near the Rome airport. Booth and Harvey (1983) proposed an integer programming formulation for minimizing traffic delays. Friedman (1986) discussed a pacing algorithm for determining the optimal velocity of aircrafts for meeting a prescribed schedule. The aforementioned models brought insight to the flight control problem and were well received by practitioners - refer to the proceedings of the Princeton Conference (1983). However, although the plausibility of these models is evident, they are far from being implementable. First, they do not consider in an integrated fashion the risk associated with a particular economic-efficient policy. Second, no attempt was made to evaluate these models using real data from large scale instances or to develop computationally tractable algorithms for their analysis. The work described here deviates from earlier research in several factors: (1) it integrates in a unified framework the congestion (as a risk surrogate) and cost components of air-traffic control; (2) emphasizes computationally efficient solution of these models by integrating advances in numerical optimization, network flow theory and supercomputing; and (3) evaluates the applicability of these models using real data from a representative sector of the Indianapolis control center over a six-hour time interval. The paper is organized as follows: Section 2 gives details of the air-traffic control models describing the network structure. The development of a prototype for the Indianapolis control sector is discussed in Section 3, followed by a presentation of the solution methodology (i.e. numerical optimization and supercomputing) in Section 4. Computational results are presented in Section 5.
167
The final section discusses possible uses of the model, identifies directions for further research and summarizes the conclusions.
2. Air-traffic control modeling Airspace, for the purpose of air-traffic control, is divided in a number of ways. At heights above 18000 feet all airspace is considered part of the ATC system. This controlled airspace is subdivided horizontally and vertically into some 500 sectors, each one handled by an air-traffic controller and an assistant. As a plane passes from one sector to another it is handled off to a different controller. At the lower altitudes the airspace over busy terminals is designated a terminal control area. All aircraft entering this space - shaped like an inverted cone - must be in contact with a controller. Below the sector airspace and around the terminal control areas general aviation aircraft operates under visual flight rules. In addition to monitoring traffic through individual sectors by the controllers, a computer network calculates the nationwide flow of traffic into some 40 key airports. If an airport is becoming clogged and more airplanes can arrive than can be accomodated, a schedule of ground delays is computed to keep planes from taking off for the congested airports. The overall flow of traffic is therefore determined by possible delays introduced by the centralized system or by alternative routing decisions made by the controllers. We describe now mathematical modeling approaches to these two components of the system. The first model looks at the problem from the controller's point of view, determining optimal altitudes and possible delays to avoid congestion. The second model studies the global problem of smooth traffic flow, prescribing optimal delays at key airports. 2.1. Integrated cost / congestion planning model
Figure 1 depicts a multi-period model for a single sector over all possible altitudes. Nodes indicate the origin or destination of flights whose paths travel through the target sector during time intervals of interest (i.e. the time horizon). Arcs indicate interconnecting paths between origin and destination through the target sector. We do not require the origin or destination airports of a
168
S.A. Zenios / Network based models for air-traffic control
Fig. 1. Integrated cost/congestion network model
flight to be in the jurisdiction of the target sector; the flight paths however should cross the sector for a flight to be included in this model. Note that an aircraft can follow a number of alternative altitudes in traveling through a sector on a predetermined path. As we extend the model for multiple time periods we introduce additional nodes in the obvious manner: one node for each airport during every time period. Flow on arcs interconnecting the multiple time periods represent delays - either at the origin or destination airports. A collector node for every destination airport is introduced to monitor terminal facilities. The flow on the arcs entering a collector node during a particular time period indicates total arrivals and is limited due to airway capacities. The use of multipliers on these arcs (generalized network model) changes units from planes to passengers and can be used to study the utilization of airport facilities. The underlying optimization model includes the dual objectives of optimizing congestion and cost. Cost is measured in terms of fuel consumption, either en route or due to delays at takeoff or landing. Other costs, like crew salaries and a surrogate for customer dissatisfaction, can also be included. Quantifying the risk component of the system is a critical and highly controversial issue. For example, Odoni and Endoh (1983) consider a probabilistic measure of risk at the micro-event
level and McDonald (1985) measures risk as a function of controllers' workload. For the purpose of our analysis we are optimizing a relative risk measure: the key idea is to compare some risk norm under different scenaria or system states. It is recognized that this measure of risk does not translate into a probabilistic notion of collisions. However, it provides a measure of a preferred quantity (less congestion is preferable than more congestion). As such it can assist in identifying an efficient state of a control sector with low congestion. Congestion, as specified through appropriate occupancy coefficients, is employed as a risk surrogate. The coefficients are defined as follows: Lt Ot(i,j),(k,I)
= (Time spent by aircrafts on route (i, j ) on same altitude L as aircrafts on route (k, 1) during period t) / ( T i m e spent by aircrafts on route (i, j ) on altitude L during period t). Occupancy coefficients indicate the fraction of time that two planes are in the vicinity of each other thus contributing to the congestion of the sector. See Section 3 for a sample calculation of the coefficients.
S.A. Zenios / Network based models for air-traffic control
The system-wide relative risk is then defined as
R(x) = ~, ~., ~ t
L,
1
Lt Lt~ 2 gt(i,jXk,l)~XijXkl ) '
L (i,j)(k,l)
where x ,L, and x[5' indicate the number of flights on paths (i, j ) and (k, l) at altitude L and time period t. We introduce some additional notation in order to formalize the model: i, j,
= Origin airport node at time period t. = Destination airport node at time period
J
= Collector node (sink) for destination airport j, for all t. = N u m b e r of aircrafts staying at origin node i during time period t to t + 1. This is a convenient notation for x66+. = N u m b e r of aircrafts delayed at destination node j during period t to t + 1. This is a convenient notation for xj,j,+. = N u m b e r of aircrafts from airport i to airport j during period t at flight altitude L. This is a convenient notation for x ttj (~(t" = N u m b e r of flights at airport j granted permission to land at time period t. This is a convenient notation for x j, j. = Fuel consumption of flights delayed on the ground at origin airport i. = Fuel consumption of flights delayed airborne at destination airport j. = Fuel consumption of flights from i to j traveling at altitude L. = N u m b e r of flights originating at airport i during period t. = Bound on the number of flights that airport i can accomodate in takeoff queue. = Bound on the number of flights that airport j can accomodate in landing queue. = Bound on the number of passengers that airport j can accomodate. = Average number of passengers per plane landing at airport j during period t.
t.
x~
x~
x~ '
x~j
CG, CAj L C&, s6 UT,
ULj
dJ
m',.,
The subscript t of a parameter or variable is dropped if there is no ambiguity in doing so. For example we write CG~ instead of CGi, since this
169
parameter is time-independent; likewise for CA j, U T i, U L j and dj. The time indices t are used to indicate a time slice of fixed duration - in our empirical investigation it is set equal to approximately ½ hour. During this interval we consider all the flights through the sector both for congestion and fuel consumption via the variables xij. Lt Of course this assumes that all the planes through the target sector will take the same amount of time to travel through. This is a simplification that assumes a circular shape for the sector and identical speeds for all the aircraft. While the first assumption is not unreasonable, the second one is in general not true. Hence we use as a time interval the m a x i m u m time tmax it takes for an aircraft in the model to go through the sector. Any aircraft that travels through the sector is less than tmaX is still assumed to be in the sector for the duration tmax. Hence the risk measure is an upper bound on the true congestion. A limitation with this choice of time interval is that airplanes are assumed to arrive at their destination at the same time (since they spent identical times to travel through the sector). As a result the impact of the schedule on the destination is not known and the model is of little value there. The purpose of the model is primarily to optimize flow of traffic and congestion at the target control sector. The impact of this model on congestion at the origin and destination airport can be limited by allowing only a small number of short delays in our prototype up to three 10-minute intervals. Hence the developed schedule will be close to a prototype macro-schedule developed to optimize congestion at airports (see the model of Section 2.2). The more flexibility is given to (OPT1) the less congested this sector will be at the expense of long delays and congestion at the origin and destination airports. With this discussion on notation and congestion function we may now formalize the model: (OPT1) minimize W1
{ ~ ~ CGix[ + ~2 ~ " t
it
t
CAx;- + Y'. ~
j,
+ w2R(x) subject to t--1
E E x L t "-~ X~ = Si-]- Xi L j,
L
~
t (i,j)
Ci9x[)t}
S.A. Zenios / Network based raodels for air-traffic control
170
for all origin airports i and periods t,
(1)
~f~xLt ~,.., ,j +xjt - 1 = x jtj + x j t L
it
for all destination airports j and periods t, (2) m jtj x j jt <~dj
for all destination airports j and periods t,
(3)
0 ~
(4)
0 ~
(5)
for all destination airports j,
W 1 "{" W2 = 1 .
(6)
The first two terms of the objective represent fuel consumption either due to delays on the ground before take off or delayed landing at the destination airport. The third term accounts for transportation costs for flights over all possible altitudes. The last term in the objective reflects the risk measure as discussed earlier. Constraints (1) are flow continuity constraints at the origin airports. All flights scheduled for departure from an airport are allowed to do so with some (possible) delay as indicated by the variable xj. Constraints (2) impose similar conditions at the destination airports. Constraints (3) account for all passengers destined for airport J and it can handle limits on airport facilities. If all multipliers are set equal to one (m~s = 1), then the model accounts for aircrafts and not passengers. This set of constraints is useful in monitoring passenger service facilities (e.g. limits imposed by waiting room capacity, baggage claim areas, customs inspection queues and so on). Constraints (4) and (5) impose a limit on the number of aircrafts a particular airport can handle in the takeoff or landing queues respectively. The parameters wt and w2 are goal programruing weights used in the context of the multicriteria optimization model developed. They are not meant to be set a priori. Instead they are used as part of a sensitivity analysis to identify an efficient frontier for the operation of the system. This issue is discussed further in a later section.
2.2. Air-traffic flow control Model (OPT1) considers flights with paths through a control sector of interest. It can be used as a real time operational model to determine the number of flights allocated to the different al-
titudes of the target sector given a schedule of flights over a short time period (for example a few hours). It also estimates (short) delays for different flights to optimize the operation of the target sector. We study now a global model that considers the flow of traffic through major airports at a nationwide level. The objective of this model is to minimize delays of departing flight while avoiding congestion at the destination airports. This target is a achieved by imposing bounds on the throughput capacity of both origin and destination airports. The optimization is achieved on a system-wide basis without any consideration of congestion at intermediate sectors. This problem is currently solved by simulation as discussed in the special issue of Spectrum (1986). As a result of this model control sectors may become highly congested when flights from different origin/destination pairs are going through. Model (OPT1) would be used to optimize the control sectors, creating of course more congestion at the airports than the congestion dictated by the model we introduce next. Integrating the two models successfully still remains an open question. A simple network model studying four airports is depicted in Figure 2. Every airport (represented by triangular nodes) has a set of inbound and a set of outgoing arcs indicating arriving and departing flights. Every outgoing arc corresponds to a particular time period - flow on these arcs indicates departing flights during every time interval. Incoming arcs represent the arrival time of flights. A node is used to represent departing flights for every time period and a second node represents arriving flights. A third set of nodes indicate the destination (or origin) airports for all flights during a particular time period. Figure 3 summarizes the attributes of the nodes modeling a single airport. Some additional notation is needed: I = Set of nodes indicating airports. Jt i = Node indicating destination of flights departing airport i in period t. K/ = Node indicating origin of flights arriving at airport i in period t. D/ - Set of destination airports for flights departing airport i in period t. O/ = Set of origin airports for flights arriving at airport i in period t. Ckt = Cost of flight from node k to node L s~j = Flights scheduled from airport i to airport j during the time horizon.
S.A. Zenios / Network based models for air-traffic control
171
ORD STL
ATL
ORD
ORD
STL PHL
ORD
ORD
PHL
ORD
/
~.
ORD
PHL
PHI.
Fig. 2. Network model for traffic flow control
Departures from airport i
J
Destination airports of flights departing
at period t, t+l, etc. J
airport i at period r This is the set Dlt.
o
n
t
Origin au'pons of flights arriving at airport i at period r This is the set 0 t .
,l~ }0~+1 period t, t+l, etc. F i g . 3. N o d e a t t r i b u t e s o f f l o w c o n t r o l m o d e l
172
S.A. Zenios / Network based models for air-traffic control
dj~
= Flights scheduled for airport i from airport j during the time horizon. Uot,i UO t~ = airport capacities for departing and arriving flights respectively during time period t for airport i. With this notation we may now formalize the model: (OPT2) i i E rZ - , c iiJ t ¢i + E E cr,,f/¢,, JiJ t
minimize
i~l
t
i¢l
t
subject to
Sig = ~_, x~j for all origin-destination ' pairs (i, j ) , j ~ J ] ,
(7)
dji = Y'~x~i for all origin-destination t pairs(j,i), j~K/,
(8)
X x,j 1~o/
for all destined airports i in period t, and originating airport j ~ K t, f/), = ][~ xj'
(9)
for all airports i in period t,
(10)
I~D/
Xk/=Xlm
for all k = J / ,
l ~ D t k,
m~O/,
(11) i 0 <~fK,i, <~ Uot,
(12)
i 0 ~
(13)
The objective function represents the penalties imposed for possible delays. A linear relation is employed which is consistent with technology considerations (i.e. fuel consumption, salaries, maintenance and so on). Constraints (7) and (8) account for flights leaving and arriving at every airport. Constraints (9), (10) and (11) ensure that every flight reaches its destination - and arrives from the anticipated origin. These are essentially multicommodity network constraints. All planes are considered identical for the purpose of determining optimal delays; thereafter every plane is scheduled for is destination. Finally, constraints (12) and (13) impose bounds on the maximum number of flights that an airport can accomodate, either for takeoff or landing.
3. A prototype model for the Indianapolis control sector
A modeling system was built for the Indianapolis control sector - based on the multiperiod model (OPT1) of Section 2. It serves as a prototype to demonstrate the plausibility of the proposed methodology and study computational issues related to the efficient analysis of these models. One of the key issues in the success of modeling real world systems is the ability to collect the required data in a timely fashion. The comprehensive National Airspace System Plan will provide a centralized database of aircraft scheduled for, or actually flying, the high altitude jet routes. It can serve as the core database for the model. For our prototype the following data were required: airport information, flight information and fuel burn data. Table 1 provides more details. Some data, like the airport coordinates, were available through sources used in the past by FAA, while other data had to be collected for the network model. The following sources were employed: • International Official Airline Guide (IOAG) tape providing information about the airports. • Flight Progress Strip data collected by the control center providing flight information. • Fuel Burn Model (FUB) developed by the FAA providing data on fuel burn rates for different types of aircrafts at different altitudes.
Table 1 Data requirements for network models Airport information (1) Airport Identification Code (2) Geographical coordinates Flight information (1) Flight Identification Code (2) Origin airport (3) Destination airport (4) Cruise altitude upon entrance into target sector (5) C~is e altitude o n exit from targe t sector (6) Time of entrance !nt0 target sector (7) Time of exit from target sector (8) Flight hemi-code specifying legitimate cruise altitudes Fuel burn data (1) Aircraft type (2) Fuel burn rate per hour for every legitimate cruise altitude (3) Fuel burn rate per nautical mile for every legitimate cruise altitude
S.A. Zenios / Networkbasedmodelsfor air-traffic control Flight progress data were collected for a high traffic period on January 9, 1985• A total of 185 aircraft crossed the sector over a six-hour period. The duration of flights through the sector ranged from 4 to 23 minutes. Five distinct cruise altitudes above 29 000 feet were selected by the planes. The model was built based on the following considerations: (1) Allow for delays at the origin airport - up to three 10-minute intervals - and similar delays at the destination airports. This time grid can be made finer by considering a larger number of progressively smaller delay intervals (for example six five-minute intervals). The added accuracy will be balanced by the larger network that has to be solved. (2) Allow for every plane to follow one of two alternative cruise altitudes besides the one currently preferred• Choice was restricted to altitudes one level above and one level below the primary altitude• Again, this restriction can be relaxed at the expense of generating larger network models• The aircrafts could be instructed to follow any one of the five legitimate cruise altitudes specified for a particular flight. This model includes the dual objectives of assessing risk and cost as discussed in Section 2. Transportation and delay costs are linear functions of fuel burn data as obtained from the FUB model. The risk analysis is based on the occupancy coefficients. Example. If a single plane travels at altitude 35 000 feet, and the plane is in the target sector during the whole interval of 30 min its contribution to the risk function is 0. If a second aircraft travels at 35 000 feet for 30 rain the corresponding risk coefficient is 1. If the second plane stays in the target sector for 15 min the risk coefficient is 0.5. The resulting nonlinear network problem has 1295 nodes, 2873 arcs and approximating 15 000 non-zero occupancy coefficients for interacting flights. Using a finer time discretization and allowing the planes to follow five alternative altitudes would result in a model with 2405 nodes and 8510 arcs. A global model of all 500 sectors of the U.S. airspace gives a network with millions of arcs and nodes and is beyond the computing capabilities of existing algorithms. A modeling system - called F A A N E T - was built to integrate the available databases with network generators and the optimi-
173
ooT GENERATE NETWORK • ALTERNATIVE CRUISE ALTITUDES • POSSIBLE DELAYS • T R A N S P O R T A T I O N COSTS
IDENTIFY
RISK
1
• I N T E R A C T I N G ARCS • FLIGHTS AT S A M E A L T I T U D E • FLIGHTS IN T A R G E T SECTOR S I M U L T A N E O U S L Y
PREPROCESSOR
l
• W R I T E N E T W O R K IN T H E F O R M A T OF THE OPTIMIZATION CODE • WRITE RISK COEFFICIENTS
1
OPTIMIZATION
CODE
Fig. 4• The FAANET modeling system zation software. The overall modeling process is shown in Figure 4.
4. Solution methodologies 4.1. Numerical optimization The size of the models studied in this paper lies beyond what is considered computationally feasible with off-the-shelf optimization packages on general purpose computers• However, it is recognized that models (OPT1) and (OPT2) have the special structure of network optimization problems with nonlinear and nonseparable objectives: (NLGN) minimize
F(x)
subject to
A x = b,
(14)
1 ~ x ~ u,
(15)
where F ( x ) is a convex (possibly nonseparable) continuously differentiable real-value function, x is the vector of decision variable, A is a generalized node-arc constraint matrix, l and u are bounds on the decision variables and b is a vector of supplies and demands. The basis of these models is a forest of l-trees (i.e. trees with one extra arc creating a loop)• This special graph structure can
174
S.A. Zenios / Network based models for air-traffic control
be exploited in designing efficient software in two ways: (1) Sparse matrix techniques reduce the data storage requirements and therefore allow us to solve very large problems. (2) The inverse of the basis B appears in products of the form B - l y and B-l~r where y is the pivot column and ~r is the first order multiplier estimates. These products can be computed efficiently using graph data structure as described in Mulvey and Zenios (1986) or Zenios (1986) and their references. Motivated by the success of several researchers in solving linear generalized networks with hundreds of thousands of nodes and arcs [see for example Glover et al. (1978), Mulvey and Zenios (1986), Brown and McBride (1986)] we developed a nonlinear programming system for the same class of problems. The software - called GENOS implements the primal truncated Newton algorithm of Dembo and Steihaug (1983) using the ideas of active set partitioning of Murtaugh and Saunders (1978). In a manner similar to most nonlinear programming procedures each truncated Newton iteration consists of two stages: a search direction routine and a step length routine. In its simplest form - presented here for completeness the algorithm proceeds as follows: Step O. Initialization. Set k = 0 and determine a starting point x k. Step 1. Compute a search direction using a quadratic approximation to the objective function. Solve inexactly Newton's equation -
Although interested readers should consult the references, some explanations are in order. The constraint matrix is partitioned into basic (B), superbasic (S) and nonbasic (N) submatrices respectively: A = [ B I S I N ] . Decision variables are partitioned in a similar fashion x = [x s I xslxN]. The partitioning is based on first order estimates of Lagrange multipliers. In Step 1 a descent direction (Ps) is computed for the superbasic variables and the basic variables are adjusted along PB to maintain feasibility. Nonbasic variables remain fixed. Matrix Z is a projection matrix that spans the null space of A,
Z=
-B-1S) I 0
•
This matrix is computed efficiently using techniques from network flow theory to form the product B-aS. I and 0 indicate the identity and null matrices respectively and they are not computed or stored explicitly The forcing sequence parameter r k determines the accuracy in solving Newton's equations in Step 1, in a fashion that enhances the efficienc3>of our implementation while preserving the convergence properties of the algorithm. The truncated Newton algorithm and its implementation is described in detail by Ahlfeld et al. (1987) and Zenios (1986) where interested readers are referred.
4.2. Supercomputing [ Z t H ( x k ) Z ] p s = - Z t g ( x k) + r k, where H(xk), g(x k) denotes the Hessian matrix and gradient vector of the objective function respectively. Compute p B -- ( - B - ~S) Ps and set p N =0.
Step 2. Perform one-dimensional search along pk = [Pal Psi PN]: minimize
F( x t + a p t ) .
Let a* denote the solution. Step 3. Take a descent step: x t + l = X k + o t * p k.
Step 4. Check for convergence, set k = k + 1 and go to Step 1.
The development of specialized algorithms that capitalize on the problems intrinsic structure is only one approach to solving the air-traffic control models. Going a step further we examined the possible use of supercomputing in solving the developed models. The software system GENOS was installed on a CRAY X-MP/24 supercomputer at Boeing Computer Services. The CRAY X-MP is a vector computer with multiple independent functional units capable of achieving peak computational performance of 210 MFLOPS (Million Floating Point Operations per Second). In practice, performance of approximately 180 MFLOPS is not uncommon. It became evident from the first attempts to use GENOS on the CRAY that the graph structure of the underlying network prob-
S.A. Zenios / Networkbasedmodelsfor air-trafficcontrol lem inhibited the use of vector functional units. The program achieve an average performance of 18 MFLOPS. A timing analysis of the program for a typical test problem revealed that the following parts of the program were the most computationally intensive: (1) Conjugate Gradient routines (CG). (2) F u n c t i o n / G r a d i e n t / H e s s i a n evaluation (CALFGH). (3) Matrix partitioning routines. (4) Basis inversion routines. (5) Superbasics data management routines. (6) General purpose sorting routine. The C G routines implement the linear conjugate gradient algorithm for solving Newton's equation in Step 1. For large problems as such as 90% of the program is spent in this step. Attention was focused on the efficient coding of this part of the algorithm for the vector capabilities of the CRAY. The advantage of C G for large scale computing is that it does not require explicit computation of the inverse reduced Hessian ZtH(xk)Z. Instead it forms a series of products of the form Ztv, H ( x k) V and Zo, where v is a vector defined iteratively during the algorithm. The sparsity of matrices Z and H(x k) inhibited the efficient vectorization of the matrix-vector products in CG. The implementation had to be changed as described below in order to achieve the required efficiency:
175
(1) Matrices Z and H(x k) were accessed by column and row respectively and not on an element by element basis as in the original implementation. Therefore the required matrix-vector products are formed as a series of inner vector products, and not as a series of scalar multiplications and additions. (2) The vector products of item (1) are sparse and would not execute in CRAY vector mode. We first have to gather the elements of the vector v corresponding to the non-zero components of Z and H(xk). The vector products are now dense; they are computed in vector mode by the add and multiply functional units of the CRAY and the results are scattered according to the sparsity structure of the matrices. (3) Sufficient care is also taken to arrange data in memory so that the required quantities would be available during every clock cycle to avoid delaying the vector function units in item (2) above. The performance of the modified GENOS code is illustrated in Table 2. In particular the last two columns of the table indicate that the vectorized program is about three times faster than the original version. Additional work was required to structure the program in taking advantage of the multiple independent CPUs of the CRAY. An additional improvement by a factor of 1.8 was observed on a CRAY X - M P / 4 8 with two CPUs dedicated to the program. The restructured program would achieve performance in excess of 160
Table 2 Performanceof the network optimization softwareGENOS Problem (nodel/arcs)
NLPNETG Solution time in CPU seconds VAX 11/750 IBM3081-D
CRAY X-MP Original
CRAY X-MP Vectorized
PTN150 (150/196)
23.86
1.98
0.328
0.009
PTN660 (666/906)
297.93
22.85
2.435
1.227
SMBANK (64/117)
21.50
2.64
0.358
0.129
BIGBANK (1116/2230)
9100.00
376.74
244.107
46.194
GROUPlac (200/500)
1652.00
218.64
18.668
4.983
176
S.A. Zenios / Network based models for air-traffic control
MFLOPS. The modified program on the CRAY achieved a thirty-fold improvement in performance over the IBM 3081 version and more than 300-fold improvement over the VAX 11/750. Further details on the use of vector computers in solving large scale network problems are given in Zenios and Mulvey (1986).
5. Computational results A series of computational experiments were conducted in order to investigate the efficiency of the developed software, demonstrate the feasibility of the proposed methodology and study the developed prototype models. The Indianapolis prototype - and a companion miniature test problem - were first solved with the general purpose optimizer MINOS of Murtaugh and Saunders (1977) on a VAX 11/750 running UNIX 4.3BSD. The performance of MINOS in solving these problems is summarized below: Test problem (28 nodes, 65 arcs): 15.22 CPU seconds, Indianapolis prototype (1295 nodes, 2873 arcs):
ficiency with which these models are solved. To study the ability of the models to give acceptable responses in managing the traffic over high altitude jet routes we developed a histogram (Figure 5) of the number of flights at 35000ft under current system conditions, and following the use of the model. We note that 35000ft is currently the most congested altitude. The model indicates how flights should be diverted to 31000ft and 39000ft in order to achieve uniform congestion over all altitudes. The diversion of flights is determined taking into consideration the cost incurred by any decision in terms of fuel consumption in the context of the multi-criteria objective function employed. For example, the solution of Figure 5 was obtained by assigning equal weights to the cost and congestion components of the objective (i.e. w1 = w2 = 0.5). To study potential tradeoffs between fuel efficient operation of the system and the even distribution of traffic on different altitudes we developed the efficient frontier of Figure 6. This is obtained by varying the weights wl, wE such that their sum is equal to one. Point A of the figure indicates the current state of the system. The
-- 31 CPU hours.
Indianapolis sector 39000 fx
The same problems were solved with the network optimizer under identical computing conditions as follows: Test problem:
1.24 CPU seconds,
Indianapolis prototype:
--- 15 CPU minutes.
Both algorithms solve the problem to the same level of accuracy and the objective function values are identical to six decimal points. G E N O S was then used to solve the problems on a CRAY X - M P / 2 4 at Boeing Computer Services. The results are summarized below: Test Problem: 0.17 CPU second, Indianapolis prototype: 3.50 CPU seconds.
,~ D
i' TT ,T 0
8
TT 12
TTTT; 16
20
. ~ Before OptJm,zauon t ~ A f l e r Optimization
12
9~i 6~3[0 0
I TT1 TIT T 16
3 1 0 0 0 ft
This set of runs demonstrates the feasibility of the proposed methodology. Large scale instances of air-traffic control models can be solved within a few seconds with network optimization algorithms on a supercomputer. This makes possible the use of the developed models as operational planning tools under real time conditions. The quality of the results obtained from the model is certainly more important that the el-
Before Optimization After Optimization
[TTTT 4
Before Optimization After Optimization
8
T TTTTTT 12
16
20
Fig. 5. Histogram of flight congestion at 31, 35 and 39 thousand feet
S.A. Zenios / Network based models for air-traffic control
177
RISK 20
A
c;
0 19
I
I
I
I
20
21
22
23
I~
24
(x 10:) C O S T (lb)
Fig. 6. Cost-congestion efficient frontier for the Indianapolis sector
efficient frontier is not meant to serve as a tool for direct tradeoffs between congestion and cost. Such choices are often difficult and highly controversial and should be made through an informed political process. The model provides essential input to such a process. It is also useful in comparing the current state of the system with possible modifications or extensions.
6. Conclusions and future research
This paper discussed an optimization planning model for the U.S. air-traffic system. Use of this model was demonstrated with real data from a representative sector. By developing specialized algorithms we were able to solve the resulting problems efficiently, thus demonstrating the feasibility of the approach for strategic planning. In addition, taking advantage of recent technological developments in supercomputer design, we were able to solve very large problems in a matter of seconds. This research has established that network models can be used as the basis for assessing some aspects of the U.S. airspace. The technology - in terms of computer systems and algorithms is available and the required data can be collected. Further research is needed, however, in the mathematical representation of the air-traffic system and in integrating the global solution of model (OPT2) with the more myopic recommendations of (OPT1). This study considers a deterministic model of risk. In general, a stochastic program
which considers systemwide control aspects would be more appropriate. While congestion was used as a risk surrogate, other criteria - like expected number of aircraft conflicts, or expected number of controller interventions - should be evaluated as alternatives. This is an area of potential interface between the optimization model described here and the probabilistic models of Odoni and Endoh (1983). While substantial progress has been made towards building and verifying the model, additional work needs to be done in model validation. Alternative strategies generated by the model have to be compared with existing methodologies to establish the correspondence of the model to the perceived reality. This is another area of potential interface between the optimization model and the general probabilistic framework presented in Powell, Mulvey and Babu (1984). Additional work is needed in developing efficient solution algorithms for model (OPT2) that capitalize on the multicommodity structure of the problem and the dynamic nature of the air-traffic system. The two models considered in this paper, namely (OPT1) and (OPT2), should finally be integrated in a modeling system in order to complete a practical tool.
Acknowledgement
The work described in this report was partially supported by NSF grant DCR84-01098, NASA
178
S.A. Zenias / Network based models for air-traffic control
g r a n t N A G - I - 5 2 0 , a n d a g r a n t f r o m the R e s e a r c h F o u n d a t i o n o f the U n i v e r s i t y o f P e n n s y l v a n i a .
References Our Burdened Skies, Special Issue of 1EEE Spectrum, November, 1986. Safety Issues in Air Traffic Systems Planning Design (1983), Princeton University, Princeton, NJ. Ahlfeld, D.P., Dembo, R., Mulvey, J.M., and Zenios, S.A. (1987), "Nonlinear programming on generalized networks", ACM Transactions on Mathematical Software 13 (4), 350367. BieUi, M., Calicchio, G., Nicoletti, B., and Ricciardelli, S. (1982), "The air-traffic flow control problem as an application of network theory", Computers and Operations Research 9 (4), 265-278. Booth, G.R., and Harvey, C.W. (1983), "Minimizing air-traffic delays by mathematical programming", in: Safety Issues in Air Traffic System Planning Design, Princeton University, Princeton, NJ. Brown, G.G., and McBride, R.D. (1984), "Solving generalized networks", Management Science 30 (12). Dembo, R.S., and Steihaug, T. (1983), "Truncated-Newton algorithms for large-scale unconstrained optimization", Mathematical Programming 26, 190-212. Ferguson, A.R., and Dantzig, G.B. (1957), "The allocation of aircraft to routes", Management Science, 3, 45-73.
Friedman, M.F. (1986), "Optimal timing of speed control in air-traffic control", Paper presented at the TIMS/ORSA National Meeting, Miami, FL. Glover, F., Hultz, J., Klingman, D., and Stutz, J. (1978), "Generalized networks: A fundamental computer-based planning tool", Management Science 24 (12), 1209-1220. McDonald, B. (1985), Private communication, 1985. Mulvey, J.M., and Zenios, S. (1985), "Solving large scale generalized networks", Journal of Information and Optimization Science 6, 94-11. Murtaugh, B., and Saunders, M. (1978), "Large-scale linearly constrained optimization", Mathematical Programming 14, 41-72. Odoni, A., and Endoh, S. (1983), "A general model for predicting the frequency of air conflicts", in: Safety Issues in Air Traffic Systems Planning Design, Princeton University, Princeton, NJ. Powell, W.B., Mulvey, J.M., and Babu, S. (1984), "Development and validation of models of air-traffic control performance", Final report, DOT DTFA 03-82-c-400046, Princeton University. Zenios, S.A. (1986), "Sequential and Parallel Algorithms for Convex Generalized Network Problems and Related Applications", Doctoral dissertation, Civil Engineering Department, Princeton University, Princeton, NJ. Zenios, S.A., and Mulvey, J.M. (1986), "Nonlinear network programming on vector supercomputers: A study on the CRAY X-MP", Operations Research 34 (5) 667-682.