INTEGRATED RISK/COST PLANNING MODELS

FOR THE U.S. AIR TRAFFIC SYSTEM

John M. Mulvey

and

Stavros A. Zenios

(NASA-CR-17727U) INTEGRATED EISK/COST N86-26312

PLANNING MODELS FOE THE OS AIE TRAFFIC

SYSTEM (Eriaceton Univ., N. J.) 41 p

HC A03/MF A01 CSCL 17G {Jnclas

G3/04 H3602

Department of Civil Engineering

Princeton University

Princeton, New Jersey 08544

Engineering - Economic Series

Report EES-85-9

Revised, December, 1985

and
Abstract

-We_descr-ibe a prototype network planning model for the U.S. Air

,• c. «>t-a

Traffic control system/ The model encompasses the dual objectives of

managing collision risks and transportation costs where traffic flows

can be related to these objectives. The underlying structure is a

network graph with nonseparable convex costs; the model is solved efficiently by capitalizing on its intrinsic characteristics. Two specialized

algorithms for solving the resulting problems are described: (1) truncated

Newton, and (2) simplicial decomposition.

The feasibility of the approach is demonstrated using data collected

from a control center in the Midwest. Computational results with different

computer systems are presented — including a vector supercomputer

(CRAY-XMP). The risk/cost model will have two primary uses: (1) as a

strategic planning tool using aggregate flight information, and (2) as an

integrated operational system for forecasting congestion and monitoring

(controlling) flow throughout the United States. In the latter case,

access to a supercomputer will be required due to the model's enormous

size.

I'D A\

A

I \UVfc

(LiT t\ rw.-i

Vt.C-H-NjiQ.ufe

C50 (Y\9 '

\
-1-

1. Introduction

By the year 2000, the U.S. Federal Aviation Administration (FAA)

plans to spend over $11 billion oh a program to upgrade and consolidate the

U.S. air-traffic system — called the National Airspace System Plan [NASP].

On an average day in 1984 approximately 17,000 aircraft traveled these

routes. It is is expected that traffic volume will increase by 62% over

the next decade. The FAA has responsibility for managing the air-traffic

system, especially those elements affecting risk. Indeed, the agency has

been able to reduce overall (and relative) risks during the past, thirty

years [ 1 ]. As new hardware and software technologies have been implemented,

the FAA has provided mechanisms for monitoring factors associated with

system risks. Recently a trend to deregulate the airline industry in the

interest of greater efficiency and competitiveness has begun. -Despite this

trend the FAA must possess efficient procedures for assessing (and ultimately

controlling) system risks.

In this report we discuss the development of a prototype network

planning model for flights on high-altitude jet routes over the U.S. airspace, in conjunction with the NASP. The model encompasses the dual objective

of assessing risk-related measures and transportation costs. The underlying

mathematical model has the special structure of networks graph. Since safety

components depend upon interacting variables, the proposed model falls in

the category of a nonlinear network with nonseparable cost functions.

The model can serve as the building block for a management information

system that can assist FAA in two basic settings:

(a) As a strategic planning tool, or

(2) As an operational planning program for air-traffic

scheduling and routing.
-2-

In the former case it is essential that the optimization can be solved

efficiently and accurately. In this way extensive sensitivity analyses

can be carried out and answers to "what if" questions obtained in a

comparative fashion and without regard to numerical instability. The

use of the model as an operational planning tool depends on the ability

to "solve" the model under real-time conditions.

By capitalizing on the special structure of the underlying network

basis (a forest of 1-trees) we have developed highly efficient solution

algorithms [2, 12, 13]. Thus large-scale problems can be handled within

minutes or hours; repeated runs can be made and the use of the model for

planning is feasible. Going a step further we have streamlined our

algorithms for the architecture of vector supercomputers, in particular

the CRAY-XMP at Boeing Computer Services. By taking advantage of the

Cray's pipeline features we were able to solve the test cases in a matter

of seconds — thus demonstrating feasibility of the network model as an

operational planning tool.

The rest of this paper is organized as follows. The mathematical

models are defined in Section 2. A general description of two solution

algorithms —' (1) Truncated Newton, and (2) simplicial decomposition —

appears in Section 3. Data requirements for the model are listed in

Section A. Next, Section 5 presents the modeling of a representative

air-traffic control sector (Indianapolis); some computational

results are given in Section 6. Finally, we discuss limitations of the

current model and directions for future research.
2. Air Traffic Control Modeling

A stochastic programming model for the air-traffic control problem

was proposed in the early work by Ferguson and Dantzig [5]. More recently,

general aspects of air-traffic systems planning and design were the focus

of a conference at Princeton University [ 1 ]. We propose here a network

formulation for the air-traffic planning problem.

The generic nonlinear network model takes the form:

[NLGN] Minimize F(x)

Subject to:

X

 - -f u

 = V '• '

ii-lV-V*

 for (i'j) e E (2'2)

where

F(x) = convex function

{N} = set of nodes

{E} = set of arcs (edges)

x.. = flow over arc (i,j)

m, . = multiplier* on arc (k,i)

b. = supply/demand for node i

x = {x..|(i,j) c E}

'

+ f'"i,j) e E}

6^ - {j|(J,i) e E) -

..(^-.) upper (lower) bound on arc (i,j)

X = {x±.|xi. satisfies constraints (2.1) and (2.2)}

*Multipliers are indicated so that airport congestions and other factors

can be modeled. These features will be descussed In Section 7.
-4-

We may rewrite [NLGN] in the following form:

[NLGN] Minimize F(x)

Subject to x £ X

where

X E {x|A • x = b,£ _< x <: u)

Since the [NLGN] basis is a "forest of 1-trees" (collection of subtrees

with one extra arc per subtree creating a single loop), efficient procedures

are available for storing and updating the basis and other aspects of the

algorithm [9, 12].

Figure 1 depicts a simplified example graph of the first planning model.

This example consists of five airports and interconnecting routes, representing

an aggregate of individual flights into and out of the designated airports.

Note that each airport has three triangular nodes, indicating net traffic,

number of incoming flights and the number of outgoing flights. Network arcs

points in the obvious direction of traffic movements.

The corresponding optimization model is shown below:

[OPT1] Minimize {w • [ I I c* z*.] + w • s(z)} = F(*)

*• / • • \ r- f> fTt ^ J * J ^ tl,JJt-A..

Subject to:

„ = Y x ik

yJ±.

(i'k)£Al

-• 1 vn 1 "n ( k'J>

 such that

 J £

J KJ Kj J A,kEM

ik

 i£A

«

 kELi

X.. = ) Z

A

 iEA

ik I z-n,

 ltA

»
-5-

r V.

ORSG&NAL PAGE fS

OF POOR QUALITY

0)

T>

O

U

•H

u

•H

V4

•H

m

x

0)

14&

(r.
-6-

ykj

^

0 2 u «R

 and

yi

 = bi i£A

V

y

j'Xik'ykj'2ij -° *'J'k'

r

W1'W2 - °

where we have used the following decision variables:

x.(y.j) = number of total outbound (inbound) flights, airport

x., = number of outbound flights from airport i to airport k

IK

z.. = number of flights over route r, airport i to airport j

y, . = number of inbound flights from airport k to airport i.

Kl

Notation:

A = set of airports

(i,j) = feasible pair (linking airport i to airport j)

R . . = all feasible routes for pair (ij)

L. « outbound destinations for airport i

M. = inbound airports for airport j

Ax = {(i,j)|lei and

c. . = marginal transportation cost on route reR. .

b. = net traffic at airport i

F(z) = convex nonseparable function, measuring risk
-7-

This single period planning model aggregates all individual flights

between airport-pairs over the planning horizon. While extending the model

to a large number of multiple periods will considerably increase the model's

size, it would be, conceptually, straightforward to depict several periods.

As such, the [OPT1] model can be used for supporting airport resource

decisions, such as opening new runways and expanding the capabilities of

control towers. The primary aims are to identify bottlenecks and to predict

imbalances for the U.S. air-traffic system, given a variety of scenarios.

In this representation the optimization model takes a bi-criteria

objective function consisting of risk and systemwide cost, with respective

weights w2 and w.. It should be stressed that the weights are used only as

part of the sensitivity analyses to identify the efficient frontier and are

not meant to be set £ priori. Transportation cost is easily quantified in

terms of traveled distance and fuel burn rates of the different aircraft

models. A monetary value can also be assigned on factors like customer

dissatisfaction due to flight delays and so on. Quantifying the risk

component of the system, however, is a critical and controversial issue.

At best we may consider optimizing a relative risk measure: the key idea

is to compare some risk norm under different scenarios or system states,

or with other similar systems. Odoni and Endoh [15] consider a probabilistic

analysis of risk. In Section 5 we touch upon the issue of modeling risk as

a function of congestion in the target sector during time intervals of

interest.
-8-

The issue of identifying an acceptable point on the efficient riskcost frontier is complex. It is obviously difficult for society to

compare human risks and congestion delays against financial criteria.

See, for example, Rowe [18]. Ultimately the solution depends upon societal

tradeoffs and can only be derived through an informed political process.

The network framework provides essential input data for this process by

tracing out the efficient frontier; see Section 5.

The aggregate model [OPT1] is unable to provide adequate details for

operational planning purposes.

Therefore, a second network model has been developed and an example

is shown in figure 2. Here, decision variables monitor possible delays and

alternative altitudes for every flight departing to or arriving at the

airports of interest. The corresponding mathematical respresentation is

defined below. In addition to optimizing transportation and delay costs

this model limits the number of flights traversing regions of interest —

e.g., control sectors — as a surrogate for minimizing risk. Including

r M*^*«A i %*.. . r* V r* Pt Pt. ) ) C^ BQ-.» / i_/^- ^ * w.

[OPT2] Minimize J J \ c^ . x^ + I I f f ^ ^^ k k '

these aspects causes [OPT2] to be considerably larger than [OPT1].

,P

fcF pcPt teTf

 f

Subject to:

d° = 1 for all

d*"1 = d* + I x** for all feF,. teT<

+ w for all aCA, k£K
-9-

FLIGHTS AIRPORTS

STL

Figure 2: Example Air-traffic Network (Operational Mode!)
-10-

-

 U for a11 acA>

I

 xf 1 L for a11

 «S. keK

x

pt

£R

kf r

x{ £Rr

where we have used the following decision variables:

x^ "« 0,1 : flight f follows route p during time period t

df = 0,1, tcT, : flight f delays its departure during period t

fl

w for keK : number of planes delayed at destination

K, .

airport a during period k

f t XV — 0.««.y««» } I V

k max

y, k = 1,2,...,K : number of planes landing at destination

airport a during period k.

Notation:

aeA : set of airports

feF : set of flights

p,tP, : set of possible routes for flight f

k = 1,2,...,K : time of periods in planning horizon

UluX

t.eT, : set of possible time periods for departure of flight f

k

R : set of flights traversing control sector r at period k

Dj* : set of flights whose destination is airport a,

and arrival time is k

c^ : marginal transportation cost for flight f on route

p during period t

c, : marginal departure delay cost for flight f during

period t
-11-

c, , : marginal landing delay cost at airport a

during period k

3 a

u, (£, ) : upper (lower) limit on the total number of

flights landing at airport a during period k

k

L : limit on the number of flights traversing control

.sector r during period k

S : set of control sectors

This model controls individual flights: for every flight it specifies

possible altitudes, routes to be followed and any departure/landing delays.

The objective function incorporates the marginal cost of transportation,

and associated delay costs.* Transportation risk is assessed by imposing a

limit on the number of flights through a particular sector during every time

interval of interest. As mentioned this formulation provides more details

than [OPT1] since it controls flights on an individual basis. The added

details result in a more complex model: the model is not only larger than

[OPT1] but also it has integer variables. It is, however, equivalent to

[OPTl]in the risk minimization aspects where congestion in an identified

region provides a surrogate for risk. By varying the limit on the number of

flights traversing a control section during every time period we may achieve

the same result as by varying the weights in the multi-objective formulation

in [OPT1], Mathematically the first formulation is more tractible: it deals

with a continuous nonlinear network model, while the second model deals vith

a linear multicommodity network problem with integer variables.

*Note that the model minimizes total transportation costs, following a

utilitarian economies. If this approach adversely impacts individual carriers,

the model can be adjusted through multiple weightings or by means of side

payments between carriers (or the FAA).
-12-

The two models are expected to serve different purposes. The

strategic planning model is intended to assess aggregate data, while the

operational planning model deals with more microscopic aspects of airtraffic control.
-13-

3. Solution Algorithms

We describe here the main features of two nonlinear programming

algorithms that have been specialized to solve [NLGN]. The methods are:

(1) A second order algorithm based on truncated Newton search directions

(TN), and (2) A first order algorithm based on the simplicial decomposition of the feasible region (SD). Both algorithms possess distinct

characteristics: the truncated Newton algorithm can identify solutions

to a high degree of accuracy; the simplicial decomposition algorithm can

quickly converge to an approximate solution. The interested reader should

refer to the papers by Ahlfeld et al. [2] and Mulvey et al. [13] for more

detailed description of these algorithms and some computational results.

3.1 Truncated Newton Algorithm

In a manner similar to most nonlinear programming procedures, each

(TN) iteration consists of two stages: (1) a search direction routine, and

(2) a step length 'routine. Table 1 depicts the overall flow, in which the

r

notation refers to the model [NLGN] presented in Section 2.

The search direction must fulfill certain essential features so that

the overall algorithm will converge and so that performance efficiencies

are attained. First, the direction must both maintain feasibility and

point downhill (in a minimization context). Defining the search direction

—k —k th

as p and given a feasible point x at the k iteration, the usual Newton

—k

method for calculating p would solve the following quadratic programming

problem:
-14- ORIGSNAL

POOR

Table 1 : The Truncated Newton Algorithm

Me«if> MI «.'

ittiv. con»tmr.t«

tt ceritrttrtt

Ceaputi (t«r *- Mti(fyin«

centitieni (er f)eb*3

cenvtr«*nc*
-15-

• \

,„„, Minimize *'1,-k.. t^-.-k-.-k , ,-kst-k

[QP] _k -^(P ) G(x )p + g(x ) p

P

Subject to:

A • p

,, k if x.

3

., k if x.

where

g(xk) gradient of F(xk) at xk

G(xk) Hessian of F(xk) at xk .

By restricting our attention to a special projected matrix 2, whose columns

form a basis for the null space of A, i.e., A • Z = 0, the problem [QP] can

be solved using the following two formulae:

where

(Z G Z). p

k

S

-k

P

(3.1.1)

-k

Pc

Z = S

N

m

s

n-s-m

-k

G semi-positive approximation to the Hessian at point x

-k -k

g gradient of objective function at point x

and where the decision variables have been partitioned into three sets:

*

 e {xb|xs|xn]

A • [B|S|N]

Igblgslgn)

p

k
-16-

-k

The benefits of the Newton direction p are greatest in the neighborhood of a solution; however, it is expensive to calculate the solution of

equation (3.1.i). In response, we adjust in a dynamic fashion the degree

{ k\ of accuracy of solving [QP]. A forcing sequence T.TI } —> 0 is employed in

this regard. Accuracy is defined according to the relative residual in

equation (3.1.1), N^-l l , in ^ich r

k

 = (Z

t

 G Z) ' p* + Z*ik and

llzi k l ! .

 s

is a vector norm in R . The minor iteration (see table 1) continues only

I I kiI ,

until the required accuracy is attained. Thus, * J ^ ri defines

llz^i!

the termination criteria for the minor iterations.

When the algorithm is far from the solution the reduced gradient —

| JZ g | | — is large and little work is required to locate a direction

satisfying the acceptance criteria. Only the basic and super basic

variables are optimized. If one of these variables hits a bound, the

constraining variable is transferred into the set of nonbasics [x ]. As

| |Z g )| is reduced the acceptance criteria becomes more restrictive and

the current solution to the direction finding problem lies closer to the

.Newton direction.

At this point, the nonbasic variables must be tested for optimality.

First order estimates for the Lagrange multipliers are computed as follows:

*J' '4 ' '-1
-17-

In this environment non-basic variables that reduce the objective function

when moving away from their bounds (i.e., if V^ < 0 and x^

 = un> or if

V*-> 0 and xj = £J) along with free non-basic variables (i.e., ^j < xJ< «>

n n n n.n "

called eligible. Eligible variables are transferred to the superbasic set

[x ] in conjunction with a maximal basis [A], and the TN algorithm continues

the next niajor iteration with the new partition.

A sizable portion of the algorithm's execution time involves computing

-k • • the search direction p . While the truncated-Newton method can use any

s

iterative method fol solving equation (3.1.1), we have chosen the linear

conjugate-gradient [CG] method. Although the reduced Hessian matrix

Z GZ is typically dense, the product required by [CG], (ZtGZ)^, is easily

computable due to the sparsity of the large-scale components.

The success of the conjugate-gradient method depends upon locating a

"good" search direction in a small number of iterations. Thus, preconditioning

the reduced Hessian by the matrix P is important so as to reduce the number of

CG iterations. Whereas the usual initial element of the CG sequence is

-k -k

g , the vector P • g becomes an initial element when preconditioning,

S S

where P is a positive-definite matrix. See [2] for further details.

3.2 Simplicial Decomposition Algorithm

While the TN algorithm is capable of solving large nonlinear network

problems, we felt that a first order approximation would be better suited

for ultra-large examples — problems with more than 10,000 nodes and

100,000 arcs. The simplicial decomposition algorithm was selected by us to

meet this goal. This algorithm is best examined in the context of general
-18-

decomposition methods. .Following Geoffrion [6] we may place the algorithm,

ir. the class of "inner linearization/restriction" methods. The SD algorithm

iterates between (1) solving a linearized subproblem and (2) solving a

nonlinear master problem on a restricted space subject to non-negativity

constraints. Table 2 depicts the overall flow. The algorithm was first

presented in this form by von Hohenbalken [12,20]. Holloway [8] proposed

the same algorithm, as an extension of the Frank-Wolfe method, and recently

Lauphongpanich and Hearn [10] devised a restricted verion for the traffic

assignment problem. We have developed a version of SD specialized to

handle [NLGN], whereby the master problem is solved inexactly [13].

The following theorem due to Caratheodory provides the necessary

theoretical foundation for the algorithm.

Theorem:

Let X e R be a non-empty convex polytope. Then every

x e X lies in the relative interior of one of a finite

number of simplices whose vertices are extreme points of X.

See [19] for a proof. The main idea behind SD is simple in principle:

(1) First solve a linearized subproblem to get the extreme points of X.

We need not generate all the extreme points of the feasible region;

this would result in a problem as complex as [NLGN]. Instead we

generate extreme points as needed — in a manner reminiscent of

Dantzing's column generation method. A (K-l) dimensional simplex

is defined from K extreme points, and the search for the optimum

is now restricted on the generated simplex.
-19-

INITIALIZATION

K-l STARTING POINT XK

SUBPROBLEM

SOLVE

SET:

r»Yur*

wf c -£

WIN Y7f(X« )

AJ- b

Ls>su

YES

MASTER PROBLEM

SOLVE

SET:

MIN

*i . l'l.....H-l

SET : >I-I

o, L-J,

NO

Table 1: The SiBpliciftl Decomposition Algorithm
-20-

(2) The master problem attempts to optimize the original function on the

simplex generated by the subproblem. We have

Minimize F(w y)

w

Subject to:

K

Ji"'-1 . ' ' .

0 <_ w. <_ 1

where y = yi»y2 ..... Vv *s tne set of generated extreme points and

w = w1,w0,...,w are associated weights. We may reduce further the

i. / K .

dimension of the master problem making use of the implicit function

theorem:

K-l

Minimize F( T w.(y, - y ))

Subject to:

0£w. i = l ..... K-l

Thus we have a nonlinear problem in dimension (K-l) subject to simple nonnegativity constraints.

Reducing the master problem to a sequence of unrestricted (K-l)

dimensional problems, this problem can be solved by any unconstrained algorithm.

Refer to the work by Mulvey et al. [13] for computational experiments. The

solution to the master problem is more important when the solution to [NLGN]

lies in the current simplex. Thus we adjust in a dynamic fashion the

accuracy in solving the master problem. Again a forcing sequence
-21-

{r, } —> 0 is employed and the master problem terminates when

i • T I • k

I |D * g| i 2.

 r> ~~ wh616

 D is

 tne derived linear basis representing the

current simplex. Once this degree of accuracy is achieved we return to

the subproblem. " • .

.Simplicial decomposition provides us with a modular algorithm for

solving [NLGN]. It converges rapidly to a good approximate solution,

represented as a linear combination of a few extreme points (K « n).

If high accuracy is requires — i.e., exact representation of the optimum

solution — then we need generate n + 1 extreme points and the master

problem becomes as difficult as the original problem. In addition the

subproblem preserves any special structure that may be present in the

original problem. For [OPT] the subproblem iterations consists of solving

a linear generalized network problem. Code LPNETG [12] was employed,

modified to allow restarting from the basis of the last subproblem. The

master problem was solved using a Quasi-Newton algorithm. Reference [13]

provides further details for the implementation of SD and accompanying

results from computational experiments.
-22-

4. Data Requirements «

One of the most important aspects of modeling real vorld systems is

the ability to collect the required data in a timely fashion. In the case

of the Models [OPT], data may be classified in two categories: static and

dynamic. By static we mean data that do not change over a long period of

time (e.g., airport locations) and by dynamic we refer to data that change

with time (e.g., flights scheduled), or with technological innovation

(e.g., aircraft fuel burn rates, navigation systems and flight management).

For the model to be useful the input data must be readily available. A key

component of the comprehensive National Airspace System Plan is a centralized

data-base of .aircraft scheduled for, or actually flying the high altitude

jet routes [7]. This data-base can serve on a real-time basis for updating

the data required for [OPT].

For the prototype model developed the following data were required:

(I) Airports information

(II) Flights information

(III) Fuel burn data.

Table 3 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:

(I) International Official Airlines Guide (IOAG) tape,

providing information about the airports

(II) Flight Progress Strip data collected by the Control

Center, providing flights information

(III) Fuel Burn Model developed by the FAA providing data

about fuel burn rates for different types of aircrafts.
-23-

PAGE IS

OF POOR O'JALITY

Airports Information

1. Airport ID code

2. Geographical coordinates

* Flight Information

1. Flight ID

2. Origin airport

3. Destination airport

4. Cruise altitude on entering target sector

5! Cruise altitude on exit from the target sector

6. Time flight enters the target sector

7. Time flight exits from the target sector

8. Flight Hemi code defining 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

Table 3 : Model data requirements
-24-

5. Modeling a Sector of the Indianapolis Control Center

A network model was built for the airspace centrj]iec by a sector of

the Indianapolis center. The purpose of the model is to serve as a

prototype to illustrate the use of the optimization algorithms described

earlier as well as the feasibility of the proposed model. Data were

collected for a high traffic period on January 9, 1985 in which a total of

185 aircrafts crossed the sector over a 6-hour period. The duration of

flight through the sector ranged from 4 to 23 minutes. Five distinct cruise

altitudes above 29,000 feet (FL 290) were selected by the planes.

The model was built as a multi-period network as shown in figure 3.

The following provisions were made in the model:

(I) 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 (e.g., six

5-minute intervals). The added accuracy will be balanced by

the larger network that has to be solved.

(II) Allow for every plane to follow one or two alternative cruise

altitudes besides the one currently followed. Choice was

restricted to the cruise 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 four

or five legitimate cruise altitudes, specified for the

particular flight.
-25

-

Ii1t!.1u.
-26-

The model includes the dual objective of assessing risk and cost,

as proposed in Section 2. Delay cost was considered as a function of

fuel burn data. Other delay costs like crew salaries and a surrogate

for customer dissatisfaction could also be included. The risk analysis

was based on occupancy rates (congestion) at the same altitude during

time intervals of interest.

We define the occupancy at level L during a period T as:

Time spent by aircraft on route i on the same

LT e altitude (L) as aircraft on route j during period T

ij Time spent by aircraft on route i

altitude L during period T

The system-wide relative risk was then defined as:

LT

II I oft (x -x )'

T L i.j

 1J 1 J

where x.,x. indicate the number of flights on arcs i and j respectively.

Summation was taken over all legitimate cruise altitudes (L), and over ten

36-minute intervals that cover the six hour planning period. Higher

interactions — between more than two planes — were ignored. Figure 4

depicts the different phases of the modeling procedure.

Example;

(1) If a single plane travels at altitude 39000 feet, and the plane

is in the target section during the whole interval of T • 30 min.,

its contribution to the risk function is 0.

(2) If a second aircraft flies at 30000 feet for T « 30 min., the

corresponding risk coefficient is 1.

(3) If the second plane stays in the target sector for 15 min.,

the risk coefficient is 0.5.
-27-

PROGRAM

FAANET

(FUEL BURNf

\MODEL \

GENERATE NETWORK

• SI.TED-MIVE CRUISE ALTITUDE

•POSSIBLE DELAYS

•TRANSPORTATION CO£TS

IDEtJTIF Y RIS K

*P.CS

i

"'Ts AT S*M£ ALTITUDE

.TE ]N URSET SECTOR

PPEPP.-OCESSOR

•Wp:TE NETWORK IN THE FORMAT OP

T HE OPTIMIZATIO N COOE

•W»JT E K1SK COEFFICIENT S

1

OPTIMIZATION CODE

*NLPNETG

*NGSD

*MINOS

Figtire 4; Modeling the Indianapolis Control Sector
-28-

This formulation provides for two degrees of freedom in the decision

making process: delay the departure/arrival times in an optimal fashion

and instruct aircraft to follow alternative altitudes to reduce congestion.

The resulting network problem has 1295 nodes, 2873 arcs and approximately

15000 non-zero coefficients describing interacting flights. The test case

was solved using code NLPNETG. By varying relative weights on the

transportation and risk function in a systematic fashion the risk/cost

efficient frontier was traced (figure 5). Again, the efficient frontier

is not meant to serve as a direct way of comparing risk with cost. Instead

it guides one in evaluating alternative modes of operation of the airtraffic control system, as generated by the model, or with currently

followed procedures.

The major advantage of this methodology is that it generates a

sequence of alternatives that are efficient; i.e., both risk and cost

values cannot be improved simultaneously. This is easier to understand

if we notice the location of point A in figure 5 — this point was obtained

by solving the optimization problem inexactly. From point A we may move to

a series of alternative solutions for which the system is better off, both

with respect to transportation cost and risk.

To study the effect of airplane congestion, we developed a histogram of

all planes flying at a particular altitude, during the ten time intervals of

interest. Figure 6 summarizes tha results for three particular altitudes,

before and after the optimization model was used. Note that as expected planes

were diverted from a highly congested altitude (35000 ft.) to less congested

routes (31000 ft. and 39000 ft.). This result was obtained with relative

weights 0.5/9.5 on both risk and transportation costs.
-29-

20

16

12

20 21 22

(

tost (it;

Figure 5 : Efficient Frontier for Risk / Cost Tradeoff*
-30-

ORIGJNAL PAGE IS

O. POOR QUALITY

4 J « I ' i

UN fffiO*

Cftiiiiitiori JHTI iftP epiiuution

Cen'.ro! Cnlr

Altitude : 3B.OK ft

1 i 3 4 .

tut ftncds

«rti«ntjor. |IITi Ktr

Mi»«pelii Cffrttl Mrtr

: Juocc ft

Figure 6; Airport Congestion Before and After Optimization
-31-

6. Computational Results

The FAA Control model can grow in size to challenge the capabilities

of general purpose algorithms and conventional mainframe computer systems.

The developed Indianapolis prototype consists of 1295 nodes and 2873 arcs.

If instead we had chosen six 10-minute intervals for possible delays, and

five alternative cruise altitudes the resulting problem would consist of

2405 nodes and 8510 arcs. If we consider simultaneously ten control sectors,

with the same number of flights as the Indianapolis center, the network

grows to approximately 20000 nodes and 80000 arcs.

To demonstrate the efficiency of the algorithms described in Section 3,

we have solved several problems arising from a wide range of applications.

Table A presents relevant information concerning the test problems. First,

eighteen test problems were solved with the general purpose code MINOS 114].

The results from this program formed a benchmark against which to compare

the truncated Newton — code NLPNETG — and the simplicial decomposition —

code NGSD — algorithms. Problems that cannot be solved efficiently with

the general purpose code can be solved with minimal computational resources

using the specialized network algorithms. Tables 5 and 6 summarize the

results. We observe that efficiency of the specialized network codes

increases with problem size.

Advances in parallel processing computers are expected to ensure the

feasibility of new applications. Specifically the air-traffic control

model will benefit from the use of supercomputers in two domains:
-32-

(i) Efficiently solve larger models that cover more than one

control center and ultimately extent to cover the whole

continental U.S. airspace, with finer time discretization

(ii) Solve the operational planning model, under real-time

conditions.

To illustrate the situation we have specialized NLPNETG for the CRAY-XMP

vector computer at Boeing Computer Services. The optimizing compiler

proved to be only marginally effective, due to the sparsity of the network problems. We had to analyze the algorithm in a way to take advantage

of the pipeline features, and the presence of multiple vector functional

units. Table 7 summarizes some of the results. We observe that much

can be gained by specializing the network algorithm for the architecture

of vector supercomputers. Comparisons with a VAX 11/750 and an IBM 3081

large mainframe are highlighted in table 8.

Finally the network model was solved for a range of relative weights,

with code NLPNETG. We observe (table 9) that a complete analysis can be

carried out within a few hours, even on a VAX minicomputer. Giving more

emphasis to the risk function (nonlinear component) causes the problem to

become, algorithmically, more difficult. This difficulty reflects the

price we have to pay in going from a linear, transportation cost minimizing

model, to a nonlinear model that incorporates the nonlinear form of risk

minimization.
-33-

7. Conclusions and Future Research

We have discussed in this paper a general optimization planning model

for the U.S. air-traffic system. The use of this model has been demonstrated

using air-traffic data from a representative sector. By developing

specialized algorithms, we were able to solve the resulting problem

efficiently, thus demonstrating the feasibility of the approach for strategic

planning. In addition, taking advantage of the latest technological developments in supercomputer design we were.able to solve very large probleems in a

matter of seconds; thus, the network model can be used as an operational

planning tool under real-time conditions.

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 modeling

the air-traffic system. We have considered a deterministic model of risk.

In general a stochastic analysis which also considers the systemwide optimization 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 intervention 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 [15].

Other aspects of the air-traffic control system can be examined within

a network optimization framework, such as the flow control problem [1], or

aircraft scheduling problems arising from aggregate solutions of the current
-34-

model. Issues like the uncertainty involved in determining projected*

demand for air transportation merit investigation, and the end effects

due to the finite planning horizon have to be examined. As another

extension, resource limitations can be imposed by introducing arc

multipliers on traffic arriving/departing from an airport, thus controlling

the total number of passengers an airport facility can handle.

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 and

its results 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 [16].
-35-

fs

POOR QUALITY

PROBLEM

PTN30

PTNJ50

PTN660

SMBANK

B1GBANK

RANBANK

STICK 1

STICK2

STICKS

ST1CK4

GROUPJac

GROUPlad

GROUPlae

MARK]

MARK2

MARKS

FAA.TEST

FAA.JND3

Size

(Nodes/Arcs)

30; 46

150/ 196

666/ 906

64 / 117

1116/2230

800/2000

209/ 454

650/1412

782/1686

832/2264

200/ 500

400/1000

800/2000

23/ 42

53/ 102

857/1710

28/ 64

1295/2873

Free arcs

at optimum

15

44

240

54

946

13

246

763

905

1433

100

119

230

4

10

35

6

787

Condition No.

of Reduced

Hessian

>104

>104

£10*

>10*

>\tf

>10*

£1CP

^lO2

^lO2

>}&

>}&

>104

^lO5

>J02

>102

> l&

>\CP

^JO2

Objective

value

-.3239322E5

-.4819730E5 .

-.2061074E6

-.7129290E7

-.4205693E7

.3018788E6

.6934392E1

.3124563EI

.1117978E2

.1566195E1

.1011792E5

.3834884E5

.22657 14E6

-.1145214E6

-.327J077E6

.1150057E8

.7725644E5

.2556325E1

I,,, Norm

of Reduced

Gradient

.0045

.0020

.0160

.0003

.0004

.4025

.0001

.0001

.0006

.0010

.0060

.0049

.4070

.0020

.2420

.0710

.0500

.0501

Description

Dallas Water

Distribution

models

Matrix balancing models

Stick perculation models,

electrical networks

Randomly generated, strictly

convex networks

Markowjtz

protfolio construction

Air traffic control model

Table 4 : Test Problems

Problem

PTN30

PTN 1 50

PTN660

SMBANK

B1GBANK

RANBANK

STICK 1

ST1CK2

STICKS

ST1CK4

GROUPlae

GROUP lad

GROUPlae

MARK1

MARK2

MARK 3

FAA.TEST

FAA.1ND3

NLPNETG Solution

Time (sec)

7.33

23.86

297.93

21.50

9100.00

245.37

19.05

103.12

73.75

166.50

1652.00

10227.00

7376.00

3.07

21.60

204.43

1.24

707.00

/ ot, norm

4.5E-3

2.0E-3

1.6E-2

3.0E-4

4.0E-4

4.0E-1

l.OE-4

l.OE-4

6.0E-4

l.OE-3

6.0E-3

4.9E-3

4.1E-1

2.0E-3

2.4E-1

7.1E-2

5.0E-2

5.0E-2

MINOS Solution

Time (sec)

44.08

305.01

9647.35

287.20

28800.00*

5477.97*

12392.00

•*

•*

**

2312.69

20347.00

16115.00

12.23

33.03

1341.53

15.22

111600.00

Inflow

8.8E-6

8.8E-6

4.0E-4

4.9E-4

4.0E-1

l.OE+2

3.9E-6

5.9E-2

5.9E2.8E1.3E1.4E1.5E3.4E-9

2.8E-4

Tcnnuotrd wbce wlttuo OO3% frota <n*in»jm

* Did not ceovcrr* •<«" • beun

Table 5 : Comparison of NLPNETG with MINOS
-36-

Problem

PTN30

PTMSO

PTN660

RANBANK

STICK 1

STICK2

STICK3

ST1CK4

GROUPlac

GROUPJad

GROUPlae

MARK!

MARK2

MARK 3

NGSD Solution

Time (sec)

13.25

83.90

590.35

245.37

50.25

517.48

1387.25

2688.72

314.85

( 606,90

395.38

3.07

21.60

204.43

% optimally

0.02

0.50

3.50

4.0E-1

0.40

0.94

0.08

0.15

1.60

0.40

0.52

2.0E-3

2.4E-1

7.1E-2

MINOS Solution

Time (sec)

44.08

305.01

9647.35

5477.97*

12392.00

»•

•*

•*

2312.69

20347.00

16115.00

12.23

33.03

1341.53

la/iorm

8.8E-6

8.8E-6

4.0E-4

l.OE+2

3.9E-6

5.9E-2

5.9E-1

2.8E-1

1.3E-1

1.4E-1

1.5E-1

Terminated wb»n within 0.02% from optimum

* Did DOI cotrvrrtr litti t boun

Table 6 : Comparison of GNSD with MINOS

Problem

PTN 1 50

PTN660

SMBANK

B1GBANK

STICK 4

GROUPlac

GROUPlae

MARK 3

NLPNETG Solution times (sec)

Without vectorization

0.328

2.435

0.358

244.107

*

18.668

•

•

Compiler vectorization

0.317

2.3)6

0.327

225.900

•

17.215

•

•

User vectorization

0.165

1.402

0.177

58.896

2.925

4.983

49.454

2.131

Problem not totvtd with Uu» opfcoo

Table 7 : Vectorization of NLPNETG on the CRAY/XMP
-37-

Problem

PTN150

PTN660

SMBANK

B1GBANK

GROUPlac

GROUP Jac

MARK3

Average

NLPNETG Solution times (sec)

IBM 308 1

1.37

36.68

2.33

884.53

218.46

1320.82

24.87

281.19(17)

VAX 11/750 (Unix)

23.86

297.93

21.50

9100.00

1652.00

10227.00

204.43

3075.25(184)

CRAY/XMP

0.165

1.402

0.177

58.896

4.983

49.454

2.131

16.744(1)

Table 8 : Testing NLPNETG on different computer systems

Weight ( M j )

J.O

0.999999

0.999997

0.999995

0:999992

0.999990

0.999975

0.0

Time (sec)

J2JJ

707

396

396

406

367

' 357

18

/ x norm

0.020

0.054

0.094

0.150

0.187

0.467

0.069

0.0

Table 9 : Tracing the efficient frontier with NLPNETG
-38-

References

[1] Proceedings of the Conference of Safety Issues in Air Traffic Systems

Planning and Design, Princeton University, Princeton, New Jersey,

September, 1983.

[2] D. P. Ahlfeld, R. S. Dembo, J. M. Mulvey and S. A. Zenios, "Nonlinear

Programming on Generalized Networks," Engineering Management Systems

Working Paper EES-85-7, Princeton University, June, 1985.

[3] P. Beck, L. Lasdon, M. Engquist, I. Balberg, N. Binnebaum and

D. P. Bertsekas, "On the Goldsten-Levitin-Polyak Gradient Projection

Method," Physical Review, Vol. 28, No. 7, pp. 3799-3812, 1983.

[4] R. S. Dembo and J. G. Klincewicz, "Resolving Degeneracy Using a

Maximal Basis," School of Organization and Management Working Paper

Series B #63, July, 1983.

15] A. R. Ferguson and G. B. Dantzig, "The Allocation of Aircraft to Routes,"

Management Science, Vol. 3, pp. 45-73, 1957.

[6] A. M. Geoffrion, "Elements of Large-Scale Mathematical Programming,"-

Management Science, Vol. 16, No. "11, pp. 652-691, 1970.

[7] J. L. Helms, "Implementing the National Airspace System Plan," Safety

Issues in Air Traffic Systems Planning and Design, Princeton University

Conference 155, September, 1983.

[8] C. A. Holloway, "An Extension of the Frank-Wolfe Method of Feasible

Directions," Mathematical Programming, Vol. 6, No. 1, pp. 14-27, 1974.

[9] J. L. Kennington and R. V. Helgason, Algorithms for Network Programming,

John Wiley and Sons, New York, 1980.

[10] S. Lawphongpanich and D. W. Hearn, "Simplicial Decomposition of the

Asymmetric Traffic Assignment Problem," Transportation Eesearch,

Vol. 18B, No. 2, pp. 123-133, 1984.

[11] B. McDonald, unpublished report to the FAA, March, 1985.

[12] J. M. Mulvey and S. Z. Zenios "Solving Large Scale Generalized Networks,"

Journal of Information and Optimization Science, Vol. 6, pp. 95-112,

1985.

[13] J. M. Mulvey, S. A. Zenios and D. P. Ahlfeld, "Simplicial Decomposition

for Convex Generalized Networks," Engineering Management Systems

Working Paper EES-85-8, Princeton University, 1985.

[14] B. Murtagh and M. Saunders, "Large-Scale Linearly Constrained Optimization,"

Mathematical Programming, Vol. [bf 14],.pp. 41-72, 1978.
-39-

[15] A. Odoni and S. Endoh, "A General Model for Predicting the Frequency

of Air Conflicts," Safety Issues in Air Traffic Systems Planning and

Design, Princeton University Conference #115, September, 1983.

[16] W. B. Powell, J. M. Mulvey and S. Babu, "Development and Validation

of Models of Air Traffic Control Performance," Report, DOT DTFA

03-82-C-400046, Princeton University, September, 1984.

[17] R. T. Rockafellar, Convex Analysis, Princeton University Press,

Princeton, New Jersey, 1970.

[18] W. D. Rowe, An Anatomy of Risk, John Wiley and Sons, New York, 1977.

[19] B. von Hohenbalken, "A Finite Algorithm to Maximize Certain

Pseudoconcave Functions on Polytopes," Mathematical Programming,

Vol. 13, pp. 189-206, 1975.

[20] B. von Hohenbalken, "Simplicial Decomposition in Nonlinear Programming

Algorithms," Mathematical Programming, Vol. 13, pp. 49-68, 1977.