FINAL REPORT

AN OPTIMIZATION MODEL FOR THE U.S. AIR-TRAFFIC SYSTEM

John M. Mulvey

(NASA-CR-177277) AN OPTIMIZATION MODEL FOR N86-27271

THE US aiR-TRfiPFIC SYS1EM Final Report

(Princeton Univ., N. ,-J.) 27 p HC A03/MF AOI

CSCL 17G Unclas

63/04 43135

School of Engineering and Applied Science

Princeton University

Princeton, New Jersey 08544

April, 1986

Grant No. NAG-1-520

Account No. 165-6210

National Aeronautics and Space Administration
An Optimization Model for the TJ. S. Air-Traffic System

John M. Mulvey

School of Engineering and Applied Science

Princeton University

Princeton, New Jersey 08544

April 1986^

1. Executive Summary

This report summarizes the research carried out as part of the project "An Optimization

Model for the U. S. Air-Traffic System (NASA* NAG-1-520). The primary research objective was to establish a systematic approach for monitoring U. S. air traffic in the context of

system-wide planning and control. Towards this end, a network optimization model with nonlinear objectives was chosen as the central element in the planning/control system. The network representation was selected because: (1) it provides a comprehensive structure for depicting essential aspects of the air traffic system, (2) it can be solved efficiently for large scale

problems, and (3) the design can be easily communicated to non-technical users through computer graphics.

Briefly, the network planning models consider the flow of traffic through a graph as the

basic structure. Nodes depict locations and time periods for either individual planes or for

aggregated groups of airplanes. Arcs define variables as actual airplanes flying through space

or as delays across time periods. As such, a special case of the network can be used to model

the so called flow control problem2.

Due the large number of interacting variables and the difficulty in subdividing the problem into relatively independent subproblems, we designed an integrated model which will

depict the entire high level ( above 29000 feet) jet route system for the 48 contiguous states in

the U. S. While the resulting model is large by today's standards, we felt that computer technology (hardware and software) is improving rapidly and a practical full-scale system should
- 2 -

become feasible over the next few years.

As a first step in demonstrating the concept's feasibility, we gathered data from the Indianapolis control sector and built a nonlinear risk/cost model for this airspace. The nonlinear

network program --NLPNETG-- was employed in solving the resulting test cases1. This

optimization program uses the Truncated-Newton method (quadratic approximation)3 for

determining the search direction at each iteration in the nonlinear algorithm.

It was shown that aircraft could be re-routed in an "optimal" fashion whenever traffic

congestion increased beyond an acceptable level, as measured by the nonlinear risk function.

An efficient cost/risk frontier can be traced out by altering the relative weights between cost

and risk. This curve can be employed within the context of strategic planning, e.g. adding airport resources, or as an operational decision tool for assisting air-traffic controllers.

Computational costs for this exercise were quite reasonable and the nonlinear optimization problems could be solved on an minicomputer ~ VAX 11/750 —at Princeton University.

In addition, we transported the optimization system to a CRAY supercomputer at Boeing

Computer Services, Seattle, Washington. Taking advantage of the CRAY vector architecture

required a substantial effort; see reference6 for a description of the steps that were taken. This

analysis demonstrated that the Trucated-Newton algorithm can be tuned for a vector computer and that substantial air-traffic control problems could be solved in real-time.

The technical details of the project are contained in three research reports. The report

titles are listed below:

(1) "Nonlinear Network Programming on Vector Supercomputers," Report EES-85-13,

Princeton University, submitted to Operations Research, 1986.

(2) "Real-Time Operational Planning for the U.S. Air-Traffic System", Report EES-86-5,

Princeton University, submitted to Applied Numerical Mathematics, 1986.

(3) "Integrated Risk/Cost Planning Models for the U. S. Air Traffic System," Report EES85-9, Princeton University, submitted to Management Science, 1985.
- 3 -

These reports are self-contained, describe the research work performed during the project,

and are provided in the Appendix.

2. Future Research

We believe that future research involving integrated air-traffic modeling should be

directed along several avenues. First, the size of the network test problems should be

expanded by including multiple control sectors and a longer time horizon -- more periods. As

the problems grow in size, every attempt must be made to keep the execution time within a

reasonable level. The runtime issues are especially pertinent in the case of the operational

planning model. Remember that the final goal is to develop a practical decision support system for air traffic planning.

Second, an interface should be designed so that non-technical personnel will be able to

effectively use the planning systems. It will be essential to understand the type and the form of

information which is needed by the air traffic controller, among others. An interactive procedure which employs graphics seems to be the most likely format. Perhaps, color and intensity could be used to identify the critical elements in the control domain. Also, the optimization procedure needs to be flexible enough to be able to handle alternative objective functions, priorities and possible intervention by the users. Criteria other than utilitarian can be

achieved by adjusting the objective function in the network models.

Third, the risk objective function which measures issues such as congestion, proximity,

and workload, must be refined and tested with empirical findings. There is much to be done

in this arena. For example, the analysis must be coordinated with the micro-level studies carried out by Odoni5

 and others. The field of risk analysis is only beginning to develop general

principles. As an approximation, one could use a congestion function which depends upon

geometry, time of day, weather, and other factors to evaluate the degree of acceptability for an

airspace as compared to delaying or re-routing planes anticipated to cross the region of

interest. See the second and third papers listed in the Appendix for further details.
- 4 -

In summary, the National Airspace Plan4 promises to greatly enhance the quality and

the quantity of information available on the state of the U.S. air traffic system. Never fhe less,

much effort will be needed to decide what information is critical to a safe and efficient airtraffic environment. In this regard, planning models will assist in important task of identifying

critical informational flows. Both the strategic and the operational network models proposed

in this project are designed for this task.
- 5 -

References

1. D. P. Ahlfeld, R. S. Dembo, J. M Mulvey, and S. A. Zenios, "Nonlinear Programming on

Generalized Networks," Report EES-85-7, submitted for publication to Transactions on

Mathematical Software, Princeton University, June 1985.

2. M. Bielli, G. Calicchio, B. Nicolette, and S. Ricciavdelli, "The Air-traffic Flow Control

Problem as an Application of Network Theory," Computations and Operations Research,

vol. 9, no. 4, pp. 265-278, 1982.

3. R. S. Dembo, "A Primal Truncated Newton Algorithm for Large-Scale Nonlinear Network Optimization," School of Organization and Management Working Paper, Series

B#72, Yale University, March 1983.

4. J. L. Helms, "Implementing the National Airspace System Plan," in Safety Issues in Air

Traffic Systems Planning and Design, Princeton University Conference, Sept. 1983.

5. A. Odoni and S. Endoh, "A General Model for Predicting the Frequency of Air Conflicts," in Safety Issues in Air Traffic Systems Planning and Design , Princeton University

,Sept. 1983.

6. S.A. Zenios and J.M Mulvey, "Nonlinear Network Programming on Vector Supercomputers," Report EES-85-13, Princeton University, Feb., 1986.
APPENDI X
This appendix contains the computer program used for generating the

Air-Traffic test cases in the project.

c* *

c* *

C''

c

v

c

c PURPOSE :

c This program reads the data for the airports, flights and

fuel burn model and generates the network model for input to the

nonlinear network optimizer NLPNETG

It compiles with F77 on the VAX 11/750 under UNIX 5.2

INPUT :

GENER

BODY

TAIL

LENGTH

c

c

c

c

c SUBROUTINES

c

c

c

c

c

c

c

c INPUT FILES

c

c

c

c

C OUTPUT FILES

c

c

c

c

c

c

c DEFINITION OF VARIABLES

c

CHARACTER Arrays

Reads in the data for Airports/Flights/Fuel burn

: Generates the network topology and identifies risk

: Writes on the output file the generated network

: Writes on the output file the risk function

Calculates the distance traveled by a flight

Logical unit

Logical unit

Logical unit

Logical unit

Logical unit

Logical unit

Fuel burn model data

Airports location data

Flights strips data

10 : Nonlinear network model

15 : Warning/Error messages from the program

16 : Scratch file for storage of interractions

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

c

airprt

fights

(nairpt)

(nflight)

REALM Arrays

burn (nflight*2,9)

cost (narc)

INTEGER Arays

iorig (nflight)

idest (nflight)

ifIght(nflight)

alti (nflight,2)

itime (nflight,2)

otime (nflight,2)

hemi

rate

(nflight)

(nflight)

lati (nairpt,3)

long (nairpt,3)

elev (nairpt)

11

jj

(narc)

(narc)

Four character airports code

Six character flight id code

Fuel burning rates at 9 legitimate cruise altitudes

Cost of generated network arcs

Pointer to origin airport for every flight

Pointer to destination airport

Pointer to fuel burn record

Entering/Exiting altitude

Entering time ( hour/min)

Exiting time ( hour/min)

Hemi code '

Burning rate indicator '

- 1 Lbs/nmi

- 2 Lbs/hr

Latitude of every airport

Longtitude of every airport

Elevation of every airport

From node of generated network arcs

To node of generated network arcs

-1-
c

c istart(9) : Pointer to array ITAG for arcs contributing

c to the risk function at different altitudes

c itag (narc ) : Thread of arcs contributing to risk

c

c... INTEGER Variables

c

c nflight : Number of flights

c itper : Number of time periods considered for flights delay

c time : Length of each delay time period

c inter : Number of time intervals considered for risk analysis

c

C**'

c

c

program faanet

c

c

c

include '../data/comdat'

open (7)

open (8)

open (9)

open (10)

open (15)

open (16)

print *, ' [H [J'

c

c Read Input data

c

write (15,900)

call input

c

c Generate and write out the Network topology

c

write (15,910)

call gener

call body

c

c Determine and write risk interractions over target sector

c

write (15,920)

call risk

call tail

c

print *, ' [H [J #3 #4 D 0 N E ! [2B'

c

c FORMAT Statements

c

900 format (/,lx,'*** READ INPUT DATA ...',/)

910 format (/,lx,'*** GENERATE NETWORK TOPOLOGY ...',/)

920 format (/,lx,'*** DETERMINE FLIGHT INTERRACTIONS ...',/)

c

c

close (7)

close (8)

close (9)

close (15)

close (16)

c

stop

end

-2-
c

subroutine input

c

£* * * * 1

C

•c PURPOSE :

c This subroutine reads the data from the three

c input files as described in the main program .

c

include '../data/comdat'

character*6 aflt

character*4 aorig, adest

c

c Prompt the user for control parameters

c

print *,' [l;4;7m #3 #4 #6 FAA NETWORK GENERATOR [Om [3B'

print *,'Number of time periods considered for flight delay ?'

read (5, *) itper

print *,'Length of flight delay time interval ?'

read (5,*) time

print *,'Number of time periods considered for risk analysis ?'

read (5,*) inter

print *,'Value of printing flag ?'

,read (5,*) mprint

print *,'Header card for generated network ( at most 72 char) ?'

read (5,890) title

c

print *,' [3B [4C [l;4;5;7m #3 #4 #6 R U N N I N G [Om'

c

c Read airports code and location

c

i = 1

10 read (8,900,end=20) airprt(i), lati (i,1),lati(i,2), lati(i,3),

* long(i,l), long(i,2),long(i,3),elev (i)

i = i + 1

go to 10

c

c Read fuel burn data

c

.20 nairpt = i - 1

1 = 1

c

30 read (7,910,end=40)fights(i),hemi(i),rate(i),(burn(i,j),j=l,9)

i = i + 1

go to 30

c

c Read flights data and set pointers to airport and fuel-burn recors

c

40 nfbrec = i - 1

i = 1

c

50 read (9,920,end=100) aflt, aorig, adest ,ialtl, ialt2 ,

* (itimeU/ j) , j=l,2) , (otime(i,j), j=l,2)

c

alti (i,l) = ialtl * 100

alti (i,2) = ialt2 * 100

c

if.\agl = 1

iflag2 = 1

c

c

do 60 j = 1, nairpt

if (aorig.eq.airprt(j)) then

iorig(i) = j

iflagl = 0

end if

-3-
if (adest.eq.airprt(j)) then

idest(i) = j

iflag2 = 0

end if

if (iflagl.eq.0.and.iflag2.eq.O) go to 70

60 continue

c

if (iflagl.ne.O) write (15,930) aorig , aflt

if (iflag2.ne.O) write (15,935) adest , aflt

c

70 do 80 j = 1, nfbrec

if (aflt.eq.flghts(j) ) then

iflght(i) = j

go to 90

end if

80 continue

c

write (15,940) aflt

c

90 i = i + 1

go to 50

100 nflight = i - 1

c

if (mprint.le.3) go to 999

write (15,950)

write (15,960)

write (15,970) (i,iorig(i),idest(i),

* (alti(i,j),itime(i,j),otime(i,j),j=l,2), i=l,nflight)

write (15,980)

write (15,990) (airprt(i), lati (i,1),lati(i,2),lati(i,3),

* long(i,1),long(i,2),long(i,3), elev(i), i=l,nairpt)

write (15,992)

write (15,994) (i,fights(i),ifIght(i),hemi(i),rate(i) ,

* burn(i,1),i=l,nfbrec)

c

c FORMAT Statements '

c

890 format (a72)

900 format (Ix,a4,2(5x,3i5),5x,15)

910 format (Ix,a6,14x,il, Ix,il,9(lx,ell.0))

920 format (Ix,a6,6x,a4,Ix,a4,5x,2(2x,i3),2(3x,i2),5x,2(3x,i2))

930 format (lx,'No information available for origin airport ',a4,

* ' of flight ',a6)

935 format (Ix,'No information available for destination airport ',

* a4,' of flight ',a6)

940 format (lx,'No information available for fuel burn rates of flight'

* Ix,a6)

950 format (Ix,' OUTPUT from subroutine INPUT ',/)

960 format (/,Ix,'No.',2x,'Orig',2x,'Dest',3x,'Altitude'

* ,2x,'In-tm',2x,'0-tm',/)

970 format (Ix, i3, Ix, i5, Ix, i5, Ix, i4,2x, i4, Ix, i2,2x, i2, Ix, i2, 2x, i2)

980 format (/,Ix,'Airprt',6x,'Latitude',8x,'Longitude',4x,'Elevation',/)

990 format (2x,a4,2x,3i5,2x,3i5,3x,i5)

992 format (Ix,'No.',Ix,' Fight',2x,'IfIght',Ix,'Hemi',Ix,'Rate',

* 5x,'Fburn',/)

994 format (Ix,i3,Ix,a6,2x,i5,3x,il,4x,il,2x,e!5.6)

c

c

999 return

end

c

C*

c

block data

c

-A-
include '../data/comdat'

data ihO /29000,33000,37000,41000,45000/

data ihl /31000,35000,39000,43000/

data nfaa/0/

end

-5-
subroutine gener

c

c PURPOSE

c

c

c

c

c

c

c

c INPUT

c

c

c

c

c OUTPUT

c

c

c

c

c

c

c

c

c

c

This subroutine generates the network topology for

the Flight/Airport information given for one control sector.

It generates alternative cruise altitudes for every flight

and determines associated flight costs as well as delay costs .

For every generated route (time/altitude) its contribution

to the risk function is determined .

The arrays read from the database files by routine INPUT

ihO(nflight) , legitimate cruise altitudes for hemi-code = 0

ihl(nflight) , legitimate cruise altitudes for hemi-code = 1

ii (narc)

jj (narc)

cost (narc)

ipntfl (narc)

iper (narc)

istart (inter*9)

itag (narc)

narc

from node of generated arc

to node of generated arc

cost of generated arc

flight number of this arc

time period of this arc

pointer to ITAG array

thread of Ifights contributing to risk

total number of generated arcs

include '../data/comdat'

narc

i

j

c

c=

c=

c=

c

c

c

Consider origin nodes for all time periods.

do 1000 1=1 , itper

do 500 j = 1, nflight

iinod = (i-1)*nflight + j

jjl = i*nflight + j

jj2 = (itper+i-1)*nflight + j

itin = itime(j,l)*60 + itime(j,2) + (i-l)*time

itout = otime(j,1)*60 + otime(j,2) + (i-1)*time

c

c Takeoff delay arcs

c

c

c En route arcs

c

ipnt = iflght (j) +1

if (rate(ipnt).ne.2) write (15,900) fights(ipnt)

costd = time * burn(ipnt,1) / 60.0

ipnt = iflght (j)

ihcode= hemi (ipnt)

ialtO = 0

ialtl = alti (j,l)

ialt2 = 0

-6-
c determine alternative cruise altitudes for this flight

c

if (ihcode.eq.1) go to 50

. alternatve routes for Hemi code 0 flights

if (alti(j,l) .eq. alti(j,2)) go to 30

ialtl = ( alti(j,l) + alti(j,2)) / 2

c

do 20 k = 1, 4

if (ialtl.It.ihO(k+l).and.ialtl.gt.ihO(k)) then

ialtl = ihO(k+l)

go to 30

end if

20 continue

30 continue

c

do 40 k = 1, 5

if (ialtl.eq.ihO(k)) go to 45

40 continue

c

45 fburnl = burn (ipnt,2*k-l)

if (fburnl.eq.0.0) write (15, 910) fights(ipnt)

c

if (k.gt.l) then

ialtO = ihO(k-l)

fburnO= burn (ipnt,2*k-3)

if (fburnO.eq.0.0) fburnO = fburnl

end if

c

if (k.lt.5) then

ialt2 = ihO(k+l)

fburn2= burn (ipnt,2*k+l)

if (fburn2.eq.0.0) fburn2 = fburnl

end if

c

go to 100

c

c

c alternatve routes for Hemi code 1 flights

c

50 continue

if (alti(j,l).eq.alti(j,2)) go to 80

ialtl = ( alti(j,l) + alti(j,2)) / 2

c

do 70 k = 1, 3

if (ialtl.It.ihl(k+1).and.ialtl.gt.ihl(k)) then

ialtl = ihl (k+1)

go to 80

endif

70 continue

80 continue

c

do 90 k = 1, 4

if (ialtl.eq.ihl(k)) go to 95

90 continue

c

95 fburnl = burn(ipnt,2*k)

if (fburnl.eq.0.0) write (15, 910) fights(ipnt)

c

if (k.gt.l) then

ialtO = ihl(k-l)

fburnO= burn (ipnt,2*k-2)

if (fburnO.eq.0.0) fburnO = fburnl

end if

-7-
if (k.lt.4) then

,ialt2

fburn2= burn (ipnt,2*k+2)

if (fburn2.eq.O.0) fburn2 = fburnl

end if

c

c determine the cost of alternative altitudes

c

100 continue

call length (j ,dist)

c

costO = fburnO * dist

costl = fburnl * dist

cost2 = fburn2 * dist

c

c== Store generated arcs. (For the last time period no delays possible)

c

c > Save Delay arcs

c

if (i.ne.itper) then

narc = narc + 1

ii (narc) = iinod

jj (narc) = jjl

cost(narc)= costd

ipntfl(narc) = j

iper (narc) = i

end if

c

c > Save flight arcs on different altitudes

c

C PRIMARY-1

if (ialtO.gt.O) then

narc = narc + 1

ii (narc) = iinod

jj (narc) = jj2

cost(narc)= costO

ipntf1(narc) = j

iper (narc) = i

end if

c

c PRIMARY

narc = narc + 1

ii (narc) = iinod

jj (narc) = jj2

cost(narc)= costl

ipntf1(narc) = j

iper (narc) = i

c

c PRIMARY+1

if (ialt2.gt.O) then

narc = narc + 1

ii (narc) = iinod

jj (narc) = jj2

cost(narc)= cost2

ipntf1(narc) = j

iper (narc) = i

end if

c

500 continue

1000 continue

c

c== Consider destination nodes for all time periods ==

-8-
itor = itper * nflight

itort2 = itor * 2.0

c

do 2000 i = 1, itper

c

do 1500 j = 1, nflight

c

iinod = itor + (i-1) *nflight + j

jjl = itor + i *nflight + j

jj2 = itort2 + j

c

c --- Landing delay arcs . . .

c

ipnt = iflght (j) +1

if ( rate (ipnt) .ne. 2) write (15,900) flghts(ipnt)

costd= time * burn ( ipnt, 1) / 60.0

c

c --- Store generated arcs. (For the last time period no delays possible)

c

if (i.ne. itper) then

narc = narc +1

ii(narc) = iinod

jj (narc) = jjl

cost (narc) = costd .

ipntfl(narc) = j

iper (narc) = i

end if

c

narc = narc + 1

ii(narc) = iinod

j j (narc) = j j2

cost (narc) = 0.0

ipntfl(narc) = j

iper (narc) = i

c

1500 continue

2000 continue

c

c --- FORMAT Statements -------------------------------------------------

c

900 format (Ix, 'Check fuel burn data for flight ',a6)

910 format (Ix, 'Fuel burn data for flight ',a6,' is not available')

c

return

end

c

subroutine length (iflt, dist)

c

c

c --- PURPOSE :

c Returns the distance between two airports - origin and

c destination - of a flight.

c NOTE : The distance is measured ignoring difference in

c the elevation of the origin-destination airports and the cruise

c altitude of the flight. The earth curvature is taken into account

c

c --- INPUT :

c iflt : Flight number indicator

c airprt (nairpt) : List of airport codes

c lati (nairpt, 3) : Latitude of all airports

c long (nairpt, 3) : Logitude of all airports

-9-
c

c

c

c

c

c

iorig (nflght)

idest (nflght)

degrad

radmil

Pointer to origin airport for flight IFLT

Pointer to destination airport for flight IFLT

PI/180 - converts degrees into radians

: 24844mls /2*PI*1.151 -converts rad into naut. miles

include '../data/comdat'

real latl, Iat2

c

c Parameter initialization

c

c

cc

c

c

cc

parameter ( degrad = 0.017453 , radmil = 3439.4 )

ipo = iorig (iflt)

ipd = idest (iflt)

Calculate latitudes of Origin/Destination airports in radians

Compute the linear distance between the airports as :

(GREAT CIRCLE DISTANCE) * (EARTH CIRCUITY FACTOR)

xl = lati(ipo,l) + lati(ipo,2) / 60.0 + lati(ipo,3) / 3600.0

yl = long(ipo,1) + long(ipo,2) / 60.0 + long(ipo,3) / 3600.0

x2 = latidpd, 1) + lati(ipd, 2) / 60.0 + lati(ipd,3) / 3600.0

y2 = long(ipd,l) + long(ipd,2) / 60.0 + long(ipd,3) / 3600.0

latl = degrad * yl

Iat2 = degrad * y2

dist = acos(sin(latl)*sin(lat2) + cos(latl)*cos(Iat2)*

cos(degrad*(xl-x2)))*radmil

if (mprint.ge.3) write (15,900) iflt, dist

c

c

c

900

c

c

FORMAT Statements

format (Ix,'Length of flight ',i3,' is ',e25.12)

return

end

-10-
c *:

c

subroutine body

c

c PURPOSE :

c This subroutine writes on the output file the body

c of the generated network.

c

c INPUT :

c The arrays generated by routine GENER

c

c

include '../data/comdat'

c

c Initialization

c

character*! a, s, d

character*2 bd

parameter (iO=0, il=l, i2=2, rO=0, rl=l)

data a/lha/, s/lhs/, d/lhd/, bd/2hbd/

. c

c Check the size of generated network

c

nodes = 2*nflight*itper + nflight

na = 3*nflight*itper + 2*nflight*(itper-1) + nflight*itper

nn = jj (narc)

c

if (na.ne.narc) write (15,900) na, narc

if (nn.ne.nodes) write (15,910) nodes, nn

c

c Write header card

c

write (10,920) title

write (10,930) nn, narc, i2, il

c

c Write arc data

c

write (10,940) a

do 100 i = 1, narc

write (10,950) bd, iO, ii(i), jj(i),rO,rl,rO,cost(i),rO,rO

100 continue

c

c Write supply data

c

write (10,940) s

do 200 i = 1, nflight

write (10,970) i , rl

200 continue

c

c Write demand data

c

write (10,940) d

itemp = 2*nflight*itper

do 300 i = 1, nflight

nod = itemp + i

write (10,970) nod, rl

300 continue

c

c FORMAT Statements

c

900 format (lx,'WARNING : Model should have ',i5,' arcs. Only ',i5,

* ' were generated.',/,

*llx,'Data for some flights may be missing. Check the databases',/)

-11-
910 format (lx,fERROR : Model should have ',13,' nodes. Only ',13,

* ' were generated.',/,

*9x,'Possible program error.',/)

920 format (Ix,a72)

930. format (415)

940 format (al,t72,' ')

950 format-(3x,a2,3x,12,215,3f10.1, 3f10 .1)

970 format (lOx,15,5x,f10.0,t72,' ') .

c

return

end

c

c

subroutine tail

*****

c

c PURPOSE :

c This subroutine writes on the output file the components

c of the risk function (ie. nonlinear and non-separable function).

c For the special case when only a single flight is crossing

c the sector the risk function is separable.

c

c INPUT : The arrays generated by routine RISK

.c

c

include '../data/comdat'

c

character*3 f , e

parameter (zero=0.0, two=2.0)

data f/3hfaa/, e/3hend/

imax = - 1

imin = 99999

rmxcst= 0.0

c

c > Determine the max cost to use as scaling factor.

c

do 5 i = 1, narc

costi = cost (i)

if (rmxcst.It.costi ) rmxcst = costi

5 continue

rmxcst = rmxcst / 10.0

c

c=== Determine the planning interval

c

do 10 i = 1, nflight

itmax = otime(i,1)*60 + otime(i,2)'

itmin = itime(i,1)*60 + itime(i,2)

if (itmax.gt.imax) imax = itmax

if (itmin.lt.imin) imin = itmin

10 continue

c

c > Length of interval used for risk evaluation

c

rintl = real (imax-imin) / real (inter)

write (15,900) imin, imax, rintl

c

c=== Consider all flights at the same cruise altitude

c

write (10,910) f

do 1000 i = 1, 9

-12-
ipi = istart(i)

if (ipi.eq.0.or.itag(ipi).eq.0) go to 1000

c

100. if (ipi.eq.O) go to 1000

ipj = ipi

200 ipj = itag (ipj)

if (ipj.eq.O) go to 500

c

iarci = ipntfl (ipi)

iarcj = ipntfl (ipj)

c

if (iarci.eq.iarcj) go to 200

c

iaddi = time * (iper(ipi) - 1)

iaddj = time * (iper(ipj) - 1)

isti = itime(iarci,1)*60 + itime(iarci,2) +.iaddi

iendi = otime(iarci,1)*60 + otime(iarci,2) + iaddi

istj = itime(iarcj,1)*60 + itime(iarcj,2) + iaddj

iendj = otime(iarcj,1)*60 + otime(iarcj,2) + iaddj

c

ist = maxO (isti ,istj )

iend = minO (iendi,iendj)

if (ist.ge.iend) go to 200

c

c > ALPHA is proportional to the overlapping time in sector

c

rtemp = iend - ist

alpha = rmxcst * (rtemp / rintl)

write (16,920) ipi, ipj, alpha, two, zero

go to 200

c

500 ipi = itag (ipi)

go to 100

c

1000 continue

c

write (16,920) narc, narc, zero, zero, zero

close (10)

call system ('sort +10 -15 < fort.16 » fort.10')

c write (10,910) e

if(mprint.ge.2) write(15,930)

* (istart(i),1=1,9) , (itag(i),i=l,narc)

c

c FORMATS

900 format (Ix,'Planning interval: From ',i3,'min to ',i4,

* lx,'min. Length of risk interval:',f6.2)

910 format (a3,t25,'1.0',t72,' ')

920 format (10x,2i5,3f10.4)

930 format (Ix,'ISTART(i)',9(lx,i5),/,Ix,'ITAG(i)',/,100(10(Ix,i5),/))

c

return

end

-13-
c

subroutine risk (num, ist)

c

c

c PURPOSE :

c This subrutine generates the arrays for the nonlinear

c part of the network model. For a given flight/time it identifies the

c altitude and time interval when a conflict may arise and creates the

c thread that links all interracting flights

c

c INPUTS : .

c num , number of risk periods while flight is in target sector

c ist , pointer to the flight altitude while in target sector

c

c OUTPUTS :

c istart, pointer to the starting location of the thread of

c interracting flights

c itag , thread of interracting flights

c

cinclude '../data/comdat'

c

c=== Detemine the flights that cross the target sector at the same altitude ===

c

c > Initialize storage arrays

c

do 12 i = 1, 9

istart (i) = 0

12 continue

do 15 i = 1, narc

itag (i) = 0

15 continue

c

do 500 i = 1, itper

do 200 j = 1, nflight

c

ipnt = iflght (j)

ihcode= hemi (ipnt)

ialtO = 0

ialtl = alti (j,l)

ialt2 = 0

c

iinod = (i-l)*nflight + j

jj2 = (itper+i-l)*nflight + j

c

c > Find all arcs for this flight

c

icount = 0

18 icount = icount + 1

if (ii(icount).eq.iinod.and.jj(icount).eq.jj2) goto 19

if (icount.It.narc) go to 18

write (15,910) fights(j), i

go to 200

c

19 continue

c

c > Determine the flight altitude

c

if (ihcode.eq.1) go to 50

c

c... Altitude for Hemi code 0 flights

c

-14-
if (alti(j,l).eq.alti(j,2)) go to 30

ialtl = (alti(j,l) + alti(j,2)) / 2

c

do 20 k = 1, 4

if (ialtl.It.ihO(k+l).and.ialtl.gt.ihO(k)) then

ialtl = ihO(k+l)

go to 30

end if

20 continue

30 continue

c

do 40 k = 1, 5

if (ialtl.eq.ihO(k)) go to 45

40 continue

45 continue

c

c > Save all the arcs corresponding to this flight at all altitudes

c

c PRIMARY - 1

if (k.ne.1.and.ii(icount).eq.iinod.and.jj(icount).eq.jj2) then

kml = k - 1

call save (kml, icount)

icount = icount + 1

end if

C PRIMARY

if (ii(icount).eq.iinod.and.jj(icount).eq.jj2) then

call save (k,icount)

icount = icount + 1

end if

c PRIMARY + 1

if (k.ne.5.and.ii(icount).eq.iinod.and.jj(icount).eq.jj2) then

kpl = k + 1

call save (kpl, icount)

end if

c

go to 200

c

c... Altitude for Hemi code 1 flights

c

50 continue

if (alti(j,l).eq.alti(j,2)) goto 80

ialtl = (alti(j,l) + alti(j,2)) / 2

c

do 70 k = 1, 3

if (ialtl.It.ihl(k+l).and.ialtl.gt.ihl(k)) then

ialtl = ihl(k+l)

go to 80

end if

70 continue

80 continue

c

do 90 k = 1, 4

if (ialtl.eq.ihl(k)) go to 95

90 continue

95 continue

c

c > Save all the arcs corresponding to this flight at all altitudes

c

c PRIMARY - 1

if (k.ne.1.and.ii(icount).eq.iinod.and.jj(icount).eq.jj2) then

kml = (k - 1) + 5

call save (kml, icount)

icount = icount + 1

end if

c PRIMARY

if (ii(icount).eq.iinod.and.jj(icount).eq.jj2) then

-15-
kp5 = k + 5

call save (kp5,icount)

icount = icount + 1

end if

c— T PRIMARY + 1 ---------------------------------- ------------------

if (k.ne . 4 .and. ii (icount) .eq.iinod.and. jj (icount) .eq. jj2) then

kpl = (k + 1) + 5

call save (kpl, icount)

end if

c

200 continue

500 continue

c

c --- FORMATS ----------------------------------------------------------

905 format (Ix, 'Flight ',a6,' at time period', i3,

* ' is in target sector. Time in:',i3,' Time out:',i3)

910 format <lx,'*** WARNING *** No arcs were generated for flight '

* a6,' at time period ',i5)

c

return

end

c

subroutine save (ipos, iarcn)

c

c

c --- PURPOSE :

c This subroutine stores the specified flight

c in the array of all flights that are crossing the target sector

c simultaneoulsy and at the same level

c

c --- INPUTS :

c iarcn , arc number for flight (time/alittude) considered

c ipos , pointer to array ISTART indicating altitude

c— OUTPUTS :

c istart , pointer to thread storing all flights at this level

c itag , thread of all flights crossing target sector at the

c same interval

c

include ' . . /data/comdat'

c

if (mprint .ge.3) write (15,900) iarcn, ipos

c

c --- Is this the first flight at this altitude ?

c

ip = istart (ipos)

if (ip.eq.O) then

istart (ipos) = iarcn

go to 9999

end if

ipml = ip

c -

c --- If more flights at this altitude exist find an empty spot to store this

c

100 ip = itag (ipml)

if (ip.eq.O) then

itag (ipml) = iarcn

go to 9999

end if

ipml = ip

go to 100

c

9999 if (mprint .ge.3) then

write (15,910) (istart (i) , i=l, 9)

-16-
write (15,920) (itag(i), i=l,narc)

end if

c

c FORMATS

90G format (lx,'Save arc No. ',L5,' in altitude position ',i3)

910 format (Ix,'ISTART(i) :',9(lx,i5))

920 format (lx,'ITAG(i) :',/,100(10(Ix,i5),/))

c

return

end

-17-
character*28 char(l)

character*48 char2(1), char3(l)

10 read (7,100,end=500) char(l)

read (8,110) char2(l)

read (9,110) char3(l)

write (10,200) char(1) (2 :8), char(1) (9 :12), char(1) (13:14)

* ,char(l) (15:15),char(l) (16 :16),char(1) (17:29)

* ,char2(l) (1:13),char2(1) (13:25),char2(1) (25 :37),char2(1) (37:49)

* ,char3(l) (1:13),char3 (1) (13:25),char3(1) (25 :37),charS(1) (37:49)

100 format (a28).

110 format (a48)

200 format (Ix,a6, 3x, 2x,a3,3x,a2,2(lx,al),9(a!2))

go to 10

c

500 print *, 'Finished'

stop

end

-18-
10

100

200

c

500

character*21 char(l)

read (9,100,end=500) char(l)

write (10,200) char<1)(1:4), char(1){5:6), char(l){7:8)

,char(l)(9:10),char(l)(11:13),char(1)(14:15)

,char(l)(16:17),char(l)(18:21)

format (a21)

format (Ix,a4,5x,3(3x,a2),5x,2x,a3,2(3x, a2),5x, Ix,a4)

go to 10

print

stop

end

'Finished'

-19-
character*29 char(l)

10 read (9,100,end=500) char(l)

write (10,200) char(1) (1:6), char(1) (8 :11), char(1) (12:15)

* ,char(l) (16:18),char(1) (19:21),char(l) (22:23)

- * ,char(l) (24:25),char(l) (26:27),char(l) (28:29)

100 format (a29)

200 format (Ix,a6,5x,2(Ix,a4),5x,2(2x,a3),2(3x,a2),5x, 2(3x,a2))

go to 10

c

•500 print *, 'Finished'

stop

end

-20-