# Macroscopic Model and Simulation Analysis of Air Traffic Flow in Airport Terminal Area

Authors: Honghai Zhang; Yan Xu; Lei Yang; Hao Liu
Journal: Discrete Dynamics in Nature and Society
 (2014)

---

## Abstract

We focus on the spatiotemporal characteristics and their evolvement law of the air traffic flow in airport terminal area to provide scientific basis for optimizing flight control processes and alleviating severe air traffic conditions. Methods in this work combine mathematical derivation and simulation analysis. Based on cell transmission model the macroscopic models of arrival and departure air traffic flow in terminal area are established. Meanwhile, the interrelationship and influential factors of the three characteristic parameters as traffic flux, density, and velocity are presented. Then according to such models, the macro emergence of traffic flow evolution is emulated with the NetLogo simulation platform, and the correlativity of basic traffic flow parameters is deduced and verified by means of sensitivity analysis. The results suggest that there are remarkable relations among the three characteristic parameters of the air traffic flow in terminal area. Moreover, such relationships evolve distinctly with the flight procedures, control separations, and ATC strategies.


---

## Body

## 1. Introduction

Air traffic management in terminal area is a knotty problem for controllers as this place is considered to be the air congestion, flight delay, and aviation accident-prone area. Researching on the fundamental operating features of terminal area traffic flow and deducing the macro emergence of traffic flow evolution may contribute to revealing the parameters with mutual relations, as well as the mechanism of spatiotemporal evolution, in terms of the traffic flow characteristic elements. By these means, we can move forward to exploring the objective law in air traffic in order to enrich the air traffic flow theory and to provide scientific basis for air traffic dispersion, which may have very important theoretical value and practical significance.Traffic flow parameters are the physical variables that represent traffic flow characteristics including qualitative and quantitative features of operating states [1]. Basic theories of vehicle traffic flow have developed for decades and many results have been made by scholars. Lighthill and Whitham proposed the simulated dynamic model of traffic flow after researching on the evolution pattern of traffic flow under high traffic density circumstances [2]. Meanwhile, Richards proposed a first order continuum model of traffic flow, which has been integrated as the LWR theory [3]. Biham studied urban traffic flow based on a two-dimensional cellular automaton [4], while Daganzo researched into dynamic traffic problems with a cellular transmission model [5–8]. Compared to the vehicle traffic, less research has been devoted to air traffic flow theory so far, not to mention that many studies focusing just on modeling. A simplified Eulerian network model of air traffic flow was proposed by Menon et al. [9, 10], and Bayen et al. studied the liner control problems derived from Eulerian network model [11–13]. Laudeman et al. noted a quantitative mathematical model on dynamic density of air traffic flow [14]. Complexity model based on traffic flow disturbance was advocated by Lee et al. [15]. Liu et al. proposed a one-dimensional cellular transmission model specifically applicable to air route [16]. Wang et al. studied the microscopic plane-following performance and built the air freeway flow model [17]. Primary discussion for the stability of air traffic flow operating system was advocated by Zhang and Wang, while some basic characteristics of air traffic flow were also involved in this paper [18]. Such research findings have made great foundations for further study on air traffic flow theory. However, no detailed studies have investigated the characteristic parameters with their objective evolution law of air traffic flow. In this paper, we will combine traditional mathematical formula derivation with modern simulation techniques. A macroscopic model of air traffic flow in airport terminal area will be proposed, with which we will simulate and analyze the interrelation and influential evolvement law of the characteristic parameters with the goal of providing theoretical basis for scientific air traffic management.


## 2. Macroscopic Model

### 2.1. Definitions

There are many kinds of definitions for air traffic flow parameters, because of different research intentions and different methods that can be used. Since we focus on airport terminal area, the definitions of velocity, density, and traffic flux on a segment of air route are as follows.Traffic flux (q) is the number of aircraft passing a reference profile of the observation segment per unit of time. q
=
N
/
T, where T is the observation time and N represents the number of aircraft passing in T.Density (ρ) is the number of aircraft per unit length of the observation segment. ρ
=
N
/
P, where P means length of the observation segment, while N is the number of aircraft in P.Velocity can be divided into micro and macro definitions. Microscopic definitions including instantaneous velocity (v) and average velocity (v
¯) focus on some point or profile of the observation segment; macroscopic ones including the space mean velocity (v
¯
s) and time mean velocity (v
¯
t) focus on some area extents or time ranges:
(1)
v
=
d
x
d
t
,
v
¯
=
1
N
∑
i
=
1
N
v
i
,
v
¯
s
=
D
·
(
1
N
∑
i
=
1
N
t
i
)
-
1
,
v
¯
t
=
1
N
∑
i
=
1
N
s
i
t
1
-
t
0
.One flight (f
i) of all (N) uses t
i time to pass the observation segment, where the actual distance is D. And s
i is the flying distance of f
i in the time period t
1
-
t
0.Historically, the first macroscopic traffic flow model is a continuity equation, called the Lighthill-Whitham-Richards (LWR) equation [2, 3]. Like the vehicle traffic flow, air traffic flux (q) and linear density (ρ) may also satisfy the corresponding relations for some given functions f
(
·
) and g
(
·
) with location x and time t, as follows:
(2)
ρ
(
x
,
t
)
=
f
(
q
(
x
,
t
)
,
x
)
,
q
(
x
,
t
)
=
g
(
ρ
(
x
,
t
)
,
x
)
,(3)
∂
q
(
x
,
t
)
∂
x
+
∂
ρ
(
x
,
t
)
∂
t
=
s
.On the right side of (3), s denotes the number of aircraft that enter or exit the observation segment. With the basic equation q
=
ρ
v, where v denotes space mean velocity v
¯
s and assuming every flight has the same flying case, we can derive
(4)
v
(
x
)
∂
q
(
x
,
t
)
∂
x
+
∂
q
(
x
,
t
)
∂
t
=
v
(
x
)
s
.There are two variables as density and velocity but only one equation in the LWR model. It cannot be solved as the differential equation is not closed. In response to this problem, a balancing velocity-density functional relationship, asv
(
x
)
=
v
e
(
ρ
(
x
,
t
)
), was introduced into LWR theory by assuming traffic flow invariably in equilibrium state. To plug this into the equation, a hyperbolic equation of density can be derived, as follows:
(5)
∂
q
e
(
ρ
(
x
,
t
)
)
∂
x
+
∂
ρ
(
x
,
t
)
∂
t
=
s
.


### 2.2. Arrival Traffic Flow Model

LWR model described the propagation characteristic of the nonlinear density wave to find the evolution rule of traffic shock wave and rarefaction wave by characteristics method or numerical simulation [19]. However, air traffic flow differs from vehicle traffic, especially in airport terminal area (Figure 1). First, the arrival and departure traffic in terminal area will follow the designed STAR/SID (standard terminal arrival route/standard instrument departure), which means that aircraft in every route position must act in accordance with the operational flight program, that is, flying within a certain scope of designed flight level and speed. It is usually a small scope and can be reassigned by controllers. As a result of that, some balancing velocity-density functional relationships in vehicle traffic may not exist in air traffic flow. Second, the density of air traffic flow is generally much lower than the vehicle traffic in real operations, which causes the interactions among aircraft to be weaker compared to vehicles. Therefore, it will be difficult to make statistical fit of the relationships according to the available radar data.Figure 1
Radar track plot in airport terminal area (source: ZGGG TMA).To solve the partial differential equation of the LWR theory, we use cell transmission model to discretize the continuity equation of macroscopic traffic flow. The method of discretization is applying a series of interconnected one-dimensional cells to denote the air route and using the difference equation of time discretization to describe aircraft passing through every adjoining cell. It should be feasible to model any air traffic flow scenarios by using interconnected cells, as shown in Figure2, putting all the aircraft flying along the arrival route while in several flight levels onto the same imaginary plane, hence the air traffic flow in one route can be approximately seen as one-dimensional continuous flow. In addition, one air route will be divided into several segments by series of unit cells. To simplify the model we assume that there is a unified variable which we will introduce later to represent different kinds of control measures in one cell, such as speed control, maneuvering actions, or circling, that is, adjusting flows by changing speed or flight path of certain aircraft in the cell.Figure 2
One-dimensional cell transmission model of air traffic flow in single direction.In the matter of arrival route, letN
i
a
(
t
) be the number of arrival aircraft in cell i at time t; then the changes in the number of arrival aircraft in one cell can be described by the following difference equation of time discretization:
(6)
N
i
a
(
t
+
1
)
=
N
i
a
(
t
)
+
τ
i
[
q
i
-
1
a
(
t
)
-
q
i
a
(
t
)
]
.In the equation above,N
i
a
(
t
+
1
) is the number of arrival aircraft in cell i at time t
+
1, q
i
-
1
a
(
t
) represents the flux of arrival aircraft entering cell i from cell i
-
1 at time t, while q
i
a
(
t
) means the flux of arrival aircraft exiting cell i at time t, and τ
i is the time step.Note. The number, flux et al. mentioned in this section, is specific to arrival aircraft, while departures will be included in following parts.Air traffic flow in air route segment still satisfies the basic equationq
i
=
ρ
i
v
i, where ρ
i and v
i stand, respectively, for segment liner density and space mean velocity of cell i. The traffic flux of cell i will be q
i
(
t
)
=
α
i
q
i. Coefficient α
i is the rate of outflow per unit time, which reflects the saturation level of cell i. From Section 2.1., we know segment liner density ρ
i
=
N
i
/
Ω
i, where Ω
i is the length of cell i. Since most of the time arrival aircraft in airport terminal area are in a deceleration process, we assume aircraft entering cell i from some certain arrival route positions with an initial velocity v
i
′ that comes from the STAR and then uniformly decelerating along the cell (route) with a rate a
i
a, so we can get the space mean velocity of cell i, as follows:
(7)
v
i
=
a
i
a
Ω
i
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
.To plug (7) into q
i
=
ρ
i
v
i, we can derive
(8)
q
i
a
(
t
)
=
α
i
a
i
a
N
i
a
(
t
)
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
.However, since the aforementioned situation is an ideal condition, it is necessary to take control measures for part of the arrival aircraft into consideration, such as speed control, maneuvering actions, or circling, because of the existence of traffic congestion, safety interval, and so forth in real conditions. For simplicity, we assume that there is a unified variable to stand for different kinds of control measures, since whichever measures controllers took, the actual effects caused by changing speed or flight path and so forth can all be seen as there to beN
i
a
ATC aircraft flying along the cell (route) with a new constant velocity v
i
ATC in cell i. To be clear, the variable v
i
ATC here represents the displacement velocity that is lower than the STAR designed velocity generally.Considering these problems, we derived (9) after integrating control measures into (8), as follows:
(9)
q
i
a
(
t
)
=
α
i
a
i
a
[
N
i
a
(
t
)
-
N
i
a
ATC
(
t
)
]
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
+
v
i
ATC
N
i
a
ATC
(
t
)
Ω
i
=
α
i
a
i
a
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
N
i
a
(
t
)
-
(
α
i
a
i
a
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
-
v
i
ATC
Ω
i
)
N
i
a
ATC
(
t
)
.After substitutingq
i
a
(
t
) in the difference equation (6) with (9), we got
(10)
N
i
a
(
t
+
1
)
=
(
1
-
α
i
a
i
a
τ
i
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
N
i
a
(
t
)
+
(
α
i
a
i
a
τ
i
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
-
v
i
ATC
τ
i
Ω
i
)
N
i
a
ATC
(
t
)
+
τ
i
q
i
-
1
a
(
t
)
.LetA
i
a
=
1
-
(
α
i
a
i
a
τ
i
/
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
); B
i
a
=
(
α
i
a
i
a
/
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
)
-
(
v
i
ATC
/
Ω
i
); C
i
a
=
α
i
a
i
a
/
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
); we can have the simplified equation sets, as follows:
(11)
N
i
a
(
t
+
1
)
=
A
i
a
N
i
a
(
t
)
+
B
i
a
τ
i
N
i
a
ATC
(
t
)
+
τ
i
q
i
-
1
a
(
t
)
,
q
i
a
(
t
)
=
C
i
a
N
i
a
(
t
)
-
B
i
a
N
i
a
ATC
(
t
)
.The aircraft in cell are taken as evenly distributed in the macroscopic traffic flow theory, so the number of aircraft should be equivalent number. Letd
i
a
(
t
) be the mean nose interval of adjacent aircraft at the same direction in cell i; therefore the number of aircraft in the cell is N
i
(
t
)
=
Ω
i
/
d
i
a
(
t
). If d
i
a
(
t
) becomes less than the control separation standard denoted by d
i
a
ATC, that is to say, aircraft density exceeds the threshold level, the exceeded aircraft should be arranged to take control measures. In addition, arrival aircraft in airport terminal area also have to obey the safety separation d
i
safe according to STAR; that is, d
i
a
(
t
)
≥
d
i
safe. The safety separation is evaluated by air traffic security department and usually gets more stringent compared to control separation as d
i
safe
<
d
i
a
ATC. The number of controlled aircraft at time t in cell i satisfies
(12)
N
i
a
ATC
(
t
)
=
{
(
1
d
i
a
(
t
)
-
1
d
i
a
ATC
)
Ω
i
d
i
a
(
t
)
<
d
i
a
ATC
0
d
i
a
(
t
)
≥
d
i
a
ATC
.We have discussed the one-dimensional cell transmission model of air traffic flow before, while there may be multiple air routes in terminal area and some of them may cross with each other, in the way, in the same imaginary plane, as shown in Figure3. So the converging and diverging situations are as shown in Figure 3.Figure 3
Converging and diverging of air traffic flow.According to the conservation of number, in converging situation the number of aircraft satisfies(13)
N
k
=
N
k
-
1
+
N
k
-
2
+
⋯
+
N
k
-
n
.On the contrary, in diverging situation we have(14)
N
k
+
1
=
k
k
+
1
N
k
,
N
k
+
2
=
k
k
+
2
N
k
,
⋮
N
k
+
n
=
k
k
+
n
N
k
.The coefficientk stands for the proportion of traffic flow to different air route segments; meanwhile 0
≤
k
n
≤
1, ∑
k
n
=
1.


### 2.3. Departure Traffic Flow Model

The model of departure flow in airport terminal area consists of two phases: taking off from runway and flying along air routes. In the take-off phase, departure air traffic flow will make the most use of runway time slots based on the arrival priority, in order to guarantee a smooth arrival and landing process. Flight departure schedules are generated in stochastic cases, while the departure flights that disagree with the operation time interval of runway will be delayed to ground holding procedures until enough time slots come [20–22]. In the flying stage, departure aircraft diverge in various directions from the runway center, which differs from the arrival flow as there would be both diverging and converging situations in arrival processes. Under the condition of fully isolation between arrival and departure routes, flying stage will dispense with flow control if the runway interval problem is solved in the taking-off stage. And under the condition of semi-isolation between arrival and departure routes, which means that some air routes may be overlapped in several segments (in the same imaginary plane), it is necessary to take both departure aircraft and arrival ones into consideration simultaneously to focus on the average nose interval of all aircraft in the overlapped segments (cells). If the average interval cannot meet safety needs, departure flow should be adjusted prior to the arrivals.LetN
ground
(
t
) be the number of aircraft holding for departure on ground in time t, and its change can be described as follows:
(15)
N
ground
(
t
+
1
)
=
N
ground
(
t
)
+
τ
i
[
q
plan
d
(
t
)
-
q
s
d
(
t
)
]
.In (15), N
ground
(
t
+
1
) is the number of aircraft holding for departure on ground in time t
+
1. The demand of departure per unit time that produced by flight schedules is represented by q
plan
d
(
t
), while q
s
d
(
t
) represents the actual number of aircraft taking off from the runway per unit time. The time step is still τ
i, and N
ground
(
t
) cannot be negative, so we can get
(16)
q
s
d
(
t
)
=
{
N
ground
(
t
)
τ
i
+
q
plan
d
(
t
)
N
ground
(
t
)
+
τ
i
[
q
plan
d
(
t
)
-
C
s
d
(
t
)
]
<
0
C
s
d
(
t
)
N
ground
(
t
)
+
τ
i
[
q
plan
d
(
t
)
-
C
s
d
(
t
)
]
≥
0
.Consider maximizing the use of runway time slots; variableC
s
d
(
t
) represents the maximum take-off rate of the runway, which depends on the time interval of runway operation. If the time interval of arrival landing aircraft gets large, influence such as wake vortex on the runway will have no effect on take-off aircraft, in which cases departure flows may take off by the standard time interval separately. Otherwise if arrival time interval gets intense, mutual interference between departures and arrivals will be strong, so that departure flows need to make use of the interspace of the time slots under the circumstances of an affected runway operation, as follows:
(17)
C
s
d
(
t
)
=
{
d
t
s
a
(
t
)
-
d
T
m
runway
d
T
m
runway
·
d
t
s
a
(
t
)
d
t
s
a
(
t
)
<
d
T
s
runway
1
d
T
s
d
d
t
s
a
(
t
)
≥
d
T
s
runway
.In (17), this d
t
s
a
(
t
) is the actual time interval of arrival landing and d
T
s
d is the standard take-off time interval of departures, while d
T
s
runway and d
T
m
runway, respectively, stand for the irrelevant runway operation time interval and the modified operation time interval according to the mutual interference between landing and taking off. Among this, the irrelevant runway operation time d
T
s
runway will be a critical point where interference occurred.Similar to the Arrival Model, we divide air routes into several segments by a series of unit-interconnected cells in the stage of flying along departure routes. LetN
i
d
(
t
) be the number of departure aircraft in cell i at time t and its changing process can be described as (18). On the left side of the equation N
i
d
(
t
+
1
) is the number of departures in cell i at time t
+
1; q
i
-
1
d
(
t
) represents the flux of departure aircraft entering cell i from cell i
-
1 at time t, while q
i
d
(
t
) means that the flux of departure aircraft exiting cell i at time t, τ
i is the time step.Since most of the time departure aircraft in airport terminal area are in an acceleration process, we assume aircraft entering celli from some certain departure route positions with an initial velocity u
′ that comes from the SID and then uniformly accelerating along the cell (route) with a rate a
i
d. So according to the same modeling principle from the Arrival Model, we can get the following similar equations:
(18)
q
i
d
(
t
)
=
β
i
a
i
d
[
N
i
d
(
t
)
-
N
i
d
ATC
(
t
)
]
u
i
′
2
+
2
a
i
d
Ω
i
-
u
i
′
+
u
i
ATC
N
i
d
ATC
(
t
)
Ω
i
=
β
i
a
i
d
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
N
i
d
(
t
)
-
(
β
i
a
i
d
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
-
u
i
ATC
Ω
i
)
N
i
d
ATC
(
t
)
,
(19)
N
i
d
(
t
+
1
)
=
(
1
-
β
i
a
i
d
τ
i
u
i
′
2
+
2
a
i
d
Ω
i
-
u
i
′
)
N
i
d
(
t
)
+
(
β
i
a
i
d
τ
i
u
i
′
2
+
2
a
i
d
Ω
i
-
u
i
′
-
u
i
ATC
τ
i
Ω
i
)
N
i
d
ATC
(
t
)
+
τ
i
q
i
-
1
d
(
t
)
.In the above equations, coefficientβ
i stands for the rate of departures outflow per unit time. The number of aircraft that need to be arranged to take departure flow control measures is denoted by the variable N
i
d
ATC
(
t
). As mentioned before, in circumstances of fully isolation between arrival and departure route segments there will be no extra control to the departures; that is, N
i
d
ATC
(
t
)
=
0. On the other hand, in the routes overlapped condition part of the departures may take extra controls. We assume that these controlled departure flows would move with a new constant ATC velocity u
i
ATC in cell i and be the same with variable v
i
ATC in the Arrival Model, in which both represent displacement velocity. Meanwhile the simplification form is
(20)
N
i
d
(
t
+
1
)
=
A
i
d
N
i
d
(
t
)
+
B
i
d
τ
i
N
i
d
ATC
(
t
)
+
τ
i
q
i
-
1
d
(
t
)
,
q
i
d
(
t
)
=
C
i
d
N
i
d
(
t
)
-
B
i
d
N
i
d
ATC
(
t
)
,
A
i
d
=
1
-
β
i
a
i
d
τ
i
u
i
′
2
+
2
a
i
d
Ω
i
-
u
i
′
;
B
i
d
=
β
i
a
i
d
u
i
′
2
+
2
a
i
d
Ω
i
-
u
i
′
-
u
i
ATC
Ω
i
;
C
i
d
=
β
i
a
i
d
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
.To determine the value of variableN
i
d
ATC
(
t
) in overlapped segments, we will bring the average nose interval (not time interval) of arrivals noted by d
i
a
(
t
), the average nose interval of departures noted by d
i
d
(
t
), and the average interval between both arrivals and departures noted by d
i
a
/
d
(
t
) in cell i all into consideration. The fundamental aim of flow control in this stage is to make the entire average interval of arrivals and departures meet the ATC separation requirement d
i
ATC by adjusting departure flows, as follows:
(21)
N
i
d
ATC
(
t
)
=
{
Ω
i
d
i
d
(
t
)
d
i
a
(
t
)
≤
d
i
ATC
(
1
d
i
a
D
E
P
(
t
)
-
1
d
i
ATC
)
·
Ω
i
d
i
a
(
t
)
>
d
i
ATC
,
d
i
a
/
d
(
t
)
<
d
i
ATC
0
d
i
a
/
d
(
t
)
≥
d
i
ATC
.The average interval of arrivals and departures isd
i
a
/
d
(
t
)
=
d
i
a
(
t
)
·
d
i
d
(
t
)
/
[
d
i
a
(
t
)
+
d
i
d
(
t
)
].


## 2.1. Definitions

There are many kinds of definitions for air traffic flow parameters, because of different research intentions and different methods that can be used. Since we focus on airport terminal area, the definitions of velocity, density, and traffic flux on a segment of air route are as follows.Traffic flux (q) is the number of aircraft passing a reference profile of the observation segment per unit of time. q
=
N
/
T, where T is the observation time and N represents the number of aircraft passing in T.Density (ρ) is the number of aircraft per unit length of the observation segment. ρ
=
N
/
P, where P means length of the observation segment, while N is the number of aircraft in P.Velocity can be divided into micro and macro definitions. Microscopic definitions including instantaneous velocity (v) and average velocity (v
¯) focus on some point or profile of the observation segment; macroscopic ones including the space mean velocity (v
¯
s) and time mean velocity (v
¯
t) focus on some area extents or time ranges:
(1)
v
=
d
x
d
t
,
v
¯
=
1
N
∑
i
=
1
N
v
i
,
v
¯
s
=
D
·
(
1
N
∑
i
=
1
N
t
i
)
-
1
,
v
¯
t
=
1
N
∑
i
=
1
N
s
i
t
1
-
t
0
.One flight (f
i) of all (N) uses t
i time to pass the observation segment, where the actual distance is D. And s
i is the flying distance of f
i in the time period t
1
-
t
0.Historically, the first macroscopic traffic flow model is a continuity equation, called the Lighthill-Whitham-Richards (LWR) equation [2, 3]. Like the vehicle traffic flow, air traffic flux (q) and linear density (ρ) may also satisfy the corresponding relations for some given functions f
(
·
) and g
(
·
) with location x and time t, as follows:
(2)
ρ
(
x
,
t
)
=
f
(
q
(
x
,
t
)
,
x
)
,
q
(
x
,
t
)
=
g
(
ρ
(
x
,
t
)
,
x
)
,(3)
∂
q
(
x
,
t
)
∂
x
+
∂
ρ
(
x
,
t
)
∂
t
=
s
.On the right side of (3), s denotes the number of aircraft that enter or exit the observation segment. With the basic equation q
=
ρ
v, where v denotes space mean velocity v
¯
s and assuming every flight has the same flying case, we can derive
(4)
v
(
x
)
∂
q
(
x
,
t
)
∂
x
+
∂
q
(
x
,
t
)
∂
t
=
v
(
x
)
s
.There are two variables as density and velocity but only one equation in the LWR model. It cannot be solved as the differential equation is not closed. In response to this problem, a balancing velocity-density functional relationship, asv
(
x
)
=
v
e
(
ρ
(
x
,
t
)
), was introduced into LWR theory by assuming traffic flow invariably in equilibrium state. To plug this into the equation, a hyperbolic equation of density can be derived, as follows:
(5)
∂
q
e
(
ρ
(
x
,
t
)
)
∂
x
+
∂
ρ
(
x
,
t
)
∂
t
=
s
.


## 2.2. Arrival Traffic Flow Model

LWR model described the propagation characteristic of the nonlinear density wave to find the evolution rule of traffic shock wave and rarefaction wave by characteristics method or numerical simulation [19]. However, air traffic flow differs from vehicle traffic, especially in airport terminal area (Figure 1). First, the arrival and departure traffic in terminal area will follow the designed STAR/SID (standard terminal arrival route/standard instrument departure), which means that aircraft in every route position must act in accordance with the operational flight program, that is, flying within a certain scope of designed flight level and speed. It is usually a small scope and can be reassigned by controllers. As a result of that, some balancing velocity-density functional relationships in vehicle traffic may not exist in air traffic flow. Second, the density of air traffic flow is generally much lower than the vehicle traffic in real operations, which causes the interactions among aircraft to be weaker compared to vehicles. Therefore, it will be difficult to make statistical fit of the relationships according to the available radar data.Figure 1
Radar track plot in airport terminal area (source: ZGGG TMA).To solve the partial differential equation of the LWR theory, we use cell transmission model to discretize the continuity equation of macroscopic traffic flow. The method of discretization is applying a series of interconnected one-dimensional cells to denote the air route and using the difference equation of time discretization to describe aircraft passing through every adjoining cell. It should be feasible to model any air traffic flow scenarios by using interconnected cells, as shown in Figure2, putting all the aircraft flying along the arrival route while in several flight levels onto the same imaginary plane, hence the air traffic flow in one route can be approximately seen as one-dimensional continuous flow. In addition, one air route will be divided into several segments by series of unit cells. To simplify the model we assume that there is a unified variable which we will introduce later to represent different kinds of control measures in one cell, such as speed control, maneuvering actions, or circling, that is, adjusting flows by changing speed or flight path of certain aircraft in the cell.Figure 2
One-dimensional cell transmission model of air traffic flow in single direction.In the matter of arrival route, letN
i
a
(
t
) be the number of arrival aircraft in cell i at time t; then the changes in the number of arrival aircraft in one cell can be described by the following difference equation of time discretization:
(6)
N
i
a
(
t
+
1
)
=
N
i
a
(
t
)
+
τ
i
[
q
i
-
1
a
(
t
)
-
q
i
a
(
t
)
]
.In the equation above,N
i
a
(
t
+
1
) is the number of arrival aircraft in cell i at time t
+
1, q
i
-
1
a
(
t
) represents the flux of arrival aircraft entering cell i from cell i
-
1 at time t, while q
i
a
(
t
) means the flux of arrival aircraft exiting cell i at time t, and τ
i is the time step.Note. The number, flux et al. mentioned in this section, is specific to arrival aircraft, while departures will be included in following parts.Air traffic flow in air route segment still satisfies the basic equationq
i
=
ρ
i
v
i, where ρ
i and v
i stand, respectively, for segment liner density and space mean velocity of cell i. The traffic flux of cell i will be q
i
(
t
)
=
α
i
q
i. Coefficient α
i is the rate of outflow per unit time, which reflects the saturation level of cell i. From Section 2.1., we know segment liner density ρ
i
=
N
i
/
Ω
i, where Ω
i is the length of cell i. Since most of the time arrival aircraft in airport terminal area are in a deceleration process, we assume aircraft entering cell i from some certain arrival route positions with an initial velocity v
i
′ that comes from the STAR and then uniformly decelerating along the cell (route) with a rate a
i
a, so we can get the space mean velocity of cell i, as follows:
(7)
v
i
=
a
i
a
Ω
i
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
.To plug (7) into q
i
=
ρ
i
v
i, we can derive
(8)
q
i
a
(
t
)
=
α
i
a
i
a
N
i
a
(
t
)
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
.However, since the aforementioned situation is an ideal condition, it is necessary to take control measures for part of the arrival aircraft into consideration, such as speed control, maneuvering actions, or circling, because of the existence of traffic congestion, safety interval, and so forth in real conditions. For simplicity, we assume that there is a unified variable to stand for different kinds of control measures, since whichever measures controllers took, the actual effects caused by changing speed or flight path and so forth can all be seen as there to beN
i
a
ATC aircraft flying along the cell (route) with a new constant velocity v
i
ATC in cell i. To be clear, the variable v
i
ATC here represents the displacement velocity that is lower than the STAR designed velocity generally.Considering these problems, we derived (9) after integrating control measures into (8), as follows:
(9)
q
i
a
(
t
)
=
α
i
a
i
a
[
N
i
a
(
t
)
-
N
i
a
ATC
(
t
)
]
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
+
v
i
ATC
N
i
a
ATC
(
t
)
Ω
i
=
α
i
a
i
a
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
N
i
a
(
t
)
-
(
α
i
a
i
a
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
-
v
i
ATC
Ω
i
)
N
i
a
ATC
(
t
)
.After substitutingq
i
a
(
t
) in the difference equation (6) with (9), we got
(10)
N
i
a
(
t
+
1
)
=
(
1
-
α
i
a
i
a
τ
i
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
N
i
a
(
t
)
+
(
α
i
a
i
a
τ
i
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
-
v
i
ATC
τ
i
Ω
i
)
N
i
a
ATC
(
t
)
+
τ
i
q
i
-
1
a
(
t
)
.LetA
i
a
=
1
-
(
α
i
a
i
a
τ
i
/
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
); B
i
a
=
(
α
i
a
i
a
/
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
)
-
(
v
i
ATC
/
Ω
i
); C
i
a
=
α
i
a
i
a
/
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
); we can have the simplified equation sets, as follows:
(11)
N
i
a
(
t
+
1
)
=
A
i
a
N
i
a
(
t
)
+
B
i
a
τ
i
N
i
a
ATC
(
t
)
+
τ
i
q
i
-
1
a
(
t
)
,
q
i
a
(
t
)
=
C
i
a
N
i
a
(
t
)
-
B
i
a
N
i
a
ATC
(
t
)
.The aircraft in cell are taken as evenly distributed in the macroscopic traffic flow theory, so the number of aircraft should be equivalent number. Letd
i
a
(
t
) be the mean nose interval of adjacent aircraft at the same direction in cell i; therefore the number of aircraft in the cell is N
i
(
t
)
=
Ω
i
/
d
i
a
(
t
). If d
i
a
(
t
) becomes less than the control separation standard denoted by d
i
a
ATC, that is to say, aircraft density exceeds the threshold level, the exceeded aircraft should be arranged to take control measures. In addition, arrival aircraft in airport terminal area also have to obey the safety separation d
i
safe according to STAR; that is, d
i
a
(
t
)
≥
d
i
safe. The safety separation is evaluated by air traffic security department and usually gets more stringent compared to control separation as d
i
safe
<
d
i
a
ATC. The number of controlled aircraft at time t in cell i satisfies
(12)
N
i
a
ATC
(
t
)
=
{
(
1
d
i
a
(
t
)
-
1
d
i
a
ATC
)
Ω
i
d
i
a
(
t
)
<
d
i
a
ATC
0
d
i
a
(
t
)
≥
d
i
a
ATC
.We have discussed the one-dimensional cell transmission model of air traffic flow before, while there may be multiple air routes in terminal area and some of them may cross with each other, in the way, in the same imaginary plane, as shown in Figure3. So the converging and diverging situations are as shown in Figure 3.Figure 3
Converging and diverging of air traffic flow.According to the conservation of number, in converging situation the number of aircraft satisfies(13)
N
k
=
N
k
-
1
+
N
k
-
2
+
⋯
+
N
k
-
n
.On the contrary, in diverging situation we have(14)
N
k
+
1
=
k
k
+
1
N
k
,
N
k
+
2
=
k
k
+
2
N
k
,
⋮
N
k
+
n
=
k
k
+
n
N
k
.The coefficientk stands for the proportion of traffic flow to different air route segments; meanwhile 0
≤
k
n
≤
1, ∑
k
n
=
1.


## 2.3. Departure Traffic Flow Model

The model of departure flow in airport terminal area consists of two phases: taking off from runway and flying along air routes. In the take-off phase, departure air traffic flow will make the most use of runway time slots based on the arrival priority, in order to guarantee a smooth arrival and landing process. Flight departure schedules are generated in stochastic cases, while the departure flights that disagree with the operation time interval of runway will be delayed to ground holding procedures until enough time slots come [20–22]. In the flying stage, departure aircraft diverge in various directions from the runway center, which differs from the arrival flow as there would be both diverging and converging situations in arrival processes. Under the condition of fully isolation between arrival and departure routes, flying stage will dispense with flow control if the runway interval problem is solved in the taking-off stage. And under the condition of semi-isolation between arrival and departure routes, which means that some air routes may be overlapped in several segments (in the same imaginary plane), it is necessary to take both departure aircraft and arrival ones into consideration simultaneously to focus on the average nose interval of all aircraft in the overlapped segments (cells). If the average interval cannot meet safety needs, departure flow should be adjusted prior to the arrivals.LetN
ground
(
t
) be the number of aircraft holding for departure on ground in time t, and its change can be described as follows:
(15)
N
ground
(
t
+
1
)
=
N
ground
(
t
)
+
τ
i
[
q
plan
d
(
t
)
-
q
s
d
(
t
)
]
.In (15), N
ground
(
t
+
1
) is the number of aircraft holding for departure on ground in time t
+
1. The demand of departure per unit time that produced by flight schedules is represented by q
plan
d
(
t
), while q
s
d
(
t
) represents the actual number of aircraft taking off from the runway per unit time. The time step is still τ
i, and N
ground
(
t
) cannot be negative, so we can get
(16)
q
s
d
(
t
)
=
{
N
ground
(
t
)
τ
i
+
q
plan
d
(
t
)
N
ground
(
t
)
+
τ
i
[
q
plan
d
(
t
)
-
C
s
d
(
t
)
]
<
0
C
s
d
(
t
)
N
ground
(
t
)
+
τ
i
[
q
plan
d
(
t
)
-
C
s
d
(
t
)
]
≥
0
.Consider maximizing the use of runway time slots; variableC
s
d
(
t
) represents the maximum take-off rate of the runway, which depends on the time interval of runway operation. If the time interval of arrival landing aircraft gets large, influence such as wake vortex on the runway will have no effect on take-off aircraft, in which cases departure flows may take off by the standard time interval separately. Otherwise if arrival time interval gets intense, mutual interference between departures and arrivals will be strong, so that departure flows need to make use of the interspace of the time slots under the circumstances of an affected runway operation, as follows:
(17)
C
s
d
(
t
)
=
{
d
t
s
a
(
t
)
-
d
T
m
runway
d
T
m
runway
·
d
t
s
a
(
t
)
d
t
s
a
(
t
)
<
d
T
s
runway
1
d
T
s
d
d
t
s
a
(
t
)
≥
d
T
s
runway
.In (17), this d
t
s
a
(
t
) is the actual time interval of arrival landing and d
T
s
d is the standard take-off time interval of departures, while d
T
s
runway and d
T
m
runway, respectively, stand for the irrelevant runway operation time interval and the modified operation time interval according to the mutual interference between landing and taking off. Among this, the irrelevant runway operation time d
T
s
runway will be a critical point where interference occurred.Similar to the Arrival Model, we divide air routes into several segments by a series of unit-interconnected cells in the stage of flying along departure routes. LetN
i
d
(
t
) be the number of departure aircraft in cell i at time t and its changing process can be described as (18). On the left side of the equation N
i
d
(
t
+
1
) is the number of departures in cell i at time t
+
1; q
i
-
1
d
(
t
) represents the flux of departure aircraft entering cell i from cell i
-
1 at time t, while q
i
d
(
t
) means that the flux of departure aircraft exiting cell i at time t, τ
i is the time step.Since most of the time departure aircraft in airport terminal area are in an acceleration process, we assume aircraft entering celli from some certain departure route positions with an initial velocity u
′ that comes from the SID and then uniformly accelerating along the cell (route) with a rate a
i
d. So according to the same modeling principle from the Arrival Model, we can get the following similar equations:
(18)
q
i
d
(
t
)
=
β
i
a
i
d
[
N
i
d
(
t
)
-
N
i
d
ATC
(
t
)
]
u
i
′
2
+
2
a
i
d
Ω
i
-
u
i
′
+
u
i
ATC
N
i
d
ATC
(
t
)
Ω
i
=
β
i
a
i
d
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
N
i
d
(
t
)
-
(
β
i
a
i
d
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
-
u
i
ATC
Ω
i
)
N
i
d
ATC
(
t
)
,
(19)
N
i
d
(
t
+
1
)
=
(
1
-
β
i
a
i
d
τ
i
u
i
′
2
+
2
a
i
d
Ω
i
-
u
i
′
)
N
i
d
(
t
)
+
(
β
i
a
i
d
τ
i
u
i
′
2
+
2
a
i
d
Ω
i
-
u
i
′
-
u
i
ATC
τ
i
Ω
i
)
N
i
d
ATC
(
t
)
+
τ
i
q
i
-
1
d
(
t
)
.In the above equations, coefficientβ
i stands for the rate of departures outflow per unit time. The number of aircraft that need to be arranged to take departure flow control measures is denoted by the variable N
i
d
ATC
(
t
). As mentioned before, in circumstances of fully isolation between arrival and departure route segments there will be no extra control to the departures; that is, N
i
d
ATC
(
t
)
=
0. On the other hand, in the routes overlapped condition part of the departures may take extra controls. We assume that these controlled departure flows would move with a new constant ATC velocity u
i
ATC in cell i and be the same with variable v
i
ATC in the Arrival Model, in which both represent displacement velocity. Meanwhile the simplification form is
(20)
N
i
d
(
t
+
1
)
=
A
i
d
N
i
d
(
t
)
+
B
i
d
τ
i
N
i
d
ATC
(
t
)
+
τ
i
q
i
-
1
d
(
t
)
,
q
i
d
(
t
)
=
C
i
d
N
i
d
(
t
)
-
B
i
d
N
i
d
ATC
(
t
)
,
A
i
d
=
1
-
β
i
a
i
d
τ
i
u
i
′
2
+
2
a
i
d
Ω
i
-
u
i
′
;
B
i
d
=
β
i
a
i
d
u
i
′
2
+
2
a
i
d
Ω
i
-
u
i
′
-
u
i
ATC
Ω
i
;
C
i
d
=
β
i
a
i
d
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
.To determine the value of variableN
i
d
ATC
(
t
) in overlapped segments, we will bring the average nose interval (not time interval) of arrivals noted by d
i
a
(
t
), the average nose interval of departures noted by d
i
d
(
t
), and the average interval between both arrivals and departures noted by d
i
a
/
d
(
t
) in cell i all into consideration. The fundamental aim of flow control in this stage is to make the entire average interval of arrivals and departures meet the ATC separation requirement d
i
ATC by adjusting departure flows, as follows:
(21)
N
i
d
ATC
(
t
)
=
{
Ω
i
d
i
d
(
t
)
d
i
a
(
t
)
≤
d
i
ATC
(
1
d
i
a
D
E
P
(
t
)
-
1
d
i
ATC
)
·
Ω
i
d
i
a
(
t
)
>
d
i
ATC
,
d
i
a
/
d
(
t
)
<
d
i
ATC
0
d
i
a
/
d
(
t
)
≥
d
i
ATC
.The average interval of arrivals and departures isd
i
a
/
d
(
t
)
=
d
i
a
(
t
)
·
d
i
d
(
t
)
/
[
d
i
a
(
t
)
+
d
i
d
(
t
)
].


## 3. Parameter Analysis

We focus on one single cell in the model as to deduce and analyze the interrelationship of air traffic flow characteristic parameters including flight flux, linear density, and traffic velocity. In the following sections, we discuss this problem from the two different aspects as arrival and departure, just like the models we had established before.


### 3.1. Arrival Routes

To plug (12) into (9), we can have
(22)
q
i
a
(
t
)
=
{
α
i
a
i
a
Ω
i
d
i
a
ATC
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
+
v
i
ATC
(
1
d
i
a
(
t
)
-
1
d
i
a
ATC
)
d
i
a
(
t
)
<
d
i
a
ATC
α
i
a
i
a
Ω
i
d
i
a
(
t
)
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
d
i
a
(
t
)
≥
d
i
a
ATC
.According toρ
i
a
(
t
)
=
1
/
d
i
a
(
t
), we get
(23)
q
i
a
(
t
)
=
{
v
i
ATC
ρ
i
a
(
t
)
+
α
i
a
i
a
Ω
i
+
v
i
ATC
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
d
i
a
ATC
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
ρ
i
a
(
t
)
>
1
d
i
a
ATC
α
i
a
i
a
Ω
i
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
ρ
i
a
(
t
)
ρ
i
a
(
t
)
≤
1
d
i
a
ATC
.From (23), we can see that the basic parameters such as flux q
i
a
(
t
) and density ρ
i
a
(
t
) form a piecewise function and the segment point is the reciprocal value of ATC separation d
i
a
ATC to the arrival aircraft. Assuming that coefficient α
i, standard deceleration a
i
a, initial speed v
i
′, and length Ω
i of one cell are all constant values, the slope and intercept of the relation curve will be determined by ATC velocity v
i
ATC and the position of inflection point will be determined by d
i
a
ATC. Let both sides of (23) be divided by ρ
i
a
(
t
):
(24)
v
i
(
t
)
=
{
α
i
a
i
a
Ω
i
+
v
i
ATC
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
d
i
a
ATC
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
·
1
ρ
i
a
(
t
)
+
v
i
ATC
ρ
i
a
(
t
)
>
1
d
i
a
ATC
α
i
a
i
a
Ω
i
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
ρ
i
a
(
t
)
≤
1
d
i
a
ATC
.If aircraft density comes lower than the critical value, the mean velocity of traffic flow in one cell will stay constant; otherwise, there will be an inverse relation between velocityv
i
(
t
) and density ρ
i
a
(
t
). With an increase of density in cell i, mean velocity will become lower gradually and the speed of reducing just gets more and more slow and eventually tends to the ATC velocity v
i
ATC. Using the basic formula of fluid q
=
ρ
v, we can get (25)
q
i
a
(
t
)
=
{
α
i
a
i
a
Ω
i
+
v
i
ATC
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
d
i
a
ATC
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
·
v
i
(
t
)
v
i
(
t
)
-
v
i
ATC
v
i
ATC
<
v
i
(
t
)
<
α
i
a
i
a
Ω
i
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
∃
q
i
a
(
t
)
∈
[
0
,
α
i
a
i
a
Ω
i
d
i
a
ATC
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
]
v
i
(
t
)
=
v
i
ATC
,
α
i
a
i
a
Ω
i
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
.From (25), we can see that the mean velocity of entire traffic flow in cell i ranges from ATC velocity v
i
ATC to standard velocity designed by STAR. If v
i
(
t
) is equal to one of the critical values flux q
i
a
(
t
) may be arbitrary-sized data taking from zero to the maximum flux value. With an increase of mean velocity in cell i, the flux value goes down gradually and the speed of reducing gets slow. When this mean velocity goes up to the peak value, it will cause a jump of traffic flux. Variables d
i
a
ATC and v
i
ATC are still the order parameters of the interrelationship between flux and velocity in cell.


### 3.2. Departure Routes

In the stage of flying along departure routes, aircraft keep accelerating and diverging along with the varying of orientation of air routes. There may beN
i
d
ATC
(
t
) aircraft taken departure flow control measures in some of the overlapped segments. As an exceptional case of this situation, flying stage will dispense with flow control in the case of fully isolation between arrival and departure routes; that is, N
i
d
ATC
(
t
)
=
0. We will focus on the overlapped segments and analyze the traffic flow characteristic parameters of departures in detail. To plug (21) into (18) we can have
(26)
q
i
d
(
t
)
=
{
u
i
ATC
d
i
d
(
t
)
d
i
a
(
t
)
≤
d
i
ATC
u
i
ATC
d
i
d
(
t
)
+
d
i
a
(
t
)
-
d
i
ATC
d
i
ATC
d
i
a
(
t
)
×
(
β
i
a
i
d
Ω
i
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
-
u
i
ATC
)
d
i
a
(
t
)
>
d
i
ATC
,
d
i
a
/
d
(
t
)
<
d
i
ATC
β
i
a
i
d
d
i
d
(
t
)
(
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
)
d
i
t
(
t
)
≥
d
i
ATC
.Sinceρ
i
d
(
t
)
=
1
/
d
i
d
(
t
), ρ
i
a
(
t
)
=
1
/
d
i
a
(
t
), substituting in (26), we get
(27)
q
i
d
(
t
)
=
{
β
i
a
i
d
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
ρ
i
d
(
t
)
ρ
i
a
(
t
)
+
ρ
i
d
(
t
)
<
1
d
i
ATC
u
i
ATC
ρ
i
d
(
t
)
+
1
-
d
i
ATC
ρ
i
a
(
t
)
d
i
ATC
×
(
β
i
a
i
d
Ω
i
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
-
u
i
ATC
)
ρ
i
a
(
t
)
<
1
d
i
ATC
,
ρ
i
a
(
t
)
+
ρ
i
d
(
t
)
>
1
d
i
ATC
u
i
ATC
ρ
i
d
(
t
)
ρ
i
a
(
t
)
≥
1
d
i
ATC
.From (27), we can see that flux q
i
d
(
t
) and density ρ
i
d
(
t
) form a piecewise function. Similar to the arrival segments (cells), we assume that coefficient β
i, standard acceleration a
i
d, initial speed u
′, and length Ω
i of one cell are all constant values. The density of arrival aircraft denoted by ρ
i
a
(
t
) in cell i changes over time, which has an effect on the intercept value of this linear equation. The ATC velocity of departures denoted by u
i
ATC shows the slope while the reciprocal of ATC separation d
i
ATC is the critical value of flow density, determining the position of inflection points. Then, letting both sides of (27) be divided by ρ
i
d
(
t
), we got
(28)
u
i
(
t
)
=
{
β
i
a
i
d
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
ρ
i
a
(
t
)
+
ρ
i
d
(
t
)
<
1
d
i
ATC
u
i
ATC
+
1
-
d
i
ATC
ρ
i
a
(
t
)
d
i
ATC
ρ
i
d
(
t
)
×
(
β
i
a
i
d
Ω
i
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
-
u
i
ATC
)
ρ
i
a
(
t
)
<
1
d
i
ATC
ρ
i
a
(
t
)
+
ρ
i
d
(
t
)
>
1
d
i
ATC
u
i
ATC
ρ
i
a
(
t
)
≥
1
d
i
ATC
.From (28), we can see that mean velocity u
i
(
t
) and density ρ
i
d
(
t
) also form a piecewise function. Two extreme cases are as follows: the entire density of arrivals and departures is lower than the ATC density that is denoted by 1
/
d
i
ATC; density only considers arrivals that have already exceeded 1
/
d
i
ATC. These two cases will cause mean velocity of departures to be the standard velocity determined by β
i, acceleration a
i
d, initial speed u
′, length Ω
i, and departure ATC velocity u
i
ATC, respectively. In a certain condition between such two extreme cases the slope of linear equation is greater than or equal to u
i
ATC, while the intercept value is influenced by both u
i
ATC and ρ
i
a
(
t
). According to the basic formula of fluid q
=
ρ
v, we get(29)
q
i
d
(
t
)
=
{
(
1
-
d
i
ATC
ρ
i
a
(
t
)
)
[
β
i
a
i
d
Ω
i
-
u
i
ATC
(
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
)
]
d
i
ATC
(
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
)
·
u
i
(
t
)
u
i
(
t
)
-
u
i
ATC
u
i
ATC
<
u
i
(
t
)
u
i
ATC
<
β
i
a
i
d
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
∃
q
i
d
(
t
)
∈
[
0
,
β
i
a
i
d
Ω
i
d
i
d
ATC
(
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
)
]
u
i
(
t
)
=
u
i
ATC
,
β
i
a
i
d
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
.The form of (29) is much like (25) from the arrival parts. If mean velocity of departures lies between the standard velocity and the departure ATC velocity u
i
ATC there will be an inverse proportional function between departure flux q
i
d
(
t
) and mean velocity u
i
(
t
). Moreover, the coefficient of this inverse proportional function is determined by β
i, acceleration a
i
d, initial speed u
′, and cell length Ω
i all together and changes in value of the ATC velocity u
i
ATC will make the function curve move horizontally. While if mean velocity of departures lies on either of the two critical positions departure flux may vary randomly.


## 3.1. Arrival Routes

To plug (12) into (9), we can have
(22)
q
i
a
(
t
)
=
{
α
i
a
i
a
Ω
i
d
i
a
ATC
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
+
v
i
ATC
(
1
d
i
a
(
t
)
-
1
d
i
a
ATC
)
d
i
a
(
t
)
<
d
i
a
ATC
α
i
a
i
a
Ω
i
d
i
a
(
t
)
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
d
i
a
(
t
)
≥
d
i
a
ATC
.According toρ
i
a
(
t
)
=
1
/
d
i
a
(
t
), we get
(23)
q
i
a
(
t
)
=
{
v
i
ATC
ρ
i
a
(
t
)
+
α
i
a
i
a
Ω
i
+
v
i
ATC
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
d
i
a
ATC
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
ρ
i
a
(
t
)
>
1
d
i
a
ATC
α
i
a
i
a
Ω
i
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
ρ
i
a
(
t
)
ρ
i
a
(
t
)
≤
1
d
i
a
ATC
.From (23), we can see that the basic parameters such as flux q
i
a
(
t
) and density ρ
i
a
(
t
) form a piecewise function and the segment point is the reciprocal value of ATC separation d
i
a
ATC to the arrival aircraft. Assuming that coefficient α
i, standard deceleration a
i
a, initial speed v
i
′, and length Ω
i of one cell are all constant values, the slope and intercept of the relation curve will be determined by ATC velocity v
i
ATC and the position of inflection point will be determined by d
i
a
ATC. Let both sides of (23) be divided by ρ
i
a
(
t
):
(24)
v
i
(
t
)
=
{
α
i
a
i
a
Ω
i
+
v
i
ATC
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
d
i
a
ATC
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
·
1
ρ
i
a
(
t
)
+
v
i
ATC
ρ
i
a
(
t
)
>
1
d
i
a
ATC
α
i
a
i
a
Ω
i
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
ρ
i
a
(
t
)
≤
1
d
i
a
ATC
.If aircraft density comes lower than the critical value, the mean velocity of traffic flow in one cell will stay constant; otherwise, there will be an inverse relation between velocityv
i
(
t
) and density ρ
i
a
(
t
). With an increase of density in cell i, mean velocity will become lower gradually and the speed of reducing just gets more and more slow and eventually tends to the ATC velocity v
i
ATC. Using the basic formula of fluid q
=
ρ
v, we can get (25)
q
i
a
(
t
)
=
{
α
i
a
i
a
Ω
i
+
v
i
ATC
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
d
i
a
ATC
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
·
v
i
(
t
)
v
i
(
t
)
-
v
i
ATC
v
i
ATC
<
v
i
(
t
)
<
α
i
a
i
a
Ω
i
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
∃
q
i
a
(
t
)
∈
[
0
,
α
i
a
i
a
Ω
i
d
i
a
ATC
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
]
v
i
(
t
)
=
v
i
ATC
,
α
i
a
i
a
Ω
i
(
v
i
′
-
v
i
′
2
-
2
a
i
a
Ω
i
)
.From (25), we can see that the mean velocity of entire traffic flow in cell i ranges from ATC velocity v
i
ATC to standard velocity designed by STAR. If v
i
(
t
) is equal to one of the critical values flux q
i
a
(
t
) may be arbitrary-sized data taking from zero to the maximum flux value. With an increase of mean velocity in cell i, the flux value goes down gradually and the speed of reducing gets slow. When this mean velocity goes up to the peak value, it will cause a jump of traffic flux. Variables d
i
a
ATC and v
i
ATC are still the order parameters of the interrelationship between flux and velocity in cell.


## 3.2. Departure Routes

In the stage of flying along departure routes, aircraft keep accelerating and diverging along with the varying of orientation of air routes. There may beN
i
d
ATC
(
t
) aircraft taken departure flow control measures in some of the overlapped segments. As an exceptional case of this situation, flying stage will dispense with flow control in the case of fully isolation between arrival and departure routes; that is, N
i
d
ATC
(
t
)
=
0. We will focus on the overlapped segments and analyze the traffic flow characteristic parameters of departures in detail. To plug (21) into (18) we can have
(26)
q
i
d
(
t
)
=
{
u
i
ATC
d
i
d
(
t
)
d
i
a
(
t
)
≤
d
i
ATC
u
i
ATC
d
i
d
(
t
)
+
d
i
a
(
t
)
-
d
i
ATC
d
i
ATC
d
i
a
(
t
)
×
(
β
i
a
i
d
Ω
i
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
-
u
i
ATC
)
d
i
a
(
t
)
>
d
i
ATC
,
d
i
a
/
d
(
t
)
<
d
i
ATC
β
i
a
i
d
d
i
d
(
t
)
(
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
)
d
i
t
(
t
)
≥
d
i
ATC
.Sinceρ
i
d
(
t
)
=
1
/
d
i
d
(
t
), ρ
i
a
(
t
)
=
1
/
d
i
a
(
t
), substituting in (26), we get
(27)
q
i
d
(
t
)
=
{
β
i
a
i
d
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
ρ
i
d
(
t
)
ρ
i
a
(
t
)
+
ρ
i
d
(
t
)
<
1
d
i
ATC
u
i
ATC
ρ
i
d
(
t
)
+
1
-
d
i
ATC
ρ
i
a
(
t
)
d
i
ATC
×
(
β
i
a
i
d
Ω
i
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
-
u
i
ATC
)
ρ
i
a
(
t
)
<
1
d
i
ATC
,
ρ
i
a
(
t
)
+
ρ
i
d
(
t
)
>
1
d
i
ATC
u
i
ATC
ρ
i
d
(
t
)
ρ
i
a
(
t
)
≥
1
d
i
ATC
.From (27), we can see that flux q
i
d
(
t
) and density ρ
i
d
(
t
) form a piecewise function. Similar to the arrival segments (cells), we assume that coefficient β
i, standard acceleration a
i
d, initial speed u
′, and length Ω
i of one cell are all constant values. The density of arrival aircraft denoted by ρ
i
a
(
t
) in cell i changes over time, which has an effect on the intercept value of this linear equation. The ATC velocity of departures denoted by u
i
ATC shows the slope while the reciprocal of ATC separation d
i
ATC is the critical value of flow density, determining the position of inflection points. Then, letting both sides of (27) be divided by ρ
i
d
(
t
), we got
(28)
u
i
(
t
)
=
{
β
i
a
i
d
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
ρ
i
a
(
t
)
+
ρ
i
d
(
t
)
<
1
d
i
ATC
u
i
ATC
+
1
-
d
i
ATC
ρ
i
a
(
t
)
d
i
ATC
ρ
i
d
(
t
)
×
(
β
i
a
i
d
Ω
i
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
-
u
i
ATC
)
ρ
i
a
(
t
)
<
1
d
i
ATC
ρ
i
a
(
t
)
+
ρ
i
d
(
t
)
>
1
d
i
ATC
u
i
ATC
ρ
i
a
(
t
)
≥
1
d
i
ATC
.From (28), we can see that mean velocity u
i
(
t
) and density ρ
i
d
(
t
) also form a piecewise function. Two extreme cases are as follows: the entire density of arrivals and departures is lower than the ATC density that is denoted by 1
/
d
i
ATC; density only considers arrivals that have already exceeded 1
/
d
i
ATC. These two cases will cause mean velocity of departures to be the standard velocity determined by β
i, acceleration a
i
d, initial speed u
′, length Ω
i, and departure ATC velocity u
i
ATC, respectively. In a certain condition between such two extreme cases the slope of linear equation is greater than or equal to u
i
ATC, while the intercept value is influenced by both u
i
ATC and ρ
i
a
(
t
). According to the basic formula of fluid q
=
ρ
v, we get(29)
q
i
d
(
t
)
=
{
(
1
-
d
i
ATC
ρ
i
a
(
t
)
)
[
β
i
a
i
d
Ω
i
-
u
i
ATC
(
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
)
]
d
i
ATC
(
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
)
·
u
i
(
t
)
u
i
(
t
)
-
u
i
ATC
u
i
ATC
<
u
i
(
t
)
u
i
ATC
<
β
i
a
i
d
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
∃
q
i
d
(
t
)
∈
[
0
,
β
i
a
i
d
Ω
i
d
i
d
ATC
(
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
)
]
u
i
(
t
)
=
u
i
ATC
,
β
i
a
i
d
u
i
′
2
+
2
a
i
d
Ω
i
-
u
′
.The form of (29) is much like (25) from the arrival parts. If mean velocity of departures lies between the standard velocity and the departure ATC velocity u
i
ATC there will be an inverse proportional function between departure flux q
i
d
(
t
) and mean velocity u
i
(
t
). Moreover, the coefficient of this inverse proportional function is determined by β
i, acceleration a
i
d, initial speed u
′, and cell length Ω
i all together and changes in value of the ATC velocity u
i
ATC will make the function curve move horizontally. While if mean velocity of departures lies on either of the two critical positions departure flux may vary randomly.


## 4. Simulation Experiment

### 4.1. Simulation Sample

Based on the simulation platform NetLogo [23–25] we focused on each cell in the macroscopic model and designed traffic inflow/outflow behaviors to be each agent, so as to take control of the inflow/outflow volumes among all the cells. Applying the STAR/SID procedures to runway 02L of Guangzhou Baiyun International Airport (as shown in Figure 4) into the Netlogo system dynamic simulator (as shown in Figure 5) we can emulate the operation and evolution process of air traffic flow in the airport terminal area.Figure 4
STAR/SID routes to RWY02L ZGGG.Figure 5
NetLogo system dynamic simulator.The network of arrival and departure routes in terminal can be separated into several single direction segments and converging or diverging segments that can be further decomposed into more unit cells. For the purpose of simplicity, we set one cell dimension to be equal to the length of the shortest segment in the network and the other air segments to be integer multiples of unit cell dimension. According to this, we made the following Tables1 and 2 of the terminal network.Table 1
Cell quantities of arrival routes to RWY02L ZGGG.
Segment
Code
Quantity
MUBEL-GYA
Segment1
3 cells
BIPOP-GYA
Segment2
2 cells
OSIKA-GYA
Segment3
3 cells
P101-GYA
Segment4
2 cells
GYA-AGVOS
Segment5
2 cells
GYA-D5.1POU
Segment6
2 cells
LONGTANG-TAN
Segment7
2 cells
DANZHU-TAN
Segment8
2 cells
TAN-AGVOS
Segment9
2 cells
AGVOS-D15.7IOO
Segment10
1 cell
AGVOS-D5.1POU
Segment11
1 cell
D5.1POU-D15.7IOO
Segment12
1 cell
D15.7IOO-RWY02L
Segment13
1 cellTable 2
Cell quantities of departure routes to RWY02L ZGGG.
Segment
Code
Quantity
RWY02L-D12.0IOO
Segment14
1 cell
D12.0IOO-YIN
Segment15
2 cells
D12.0IOO-TAN
Segment16
1 cell
TAN-D18.2TAN
Segment17
2 cells
D18.2TAN-POU
Segment18
3 cells
TAN-AGVOS
Segment19
2 cells
AGVOS-D5.1POU
Segment20
1 cell
D5.1POU-POU
Segment21
2 cellsEach yellow rectangle in Figure4 represents one cell, which means a “stock” of fluid in the simulator system. Each grey pipe represents flowing between adjacent cells, which means “flow” in system, and the direction of arrow is the same as the flow direction. Each black valve represents control measures, which means “strategy” in system, and the inflow and outflow of air traffic will be under control by giving the certain valves some corresponding rules. Take typical overlapped segments “TAN=AGVOS (Segment9&19)” as an example, the cell transmission model of this is like (30). In the model, N
i
/
j
a
A
i
/
j
a
B
i
/
j
a
C
i
/
j
a
, respectively, represents the number of arrival aircraft and the simplified coefficients of cell j in Segment i. N
i
/
j
a
ATC denotes the number of arrival aircraft taking control measures in cell j of Segment i. q
i
/
j
a
(
t
) means the number of arrival aircraft that exit cell j in Segment i per unit time. N
i
/
j
d
A
i
/
j
d
B
i
/
j
d
C
i
/
j
d
N
i
/
j
d
ATC
q
i
/
j
d
(
t
)
, respectively, represents the similar corresponding ones but specially for departures. Coefficient k stands for the proportion of traffic flow into some segment when diverging occurs in air routes. τ is the designed time step for the simulator system:
(30)
N
9
/
1
a
(
t
+
1
)
=
A
9
/
1
a
N
9
/
1
a
(
t
)
+
B
9
/
1
a
τ
N
9
/
1
a
ATC
(
t
)
+
τ
[
q
7
/
2
a
(
t
)
+
q
8
/
2
a
(
t
)
]
,
q
9
/
1
a
(
t
)
=
C
9
/
1
a
N
9
/
1
a
(
t
)
-
B
9
/
1
a
N
9
/
1
a
ATC
(
t
)
,
N
9
/
2
a
(
t
+
1
)
=
A
9
/
2
a
N
9
/
2
a
(
t
)
+
B
9
/
2
a
τ
N
9
/
2
a
ATC
(
t
)
+
τ
q
9
/
1
a
(
t
)
,
q
9
/
2
-
10
/
1
a
(
t
)
=
k
9
/
2
-
10
/
1
[
C
9
/
2
a
N
9
/
2
a
(
t
)
-
B
9
/
2
N
9
/
2
a
ATC
(
t
)
]
,
q
9
/
2
-
11
/
1
a
(
t
)
=
k
9
/
2
-
11
/
1
[
C
9
/
2
a
N
9
/
2
a
(
t
)
-
B
9
/
2
N
9
/
2
a
ATC
(
t
)
]
,
(31)
N
19
/
1
d
(
t
+
1
)
=
A
19
/
1
d
N
19
/
1
d
(
t
)
+
B
19
/
1
d
τ
N
19
/
1
d
ATC
(
t
)
+
τ
q
16
/
1
d
(
t
)
,
q
19
/
1
d
(
t
)
=
C
19
/
1
d
N
19
/
1
d
(
t
)
-
B
19
/
1
d
N
19
/
1
d
ATC
(
t
)
,
N
19
/
2
d
(
t
+
1
)
=
A
19
/
2
d
N
19
/
2
d
(
t
)
+
B
19
/
2
d
τ
N
19
/
2
d
ATC
(
t
)
+
τ
q
16
/
2
d
(
t
)
,
q
19
/
2
d
(
t
)
=
C
19
/
2
d
N
19
/
2
d
(
t
)
-
B
19
/
2
d
N
19
/
2
d
ATC
(
t
)
.


### 4.2. Simulation Design

In this simulation sample there are 6 entry points for arrivals of the terminal network: MUBEL, OSIKA, BIPOP, P101, LONGTANG, and DANZHU. We assume that the arrival rate of each entry point obeys the negative exponential distribution [26]. To plug the average arrival rate that comes from flight historical statistics as expected value into distribution functions we can obtain the changes of traffic flux over time at each entry point as shown in Figure 6. It should be noted that we use equivalent traffic flow in this simulation and take each simulation time step as 1 minute.Figure 6
Changes of traffic flow over time at entry points.It can be noticed that traffic flux changes significantly over time at each entry point. When aircraft keep going to convergent points, the peak of traffic wave may happen to meet another one, in which situation traffic density nearby will be too high to meet separation requirements and the probability of unsafe events will increase accordingly. Instead, the trough of traffic wave may also meet another trough that may lead to a low traffic density and a large aircraft interval at convergent points, thus reducing the time/space utilization of limited airspace. Based on the situation we designed a “valve”. Agent aims at balancing traffic flow in the simulator system to manage the inflow/outflow of adjacent cells. In addition, this is much like the principle of “Cutting peak and filling valley” in real air traffic flow management [27, 28].The basic strategy of arrival “valve” control is as follows: first determine whether the mean arrival aircraft interval is lower than the arrival ATC separationd
i
a
ATC in one cell at any time, which means whether the stock of aircraft in one cell exceeds. If yes then the excessive number of aircraft should operate with the assigned ATC velocity v
i
a
ATC and the rest do not change their speed, while if no then all aircraft can still follow the STAR. Since v
i
a
ATC is lower than normal velocity, the traffic flux will also be lower compared to no control situations. With outflow decreasing, the stock of aircraft in one cell will increase accordingly. At next time step inflow will be added to the original cell stock. If this total value still exceeds the standards, control measures should be taken like before. These cyclic steps keep going until sometime there comes a small inflow and the total value added with cell stock goes lower than the standards. In this circumstance all aircraft in the cell can operate with STAR and aircraft controlled before can exit cell normally. In the whole process, the total number of aircraft in one cell cannot exceed the safety value at any time; that is, mean interval should never be smaller than the safety separation.Take typical arrival segment “GYA-AGVOS (Segment5)” as an example and making a comparison of air traffic flux and stock in cells between before and after the “valve” control, it is easy to find that the inflow and outflow in this segment become smooth and steady when there is the “valve” control and the stock of aircraft in each cell always meets the safety requirement as shown in Figures7(a) and 7(b). The unit is flights/minute.Comparison of traffic flux and stock over time before and after control measures.
(a)
Before arrival control
(b)
After arrival control
(c)
Before departure control
(d)
After departure controlThe departure in terminal area consists of two major parts: taking off from runway and flying along departure routes. We designed two successive processes in this simulation accordingly called the airport surface part and the air routes part, which are complementary in departure process. Specifically, the airport surface part mainly consists of departure flights schedule generation, runway time occupation, take-off slots allocation, and aircraft ground holding. We assume that the generation of departure flights schedule obeys a negative exponential distribution. The number of scheduled departure flights per unit time is as shown in Figure8(a). The runway time occupation and take-off slots allocation are actually how the landing aircraft and taking-off aircraft make the most use of runway slot resources within the limited and dynamic runway operation capacity. Changes in the number of landing and taking-off aircraft per unit time are as shown in Figures 8(b) and 8(c). Since landing aircraft have the priority, taking-off aircraft that are unable to go by flight schedules will be postponed to take ground-holding processes. Changes in the number of ground holding aircraft over time are as shown in Figure 8(d).Changes of traffic flux and stock over time in airport surface part simulation.
(a)
Scheduling
(b)
Landing
(c)
Taking off
(d)
G-holdingIn the process of flying along departure routes, there is no need to take control measures under the condition of full isolation between arrival and departure routes. The comparison of departure flux and stock over time before and after control measures is as shown in Figures7(c) and 7(d). The basic strategy of departure “valve” control is as follows: first determine whether the mean arrival aircraft interval is lower than the ATC separation d
i
ATC in one cell at any time. If yes then all the departures in this cell should take control measures, which means passing the segment with constant ATC velocity u
i
ATC. If no then determine whether the mean interval of both arrivals and departures is lower than ATC separation d
i
ATC in one cell, and if yes departures must be adjusted to increase the entire mean interval until it goes above d
i
ATC.


### 4.3. Result Analysis

According to the export data derived from NetLogo simulation platform we got series of parameter scatter diagrams that could reflect the basic traffic flow characteristics through statistic and analysis, as shown in Figure9. Based on this, we analyzed the mutual influence relationship among traffic flux q, density k, and velocity v of the arrival and departure routes in airport terminal area. Limited by space this paper only made detailed discussion on q
-
k relation of arrivals and departures; meanwhile the rest of v
-
k and q
-
v relations were listed in the form of statistical diagrams for reference.Scatter diagrams of basic arrival/departure flow characteristic parameters in terminal area.
(a)
(b)Taking typical arrival segment, Segment5, as an example, we derived the basic tendency of traffic flux and density relationship for arrival routes in terminal area as shown in Figure9(a1). From the diagram we can find that the relationship tendency consists of three main stages.Stage I is the free flow state in which the number of aircraft in segment is very low and the mean interval exceeds arrival ATC separationd
5
a
ATC, which means all the aircraft can follow STAR. The traffic flux is directly proportional to density in segment and the proportionality coefficient is equal to the mean velocity of segment. Stage II is the congestion flow state in which the number of aircraft in segment increases and then the mean interval is under arrival ATC separation d
5
a
ATC, which means part of the aircraft should be assigned to take control measures including decelerating, maneuvering or holding, and so forth. The relationship between traffic flux and density occurs a inflection point. The traffic flux in segment is still directly proportional to density; however the new proportionality coefficient is equal to the space mean value of the controlled aircraft velocity and uncontrolled aircraft velocity; that is,D
N
5
a
(
N
5
a
ATC
¯
t
5
s
+
∑
i
=
1
N
5
a
ATC
t
i
)
-
1, see details in Section 2.1(4).All of the control measures in the arrival process, including deceleration, maneuvering and holding, will lead to a decline in the displacement velocity. In the diagram it made the tendency of flux and density relationship to be leveling off. While the traffic flux will still increase with density, which differs from normal vehicle traffic since after inflection point the vehicle traffic flux decreases with density. The reason why this difference exists is that the car-following behavior between adjacent vehicles has significant influences when congestion occurs on the road while air traffic normally maintains a larger safe separation. Apart from this, some control measures as hold pattern make aircraft deviate from original air routes, which does not affect the other aircraft. Therefore, after the inflection point traffic flux will not decrease but increase with a low slope. Stage III is the block flow state in which the number of aircraft in segment exceeds safe value with the control adjustment in congestion state and large amount of traffic converging continuously. The mean interval is under the safe separationd
5
safe and insecurity factors surge that should be avoided as much as possible.From the above theoretical derivation we can find that variables related to operational flight program including initial velocityv
i
′, acceleration a
i, and length Ω
i of segment (cell) become constant values as the STAR/SID are established. This paper focuses on the rest of the variables especially flow control variables including ATC separation d
i
ATC and ATC velocity v
i
ATC. Such variables will become order parameters that influence the mutual relationships among three basic air traffic flow parameters.As shown in Figure9(a2), when arrival ATC separation d
5
a
ATC rises to 25% the inflection point of traffic flux and density relationship moves forward and actually the horizontal axis value equals 1
/
d
5
a
ATC. The rise of ATC standard will lead to a decrease of free flow state in Stage I. More aircraft need to take control measures; meanwhile the frequency of high traffic density even block flow in Stage III also increases. But in whatever Stage I or Stage III, the slope of traffic flux and density relationship maintains constant; that is, before inflection point the increase tendency coincides and after inflection point the increase tendency still keeps parallel. As shown in Figure 9(a3), when arrival ATC separation d
5
a
ATC remains the same and arrival ATC velocity v
5
ATC turns down, the inflection point stays the same but aircraft pass the segment with a lower mean velocity after this point. Therefore, traffic flux has a decelerated growth with the density. Conversely, if it needs to let traffic flux reach the level before changing of ATC velocity v
5
ATC, the density of air segment should be higher, which may lead to an early arrival of the block flow state in Stage III. The free flow state in Stage I has no change now since it is not affected by control measures.Taking typical departure segment, Segment19, as an example, we derived the basic tendency of traffic flux and density relationship for departure routes in terminal area as shown in Figure9(b1). The relationship tendency also consists of three main stages. Similar to arrival parts, Stage I is the free flow state in which the traffic flux of departures is directly proportional to density in segment and the proportionality coefficient is equal to the mean velocity of segment from SID. Stage II is the congestion flow state which is unlike arrival parts, since in departure routes we focus on the entire mean interval of both arrivals and departures, not just departures. With the density of departures rising up in segment, the entire mean interval can still meet the ATC separation d
19
ATC if arrivals density is small enough. Then the tendency of traffic flux and density relationship lies on the extension line of Stage I. Stage III is the block flow state in which the number of arrivals and the number of departures increase simultaneously. Accordingly the entire mean interval becomes lower than the ATC separation d
19
ATC. Considering the principle of arrivals priority part of the departures will accelerate to leave the heavy-traffic segment in order to release the space resources for arrivals. The new slope of traffic flux and density relationship equals the departure ATC velocity u
19
ATC. The number of adjusted departures is determined by arrivals, which is shown as a series of scatter values between the tendency line of Stage III and extension line of Stage II. The vertical distance from tendency line of Stage III to the scatter values varies inversely to the density of arrivals, that is, (
(
1
-
d
19
ATC
ρ
19
a
(
t
)
)
/
d
19
ATC
)
K (K is constant coefficient); see details in (27).As shown in Figure9(b2), after raising up 25% of the ATC separation d
19
ATC, there appear more departure aircraft that need to take control measures. Thus the Stage III also appears earlier while the slopes of tendency lines in congestion and block flow states have no change. However, there are more scatter values between these two lines. The Stage II congestion flow state in which entire mean interval (high departure, low arrival) still exceeds the ATC separation will reduce its frequency-of-occurrence. More scatter values lie in the Stage III block flow state. Conversely reducing the ATC separation d
19
ATC departures may not need to be adjusted in most cases, which tend to be the Stage I free flow state. More scatter values of traffic flux and density lie more on the extension line of free flow state.Keeping the ATC separationd
19
ATC unchanged and changing the departure ATC velocity u
19
ATC we can get new relationship tendencies as shown in Figure 9(b3). The Stage I free flow state keeps the same. But in Stage III block flow state the variation rate of traffic flux with density increases with the departure ATC velocity u
19
ATC, which means departure flow has raised its outflow rate of heavy-traffic segment. Meanwhile the impact of arrivals increases by T
ρ
19
a
(
t
)
u
19
ATC (T is constant coefficient); see details in (27), which is shown as the scatter values spread a larger scope from center line of the free flow state extension direction.


## 4.1. Simulation Sample

Based on the simulation platform NetLogo [23–25] we focused on each cell in the macroscopic model and designed traffic inflow/outflow behaviors to be each agent, so as to take control of the inflow/outflow volumes among all the cells. Applying the STAR/SID procedures to runway 02L of Guangzhou Baiyun International Airport (as shown in Figure 4) into the Netlogo system dynamic simulator (as shown in Figure 5) we can emulate the operation and evolution process of air traffic flow in the airport terminal area.Figure 4
STAR/SID routes to RWY02L ZGGG.Figure 5
NetLogo system dynamic simulator.The network of arrival and departure routes in terminal can be separated into several single direction segments and converging or diverging segments that can be further decomposed into more unit cells. For the purpose of simplicity, we set one cell dimension to be equal to the length of the shortest segment in the network and the other air segments to be integer multiples of unit cell dimension. According to this, we made the following Tables1 and 2 of the terminal network.Table 1
Cell quantities of arrival routes to RWY02L ZGGG.
Segment
Code
Quantity
MUBEL-GYA
Segment1
3 cells
BIPOP-GYA
Segment2
2 cells
OSIKA-GYA
Segment3
3 cells
P101-GYA
Segment4
2 cells
GYA-AGVOS
Segment5
2 cells
GYA-D5.1POU
Segment6
2 cells
LONGTANG-TAN
Segment7
2 cells
DANZHU-TAN
Segment8
2 cells
TAN-AGVOS
Segment9
2 cells
AGVOS-D15.7IOO
Segment10
1 cell
AGVOS-D5.1POU
Segment11
1 cell
D5.1POU-D15.7IOO
Segment12
1 cell
D15.7IOO-RWY02L
Segment13
1 cellTable 2
Cell quantities of departure routes to RWY02L ZGGG.
Segment
Code
Quantity
RWY02L-D12.0IOO
Segment14
1 cell
D12.0IOO-YIN
Segment15
2 cells
D12.0IOO-TAN
Segment16
1 cell
TAN-D18.2TAN
Segment17
2 cells
D18.2TAN-POU
Segment18
3 cells
TAN-AGVOS
Segment19
2 cells
AGVOS-D5.1POU
Segment20
1 cell
D5.1POU-POU
Segment21
2 cellsEach yellow rectangle in Figure4 represents one cell, which means a “stock” of fluid in the simulator system. Each grey pipe represents flowing between adjacent cells, which means “flow” in system, and the direction of arrow is the same as the flow direction. Each black valve represents control measures, which means “strategy” in system, and the inflow and outflow of air traffic will be under control by giving the certain valves some corresponding rules. Take typical overlapped segments “TAN=AGVOS (Segment9&19)” as an example, the cell transmission model of this is like (30). In the model, N
i
/
j
a
A
i
/
j
a
B
i
/
j
a
C
i
/
j
a
, respectively, represents the number of arrival aircraft and the simplified coefficients of cell j in Segment i. N
i
/
j
a
ATC denotes the number of arrival aircraft taking control measures in cell j of Segment i. q
i
/
j
a
(
t
) means the number of arrival aircraft that exit cell j in Segment i per unit time. N
i
/
j
d
A
i
/
j
d
B
i
/
j
d
C
i
/
j
d
N
i
/
j
d
ATC
q
i
/
j
d
(
t
)
, respectively, represents the similar corresponding ones but specially for departures. Coefficient k stands for the proportion of traffic flow into some segment when diverging occurs in air routes. τ is the designed time step for the simulator system:
(30)
N
9
/
1
a
(
t
+
1
)
=
A
9
/
1
a
N
9
/
1
a
(
t
)
+
B
9
/
1
a
τ
N
9
/
1
a
ATC
(
t
)
+
τ
[
q
7
/
2
a
(
t
)
+
q
8
/
2
a
(
t
)
]
,
q
9
/
1
a
(
t
)
=
C
9
/
1
a
N
9
/
1
a
(
t
)
-
B
9
/
1
a
N
9
/
1
a
ATC
(
t
)
,
N
9
/
2
a
(
t
+
1
)
=
A
9
/
2
a
N
9
/
2
a
(
t
)
+
B
9
/
2
a
τ
N
9
/
2
a
ATC
(
t
)
+
τ
q
9
/
1
a
(
t
)
,
q
9
/
2
-
10
/
1
a
(
t
)
=
k
9
/
2
-
10
/
1
[
C
9
/
2
a
N
9
/
2
a
(
t
)
-
B
9
/
2
N
9
/
2
a
ATC
(
t
)
]
,
q
9
/
2
-
11
/
1
a
(
t
)
=
k
9
/
2
-
11
/
1
[
C
9
/
2
a
N
9
/
2
a
(
t
)
-
B
9
/
2
N
9
/
2
a
ATC
(
t
)
]
,
(31)
N
19
/
1
d
(
t
+
1
)
=
A
19
/
1
d
N
19
/
1
d
(
t
)
+
B
19
/
1
d
τ
N
19
/
1
d
ATC
(
t
)
+
τ
q
16
/
1
d
(
t
)
,
q
19
/
1
d
(
t
)
=
C
19
/
1
d
N
19
/
1
d
(
t
)
-
B
19
/
1
d
N
19
/
1
d
ATC
(
t
)
,
N
19
/
2
d
(
t
+
1
)
=
A
19
/
2
d
N
19
/
2
d
(
t
)
+
B
19
/
2
d
τ
N
19
/
2
d
ATC
(
t
)
+
τ
q
16
/
2
d
(
t
)
,
q
19
/
2
d
(
t
)
=
C
19
/
2
d
N
19
/
2
d
(
t
)
-
B
19
/
2
d
N
19
/
2
d
ATC
(
t
)
.


## 4.2. Simulation Design

In this simulation sample there are 6 entry points for arrivals of the terminal network: MUBEL, OSIKA, BIPOP, P101, LONGTANG, and DANZHU. We assume that the arrival rate of each entry point obeys the negative exponential distribution [26]. To plug the average arrival rate that comes from flight historical statistics as expected value into distribution functions we can obtain the changes of traffic flux over time at each entry point as shown in Figure 6. It should be noted that we use equivalent traffic flow in this simulation and take each simulation time step as 1 minute.Figure 6
Changes of traffic flow over time at entry points.It can be noticed that traffic flux changes significantly over time at each entry point. When aircraft keep going to convergent points, the peak of traffic wave may happen to meet another one, in which situation traffic density nearby will be too high to meet separation requirements and the probability of unsafe events will increase accordingly. Instead, the trough of traffic wave may also meet another trough that may lead to a low traffic density and a large aircraft interval at convergent points, thus reducing the time/space utilization of limited airspace. Based on the situation we designed a “valve”. Agent aims at balancing traffic flow in the simulator system to manage the inflow/outflow of adjacent cells. In addition, this is much like the principle of “Cutting peak and filling valley” in real air traffic flow management [27, 28].The basic strategy of arrival “valve” control is as follows: first determine whether the mean arrival aircraft interval is lower than the arrival ATC separationd
i
a
ATC in one cell at any time, which means whether the stock of aircraft in one cell exceeds. If yes then the excessive number of aircraft should operate with the assigned ATC velocity v
i
a
ATC and the rest do not change their speed, while if no then all aircraft can still follow the STAR. Since v
i
a
ATC is lower than normal velocity, the traffic flux will also be lower compared to no control situations. With outflow decreasing, the stock of aircraft in one cell will increase accordingly. At next time step inflow will be added to the original cell stock. If this total value still exceeds the standards, control measures should be taken like before. These cyclic steps keep going until sometime there comes a small inflow and the total value added with cell stock goes lower than the standards. In this circumstance all aircraft in the cell can operate with STAR and aircraft controlled before can exit cell normally. In the whole process, the total number of aircraft in one cell cannot exceed the safety value at any time; that is, mean interval should never be smaller than the safety separation.Take typical arrival segment “GYA-AGVOS (Segment5)” as an example and making a comparison of air traffic flux and stock in cells between before and after the “valve” control, it is easy to find that the inflow and outflow in this segment become smooth and steady when there is the “valve” control and the stock of aircraft in each cell always meets the safety requirement as shown in Figures7(a) and 7(b). The unit is flights/minute.Comparison of traffic flux and stock over time before and after control measures.
(a)
Before arrival control
(b)
After arrival control
(c)
Before departure control
(d)
After departure controlThe departure in terminal area consists of two major parts: taking off from runway and flying along departure routes. We designed two successive processes in this simulation accordingly called the airport surface part and the air routes part, which are complementary in departure process. Specifically, the airport surface part mainly consists of departure flights schedule generation, runway time occupation, take-off slots allocation, and aircraft ground holding. We assume that the generation of departure flights schedule obeys a negative exponential distribution. The number of scheduled departure flights per unit time is as shown in Figure8(a). The runway time occupation and take-off slots allocation are actually how the landing aircraft and taking-off aircraft make the most use of runway slot resources within the limited and dynamic runway operation capacity. Changes in the number of landing and taking-off aircraft per unit time are as shown in Figures 8(b) and 8(c). Since landing aircraft have the priority, taking-off aircraft that are unable to go by flight schedules will be postponed to take ground-holding processes. Changes in the number of ground holding aircraft over time are as shown in Figure 8(d).Changes of traffic flux and stock over time in airport surface part simulation.
(a)
Scheduling
(b)
Landing
(c)
Taking off
(d)
G-holdingIn the process of flying along departure routes, there is no need to take control measures under the condition of full isolation between arrival and departure routes. The comparison of departure flux and stock over time before and after control measures is as shown in Figures7(c) and 7(d). The basic strategy of departure “valve” control is as follows: first determine whether the mean arrival aircraft interval is lower than the ATC separation d
i
ATC in one cell at any time. If yes then all the departures in this cell should take control measures, which means passing the segment with constant ATC velocity u
i
ATC. If no then determine whether the mean interval of both arrivals and departures is lower than ATC separation d
i
ATC in one cell, and if yes departures must be adjusted to increase the entire mean interval until it goes above d
i
ATC.


## 4.3. Result Analysis

According to the export data derived from NetLogo simulation platform we got series of parameter scatter diagrams that could reflect the basic traffic flow characteristics through statistic and analysis, as shown in Figure9. Based on this, we analyzed the mutual influence relationship among traffic flux q, density k, and velocity v of the arrival and departure routes in airport terminal area. Limited by space this paper only made detailed discussion on q
-
k relation of arrivals and departures; meanwhile the rest of v
-
k and q
-
v relations were listed in the form of statistical diagrams for reference.Scatter diagrams of basic arrival/departure flow characteristic parameters in terminal area.
(a)
(b)Taking typical arrival segment, Segment5, as an example, we derived the basic tendency of traffic flux and density relationship for arrival routes in terminal area as shown in Figure9(a1). From the diagram we can find that the relationship tendency consists of three main stages.Stage I is the free flow state in which the number of aircraft in segment is very low and the mean interval exceeds arrival ATC separationd
5
a
ATC, which means all the aircraft can follow STAR. The traffic flux is directly proportional to density in segment and the proportionality coefficient is equal to the mean velocity of segment. Stage II is the congestion flow state in which the number of aircraft in segment increases and then the mean interval is under arrival ATC separation d
5
a
ATC, which means part of the aircraft should be assigned to take control measures including decelerating, maneuvering or holding, and so forth. The relationship between traffic flux and density occurs a inflection point. The traffic flux in segment is still directly proportional to density; however the new proportionality coefficient is equal to the space mean value of the controlled aircraft velocity and uncontrolled aircraft velocity; that is,D
N
5
a
(
N
5
a
ATC
¯
t
5
s
+
∑
i
=
1
N
5
a
ATC
t
i
)
-
1, see details in Section 2.1(4).All of the control measures in the arrival process, including deceleration, maneuvering and holding, will lead to a decline in the displacement velocity. In the diagram it made the tendency of flux and density relationship to be leveling off. While the traffic flux will still increase with density, which differs from normal vehicle traffic since after inflection point the vehicle traffic flux decreases with density. The reason why this difference exists is that the car-following behavior between adjacent vehicles has significant influences when congestion occurs on the road while air traffic normally maintains a larger safe separation. Apart from this, some control measures as hold pattern make aircraft deviate from original air routes, which does not affect the other aircraft. Therefore, after the inflection point traffic flux will not decrease but increase with a low slope. Stage III is the block flow state in which the number of aircraft in segment exceeds safe value with the control adjustment in congestion state and large amount of traffic converging continuously. The mean interval is under the safe separationd
5
safe and insecurity factors surge that should be avoided as much as possible.From the above theoretical derivation we can find that variables related to operational flight program including initial velocityv
i
′, acceleration a
i, and length Ω
i of segment (cell) become constant values as the STAR/SID are established. This paper focuses on the rest of the variables especially flow control variables including ATC separation d
i
ATC and ATC velocity v
i
ATC. Such variables will become order parameters that influence the mutual relationships among three basic air traffic flow parameters.As shown in Figure9(a2), when arrival ATC separation d
5
a
ATC rises to 25% the inflection point of traffic flux and density relationship moves forward and actually the horizontal axis value equals 1
/
d
5
a
ATC. The rise of ATC standard will lead to a decrease of free flow state in Stage I. More aircraft need to take control measures; meanwhile the frequency of high traffic density even block flow in Stage III also increases. But in whatever Stage I or Stage III, the slope of traffic flux and density relationship maintains constant; that is, before inflection point the increase tendency coincides and after inflection point the increase tendency still keeps parallel. As shown in Figure 9(a3), when arrival ATC separation d
5
a
ATC remains the same and arrival ATC velocity v
5
ATC turns down, the inflection point stays the same but aircraft pass the segment with a lower mean velocity after this point. Therefore, traffic flux has a decelerated growth with the density. Conversely, if it needs to let traffic flux reach the level before changing of ATC velocity v
5
ATC, the density of air segment should be higher, which may lead to an early arrival of the block flow state in Stage III. The free flow state in Stage I has no change now since it is not affected by control measures.Taking typical departure segment, Segment19, as an example, we derived the basic tendency of traffic flux and density relationship for departure routes in terminal area as shown in Figure9(b1). The relationship tendency also consists of three main stages. Similar to arrival parts, Stage I is the free flow state in which the traffic flux of departures is directly proportional to density in segment and the proportionality coefficient is equal to the mean velocity of segment from SID. Stage II is the congestion flow state which is unlike arrival parts, since in departure routes we focus on the entire mean interval of both arrivals and departures, not just departures. With the density of departures rising up in segment, the entire mean interval can still meet the ATC separation d
19
ATC if arrivals density is small enough. Then the tendency of traffic flux and density relationship lies on the extension line of Stage I. Stage III is the block flow state in which the number of arrivals and the number of departures increase simultaneously. Accordingly the entire mean interval becomes lower than the ATC separation d
19
ATC. Considering the principle of arrivals priority part of the departures will accelerate to leave the heavy-traffic segment in order to release the space resources for arrivals. The new slope of traffic flux and density relationship equals the departure ATC velocity u
19
ATC. The number of adjusted departures is determined by arrivals, which is shown as a series of scatter values between the tendency line of Stage III and extension line of Stage II. The vertical distance from tendency line of Stage III to the scatter values varies inversely to the density of arrivals, that is, (
(
1
-
d
19
ATC
ρ
19
a
(
t
)
)
/
d
19
ATC
)
K (K is constant coefficient); see details in (27).As shown in Figure9(b2), after raising up 25% of the ATC separation d
19
ATC, there appear more departure aircraft that need to take control measures. Thus the Stage III also appears earlier while the slopes of tendency lines in congestion and block flow states have no change. However, there are more scatter values between these two lines. The Stage II congestion flow state in which entire mean interval (high departure, low arrival) still exceeds the ATC separation will reduce its frequency-of-occurrence. More scatter values lie in the Stage III block flow state. Conversely reducing the ATC separation d
19
ATC departures may not need to be adjusted in most cases, which tend to be the Stage I free flow state. More scatter values of traffic flux and density lie more on the extension line of free flow state.Keeping the ATC separationd
19
ATC unchanged and changing the departure ATC velocity u
19
ATC we can get new relationship tendencies as shown in Figure 9(b3). The Stage I free flow state keeps the same. But in Stage III block flow state the variation rate of traffic flux with density increases with the departure ATC velocity u
19
ATC, which means departure flow has raised its outflow rate of heavy-traffic segment. Meanwhile the impact of arrivals increases by T
ρ
19
a
(
t
)
u
19
ATC (T is constant coefficient); see details in (27), which is shown as the scatter values spread a larger scope from center line of the free flow state extension direction.


## 5. Conclusions

Based on CTM we have proposed the macroscopic model of air traffic flow in airport terminal area and carried out a series of simulation experiments with the NetLogo platform. Through both of the theoretical and practical discussions, we could generally reveal the basic interrelationships and influential factors of air traffic flow characteristic parameters. The research findings are as follows.(1)
The CTM could accurately reflect the macroevolution laws of air traffic flow in terminal air route network. Meanwhile, it may also be applied to different kinds of air traffic scenes including airways, sectors, or airspaces by modifying conditions correspondingly. Moreover, it has a remarkable operational efficiency in multiagent simulations.
(2)
There are obvious relationships among the three characteristic parameters as flux, density, and velocity of the air traffic flow in terminal area. And such relationships evolve distinctly with the flight procedures, control separations, and ATC strategies. The air traffic flow characteristics may take the specific changes through flight procedure optimization, control separation modification, or ATC strategy regulation, which could be part of the scientific basis for air traffic management in airport terminal area.
(3)
The default parameters we used in simulation experiments are from practical ATC rules. Automatic optimization of these parameters for desired traffic flow characteristics should be taken into consideration in further studies. In addition, discussions in this paper focus on macro perspectives, thus the research results seem rough anyway. To obtain more detailed and elaborate traffic flow characteristics it is necessary to combine such macro studies with micro perspectives that give full expressions to the following: overflying or turning and so forth of individual behaviors and interactive effects. It also should be an important direction for further study.



---

*Source: 741654-2014-08-25.xml*