Traffic system operation optimization incorporating buffer size

JID:AESCTE AID:3948 /FLA

[m5G; v1.208; Prn:16/03/2017; 10:43] P.1 (1-12)

Aerospace Science and Technology ••• (••••) •••–•••

1

Contents lists available at ScienceDirect

2

67 68

3

Aerospace Science and Technology

4 5

69 70 71

6

72

www.elsevier.com/locate/aescte

7

73

8

74

9

75

10

76

11 12 13 14 15 16

Traffic system operation optimization incorporating buffer size

b

78 79

Yun-xiang Han a , Xiao-qiong Huang b , Yan Zhang a a

77

80

School of Automobile and Traffic Engineering, Jiangsu University of Technology, Changzhou 213000, China School of business, Jiangsu University of Technology, Changzhou 213000, China

81 82

17

83

18

84

19

a r t i c l e

i n f o

85

a b s t r a c t

20 21 22 23 24

86

Article history: Received 21 June 2016 Received in revised form 4 January 2017 Accepted 7 March 2017 Available online xxxx

25 26 27 28 29 30

Keywords: Air transportation Air traffic control Mathematical models Buffer size Discrete event system

31

This paper presents a general framework for resolving resource utilizations conflicts of air traffic system, expressed in the form of a max-plus linear model. The general air traffic system optimization model is presented taking into account buffer size in view of different purposes. The proposed air traffic system model is very flexible and it is extensible focusing on distinct circumstances. The dynamics of air traffic systems is characterized through the occurrence of discrete events such as aircrafts entering or leaving sub-segments. We focus on input control of air traffic system modeled by max-plus algebra by controlling the inflow of aircraft into the system. The constraints between input variables, state variables and output variables were obtained. The proposed method aims to minimize system total delay to meet demand specifications, allowing the air traffic controller to obtain the best control policy, which delays the occurrence of input events or varying input rate of system resources. Through simulations we verify the performance of the proposed max-plus linear model control scheme. © 2017 Elsevier Masson SAS. All rights reserved.

87 88 89 90 91 92 93 94 95 96 97

32

98

33

99

34

100

35 36

101

1. Introduction

37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59

The steady growth of air traffic all over the world leads to complicated air traffic system situations where the aircrafts may infringe safety separation standards in airspace with dense traffic flow. Air traffic management can provide effective decision tool for airspace users according to the preferred flight profiles. Specifically, current approach consists of two main activities aimed at increasing airspace capacity and decreasing air traffic controller workloads, that is, air traffic control and air traffic flow management. Additionally, the increase of air traffic volume urges to improve the efficiency of air traffic system, which was attested by the Next Generation Air Transportation System (NGATS) and the Single European Sky Air Traffic Research System (SESAR) aimed at highlighting the importance of trajectory generation in achieving the targets of increasing capacity and safety while improving flight efficiency [31,38]. The conflict detection and resolution has been a major topic in air traffic management for decades. Future position of aircraft in the airspace can be obtained using aircraft state information available such as position and altitude, as well as trend vector so that potential flight conflicts can be given by predefining metrics. The metric usually include a sole parameter such as distance or a combination of several parameters. The term ‘trajectory optimization’ is referred as actions taken when the aircraft depar-

60 61 62 63 64 65 66

E-mail addresses: [email protected] (Y.-x. Han), [email protected] (X.-q. Huang). http://dx.doi.org/10.1016/j.ast.2017.03.012 1270-9638/© 2017 Elsevier Masson SAS. All rights reserved.

ture time is known or even after the flight is airborne but with sufficient time to allow a collaborative decision-making process [32]. This term excludes instructions and clearances that require an immediate response. Over the past years, various optimization models for resolving flight conflict have been put forward. Typically, these methods can be categorized into nominal model, worst case model and probabilistic model [23]. For the nominal method, current aircraft states are projected directly into future [21]. The worst case assumes that an aircraft will perform many kinds of maneuvers. If any one of these maneuvers cause a flight conflict, then a potential conflict is predicted [11]. For the probabilistic method, uncertainties are modeled to describe trajectory uncertainties [17,36]. In addition, aircrafts can use a set of optimization maneuvers stem from optimal control model, heuristics model or multi-agent model in order to solve a flight conflict. Specifically, the most realistic models are developed in the theoretical framework of optimal control [4,7,12, 17,20,22,26,35,36,41–43]. As a consequence, the continuous-time nature of the flight conflict resolution problem is conserved, but analytical solutions can be found only for the special case. Another interesting approach for en-route flight conflict resolution is heuristic algorithms [2,27,39,40]. The genetic optimization algorithm is developed based on genetic encoding, where crossover and mutation are introduced in each generation. Further developments are presented combining speed and altitude changes using mixed-integer linear optimization model [1,19, 30]. The logical constraints are used to denote the configuration requirements for maneuvers and the collision avoidance problems

102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

2

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41

AID:3948 /FLA

[m5G; v1.208; Prn:16/03/2017; 10:43] P.2 (1-12)

Y.-x. Han et al. / Aerospace Science and Technology ••• (••••) •••–•••

can be combined into a single mixed-integer linear optimization model. Remaining optimization methods for flight conflict resolution rely on game theory [4,8,11,40]. This approach translates safety specifications into restrictions on the system state reachable sets and the Hamilton–Jacobi equations whose solutions describe the boundaries of reachable sets can be acquired using game theory. In other contexts, air traffic management system could also be supposed as a multi-agent system in which different aircrafts working together in order to keep safe. Consequently, distributed artificial intelligence methods and multi-agent system technologies for flight conflict resolution are also put forward. The presented methods involve a set of intelligent agents aiming at achieving some specific goals through negotiation with each other [5,18,24, 25,34,37]. High performance flight conflict planning methodologies can increase both airspace efficiency and safety by enforcing the separation activities between aircrafts, making it possible to integrate strategic approach into air traffic system instead of tactical process [33]. Thus, there is an urgent need to develop a novel strategic planning method for air traffic system in order to utilize the air traffic system resource efficiently. To achieve this, air traffic system operation is planned at different time frames [9,10]. Strategic planning is done several days before take-off and consists of assigning flight plans for a whole day. Our work focuses on air traffic system strategic planning. In particular, we studied the maxplus linear systems with time-invariant transition structures. For this class of systems, we build a multiple aircraft trajectory optimization model based on the dynamic of air traffic flow passing through the airspace network. The outline of this paper is as follows. Firstly, the max-plus theory is introduced and air traffic system single jet route max-plus model with finite buffer size is presented in section 2. Next, air traffic system general linear state-space representation in max-plus algebra is put forward. Section 3 considers the constraints of air traffic system with multiple jet routes. The multi-aircraft trajectory optimization cost function and case studies are demonstrated in section 4 and section 5, respectively. Finally, section 6 gives some conclusions.

44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

In general, discrete event systems that model only synchronization aspects are called max-plus-linear systems consisting of manmade systems that contain a finite number of resources, which are shared by several users all of which contribute to the achievement of some common goals [15]. Air traffic system is a typical example of discrete event system essentially consists of air traffic controller, aircraft and airspace network, which changes due to the occurrence of discrete events such as air traffic controllers’ instructions so that its behavior is governed by the progression clock ticks. The techniques available for conventional linear systems cannot be extended directly to discrete event systems. This constitutes the motivation for the work reported in this paper [28]. It is well known that max-plus theory play a key role among the modeling methodologies for discrete event system, since they are able to capture the precedence relations and interactions among the concurrent and asynchronous events that are typical in discrete event system. The linear properties of max-plus models that describe air traffic flow make control policies for the airspace very attractive. Attempts like this have been successfully made in many areas such as queuing system, manufacturing systems, telecommunication networks and railway networks [3,6,16,29]. In this paper, we will develop a max-plus-linear framework with finite buffer size for air traffic system and the proposed model can be used to control all resources of the airspace.

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97

x ⊕ y = max(x, y )

(1)

98 99 100

x⊗ y=x+ y

(2)

101

for numbers x, y ∈ R ε and we extend the (max, +) algebra operations to matrices in the following way:

102

[ A ⊕ B ]i j = ai j ⊕ bi j = max(ai j , bi j )

105

2. Modeling of the single jet route with finite buffer size

42 43

Air traffic system is a complicated dynamic system. In particular, the terminology buffer refers to en-route holding patterns or holding procedures notified in relation to an aerodrome. In addition, there are also some connections between buffer and segment or sector capacity for the air traffic system. It is well known that flight levels are used to ensure safe vertical separation between aircrafts, despite natural local variations in atmospheric air pressure. Flight levels solve this problem by defining altitudes based on a standard air pressure at sea-level. All aircrafts operating on flight levels calibrate to this setting regardless of the actual sea level pressure. The number of available flight levels corresponding to holding patterns or holding procedures is called finite buffer size, which will inevitably influence the aircraft trajectory optimization as indicated below. In the following section, we will present the modeling approach focusing on air traffic system with finite buffer size by max-plus algebra theory. We firstly give the basic definition of the max-plus algebra and present some results on a class of (max, +) functions. We adopt the convention that for all x ∈ R , max(x, −∞) = max(−∞, x) = x and x + (−∞) = −∞ + x = −∞. The two basic operations in the max-plus algebra can be represented by ⊕ and ⊗, respectively. In particular, ε is the neutral element for the operation ⊕ and absorbing for ⊗, that is, for all x ∈ R ε , x ⊗ ε = ε , ε = −∞, R ε = R ∪ {ε }. The identity element for the max-plus algebra is e = 0. Accordingly, the null matrix  and identity matrix E are defined  0, i = j as: φi j = ε (i = 1, 2, . . . , n; j = 1, 2, . . . , n) and E i j = ε, i = j (i = 1, 2, . . . , n; j = 1, 2, . . . , n), respectively. Furthermore, the basic operations addition (⊕) and multiplication (⊗) are defined as follows [15]:

[ A ⊗ C ]i j =

n  k =1

(3)

103 104 106 107

aik ⊗ ckj =

max (aik + ckj )

k=1,2,...,n

n× p

(4)

×n and C ∈ R for matrices A , B ∈ R m ε . ε We now show by an example how the single jet route with finite buffer size, characterized only by synchronization, can be modeled using max-plus algebra. To obtain the max-plus system model of air traffic system, we firstly decide the characteristic locations through which flight flow rates need to be determined. These characteristic locations consist of air route split or intersection point and airspace fixes. In general, aircraft climb profiles are defined in terms of a constant calibrated airspeed (CAS) segment and a constant Mach segment. Consider a serial air traffic system of multiple sub-segments with finite buffer sizes, aircrafts have to pass through the serial characteristic locations consecutively so as to receive service at each sub-segment. We now consider, for the sake of simplicity, a typical CAS/Mach climbing/descending jet route with finite buffer size as illustrated, for example, by A B / A  B  in Fig. 1. At a coarse grain, we could consider each sub-segment as a ‘machine’ on which to process multiple tasks with different process time (depending on the aircraft property). Let n be the number of sub-segment M 1 , M 2 , · · · , M n and m the number of aircraft P 1 , P 2 , . . . , P m , respectively. In addition, P 1 , P 2 , . . . , P m also represents an external arrival stream of aircrafts. Each sub-segment of the jet route has m aircrafts allocated

108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:3948 /FLA

[m5G; v1.208; Prn:16/03/2017; 10:43] P.3 (1-12)

Y.-x. Han et al. / Aerospace Science and Technology ••• (••••) •••–•••

3

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

90

Fig. 1. Division of single climbing/descending jet route.

to it. In other words, each sub-segment M j ( j = 1, 2, . . . , n) provides sequential service for each aircraft P i (i = 1, 2, . . . , m) and each aircraft P i (i = 1, 2, . . . , m) receives service provided by each sub-segment M j ( j = 1, 2, . . . , n). For simplicity, we assume that initially the system starts without a single aircraft. The way in which the aircrafts are processed at the sub-segments is called service activity. The most prominent example is the first come first served (FCFS) discipline and the internal mechanism of air traffic system can be forced by this discipline. Specifically, aircrafts are served in the order of their arrival using FCFS basis and this also implies that the aircrafts leave the sub-segments in the same order as they arrive at them. For the circumstances that aircraft starts to enter into a certain sub-segment or sub-segment starts to provide service, this is referred to as resource input. In contrast, we call the aircraft flows out of a particular sub-segment or sub-segment completes the service activity resource output. For the air traffic system modeling, the meaning of each variable is defined as follows: Definition xi j : the earliest time instant at which sub-segment M j ( j = 1, 2, . . . , n) starts to provide service for aircraft P i (i = 1, 2, . . . , m); Definition ul : time instant at which the lth (l = 1, 2, . . . , n + m) resource involved in the first service activity; Definition yl : the earliest time instant at which the lth (l = 1, 2, . . . , n + m) resource releases from the service; Definition t i j : process time during which sub-segment M j ( j = 1, 2, . . . , n) provides service for aircraft P i (i = 1, 2, . . . , m). As can be seen from Fig. 1, each sub-segment has one finite buffer b allocated to it and the maximum number of aircraft that each buffer can ‘store’ is finite, that is, C (b i ) (i = 1, 2, . . . , n) is finite. In the following section, we will assume that each buffer has the same size aimed at decreasing the complexity of system modeling and analysis, that is, C (b i ) = b (i = 1, 2, . . . , n). Moreover, we can associate the service activity with system variables or parameters using simple element correspondences as listed in Table 1 and Table 2. In this way, the state variables constraints of air traffic system with finite buffer size are established for four distinct scenarios as described below [13]. Besides, we give some fundamental preconditions relevant to system service before establishing the system model. One common

91 92

Table 1 Correspondences between system variables and service activities.

un+1 ( P 1 ) un+2 ( P 2 )

u1 (M1 )

u2 (M2 )

···

un ( M n )

x11 x21

x12 x22

··· ··· ··· ··· ···

x1n x2n

···

···

···

un+m ( P m )

xm1 y1 (M1 )

xm2 y2 (M2 )

93 94 95

yn+1 ( P 1 ) yn+2 ( P 2 )

···

···

xmn yn ( M n )

yn+m ( P m )

96 97 98 99 100 101

Table 2 Correspondences between system parameters and system resources.

102 103

M1

M2

···

Mn

104

P1 P2

t 11 t 21

t 12 t 22

t 1n t 2n

105

Pm

tm1

tm2

··· ··· ··· ···

···

···

···

···

tmn

106 107 108 109

feature of the jet route is that each sub-segment cannot start a fresh service activity until certain preceding sub-segments have all completed their service activities. The interactions between the air traffic system resources are governed by the following constraints: (1) For i = 1 and j = 1, it corresponds to sub-segment j is available and aircraft i is located at the entry of the sub-segment j. For the serial sub-segment j with buffer size b in the kth batch, x11 (k) which is associated with the following parameters corresponding to sub-segment M 1 (k) and aircraft P 1 (k), that is

110 111 112 113 114 115 116 117 118 119

(a) un+1 (k); (b) u 1 (k); (c) Time instant t α that the last aircraft P m (k − 1) in the previous batch departs from sub-segment M 1 (k − 1). As is well known, t α depends on the corresponding buffer status, namely, whether buffer of M 2 (k − 1) is full or not. If the buffer of M 2 (k − 1) is not filled, t α is equivalent to the departure time instant of aircraft P m (k − 1) from sub-segment M 1 (k − 1), that is, t α = xm1 (k − 1) + tm1 , where tm1 is the process time during which M 1 provides service for P m . Otherwise, t α is equivalent to the time instant that aircraft P m−α (k − 1) departs from buffer b2 and then enters to sub-segment M 2 (k − 1), as well as the time instant that aircraft P m (k − 1) enters to buffer b2 .

120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

AID:3948 /FLA

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

[m5G; v1.208; Prn:16/03/2017; 10:43] P.4 (1-12)

Y.-x. Han et al. / Aerospace Science and Technology ••• (••••) •••–•••

4

In other words, P 1 (k) cannot be served before the (m − α )th departure has taken place. In such cases, the total number of aircrafts { P m (k − 1), P m−1 (k − 1), . . . , P m−α (k − 1)} in buffer b2 is b and α = b − 1. Therefore, t α = xm−b+1,2 (k − 1). So, in summary, P 1 (k) can only be served if the following constraint are fulfilled:

20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41



 + tm1 , xm−b+1,2 (k − 1)

(5)

(2) For i = 1 and j = 1, it corresponds to sub-segment j is available and the completion of service for aircraft i provided by subsegment ( j − 1). For the serial sub-segment j with buffer size b in the kth batch, x j1 (k) ( j = 2, 3, . . . , m) which is associated with the following parameters corresponding to sub-segment M 1 (k) and aircraft P j (k), that is (a) un+ j (k); (b) Time instant t β that aircraft P j −1 (k) in the kth batch departs from sub-segment M 1 (k). As is well known, t β depends on the corresponding buffer status, namely, whether buffer of M 2 (k) is full or not. If the buffer of M 2 (k) is not filled, t β is equivalent to the departure time instant of aircraft P j −1 (k), that is, t β = x j −1,1 (k) + t j −1,1 , where t j −1,1 is the process time during which M 1 provides service for P j −1 . Otherwise, t β is equivalent to the time instant that aircraft P j −β (k − 1) departs from buffer b2 and then enters to sub-segment M 2 (k), as well as the time instant that aircraft P j −1 (k) enters to buffer b2 . In other words, P j (k) cannot be served before the ( j − β)th departure has taken place. In such cases, the total number of aircrafts { P j −1 , P j −2 , . . . , P j −β } in buffer b2 is b and β = b. Indeed, for j > b, P j −β = P j −b (k), t β = x j −b,2 (k), otherwise, P j −β = P m+ j −b (k − 1), t β = xm+ j −b,2 (k − 1). So, in summary, P j (k) can only be served if the following constraints are fulfilled:

 ⎧ x j1 (k) = max un+ j (k), x j −1,1 (k)  ⎪ ⎪ ⎨ + t j −1,1 , xm−b+ j ,2 (k − 1) , j ≤ b  ⎪ ⎪ ⎩ x j1 (k) = max un+ j (k), x j −1,1 (k) + t j −1,1 , x j −b,2 (k) , j > b

43 44 45 46 47 48

(3) For i = 1 and j = 1, it corresponds to the completion of service for aircraft (i − 1) provided by sub-segment j and aircraft i is located at the entry of the sub-segment j. For the serial subsegment j with buffer size b in the kth batch, x1i (k) which is associated with the following parameters corresponding to subsegment M i (k) and aircraft P 1 (k), that is

49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

(a) u i (k); (b) Time instant of P 1 (k) departs from sub-segment M i −1 (k), that is, x1,i −1 (k) + t 1,i −1 ; (c) Time instant t μ that the last aircraft P m (k − 1) in the previous batch departs from sub-segment M i (k − 1). As is well known, t μ depends on the corresponding buffer status, namely, whether the buffer of M i +1 (k − 1) is full or not. If the buffer of M i +1 (k − 1) is not filled, t μ is equivalent to the departure time instant of aircraft P m (k − 1), that is, t μ = xmi (k − 1) + tmi , where tmi is the process time during which M i provides service for P m . Otherwise, t μ is equivalent to the time instant that aircraft P m−μ (k − 1) departs from buffer B i +1 and then enters to sub-segment M i +1 (k − 1), as well as the time instant that aircraft P m (k − 1) enters to buffer b i +1 . In other words, P 1 (k) cannot be served before the (m − μ)th departure has taken place. In such cases, the total number of aircrafts { P m (k − 1), P m−1 (k − 1), . . . , P m−μ (k − 1)} in buffer

(7)

(4) For i = 1 and j = 1, it corresponds to the completion of service for aircraft (i − 1) provided by sub-segment j as well as the service for aircraft i provided by sub-segment ( j − 1). For the serial sub-segment j with buffer size b in the kth batch, x ji (k) which is associated with the following parameters corresponding to subsegment M i (k) and aircraft P j (k), that is (a) Time instant of P j (k) departs from sub-segment M i −1 (k), that is, x j ,i −1 (k) + t j ,i −1 ; (b) Time instant t ρ that aircraft P j −1 (k) in the previous batch departs from sub-segment M i (k). As is well known, t ρ depends on the corresponding buffer status, namely, whether buffer of M i +1 (k) is full or not. If the buffer of M i +1 (k) is not filled, t ρ is equivalent to the departure time instant of aircraft P j −1 (k), that is, t ρ = x j −1,i (k) + t j −1,i , where t j −1,i is the process time during which M i provides service for P j −1 . Otherwise, t ρ is equivalent to the time instant that aircraft P m−ρ (k − 1) departs from buffer b i +1 and then enters to sub-segment M i +1 (k − 1), as well as the time instant that aircraft P j −1 (k) enters to buffer b i +1 . In other words, P j (k) cannot be served before the (m − ρ )th departure has taken place. In such cases, the total number of aircrafts { P j −1 , P j −2 , . . . , P j −ρ } in buffer b i +1 is b and ρ = b. Indeed, for j > b, P j −ρ = P j −b (k), t ρ = x j −b,i +1 (k), otherwise, P j −ρ = P m+ j −b (k − 1), t ρ = xm+ j −b,i +1 (k − 1). So, in summary, P j (k) can only be served if the following constraints are fulfilled:

 ⎧ x ji (k) = max x j ,i −1 (k) + t j ,i −1 , x j −1,i (k) ⎪ ⎪ ⎨ + t j −1,i , xm+ j −b,i +1 (k − 1) , j ≤ b  ⎪ ⎪ ⎩ x ji (k) = max x j ,i −1 (k) + t j ,i −1 , x j −1,i (k) + t j −1,i , x j −b,i +1 (k) , j > b

72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 102

(8)

103 104 105

⎧ ⎨ y i (k) = tmi ⊗ xmi (k) ⊕ 0 ⊗ xm−b+1,i +1 (k) ( i = 1, 2, . . . , n ) ⎩ yn+ j (k) = t jn ⊗ x jn (k) ( j = 1, 2, . . . , m)

106 107 108 109

(9)

110 111 112

In the above descriptions, we associated input variables and output variables, as well as control variables with each system resource for the air traffic system service. If we write these dynamical equations described above in max-plus matrix notation, we can obtain the time evolution of the air traffic system. Combining equations (5)–(9) yields the following general air traffic system recursion model given below for the departure or arrival aircraft queues and it will provide the basis for our further analysis.



69

101

Moreover, from what has been discussed above, the air traffic system output variables can be given by:

42

68

71

 + tmi , xm−b+1,i +1 (k − 1)

(6)

67

70

x1i (k) = max u i (k), x1,i −1 (k) + t 1,i −1 , xmi (k − 1)



x11 (k) = max u 1 (k), un+1 (k), xm1 (k − 1)

18 19

b i +1 is b and μ = b − 1. Therefore, t μ = xm−b+1,i +1 (k − 1). So, in summary, P 1 (k) can only be served if the following constraint are fulfilled:

113 114 115 116 117 118 119 120 121 122

x(k) = Ax(k) ⊕ F x(k − 1) ⊕ Bu (k)

(10)

y (k) = C x(k)

123 124 125

for

126 127

T

u = [u 1 , u 2 , . . . , un , un+1 , un+2 , . . . , un+m ]

128

y = [ y 1 , y 2 , . . . , yn , yn+1 , yn+2 , . . . , yn+m ]T

129

x = [x11 , x12 , . . . , x1n , x21 , x22 , . . . , x2n , · · · , xm1 , xm2 , . . . , xmn ]

130

and

132

T

131

JID:AESCTE AID:3948 /FLA

[m5G; v1.208; Prn:16/03/2017; 10:43] P.5 (1-12)

Y.-x. Han et al. / Aerospace Science and Technology ••• (••••) •••–•••

1 2 3 4 5 6 7 8

⎡

⎢ A 22 ⎢ ⎢ .. .. ⎢ . . ⎢ A A bb A=⎢ b , b −1 ⎢ ⎢ A b+1,1 A b+1,b A b+1,b+1 ⎢ ⎢ .. .. ⎣ . .

9

⎡

10 11 12 13 14

⎢ ⎢ ⎢ F =⎢ ⎢ ⎣

15

⎡

16 17 18

⎡

20 22 23 24

⎢ ⎢ B =⎢ ⎣

25

⎡

26 27 28 29 30 31 32 33 34 35 36

A m,m−b

 F 1,m−b+1

F 1m F 2,m−b+2

..

 .. . 

.

⎤

⎤

..

0

ε 0 ⎢ ⎢ .. .. ⎢ ⎢ . . ⎢ ⎢ . .. 0 ⎢ ⎢ ε C =⎢ ⎢ ⎢ t 1n ε ⎢ ⎢ ⎢ t 2n · · · · · · · · · · · · · · · ⎢ ⎢ ⎣

⎤

tm1 tm2

41 42

⎥ ⎥ .. ⎥ ⎥ . ⎥ ⎥ .. ⎥ . ⎥ tmn ⎥ ⎥ ε ⎥ ⎥ .. ⎥ . ⎥ ⎥ .. ⎥ . ⎦ tmn

The constraints previously mentioned represent the time sequences of departure/arrival time instants of aircraft queues via a max-plus algebra. Another interesting application of the system model with finite buffer size is the conclusions given below:

43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61

Remark 1. The proposed air traffic system model reveals the inherent relationship various variables in the serial service. Furthermore, we can acquire system feedback control variable u (k) = K y (k − 1) in which K is feedback constant and the elements of K satisfy:

k = diag{k11 , k22 , . . . , kqq} (kii = ε , q = m + n)

(11)

That is, each sub-segment can be involved in the first service activity of the kth batch only after it is released from the (k − 1)th batch. In other words, the process for the next batch of aircrafts can only take place immediately after the completion of prior batch for a specified sub-segment. As a consequence, the closed loop linear model of the serial jet route can be expressed by:



x(k) = N x(k − 1) y (k) = N f x(k − 1) ⊕ M y (k − 1)

(12)

where N = A ∗ ( F ⊕ B K C ) ∈ R p × p , N f = C A ∗ F ∈ R q× p , M = C A ∗ B K ∈ R q×q , p = mn = dim(x), q = m + n = dim( y ) = dim(u ).

62 63 64 65 66

Remark 2. The variables t i j and xi j should be random variables when considering the influence of all kinds of random factors related to aircraft flight. In this way, we can transform the deterministic model into stochastic ones.

⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ E =⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎣

75 76 77 78

mn−1

⎤

0

..

.

ε 0 0

..

. 0

ε

..

. 0

..

79 80

⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎦

. 0

81 82 83 84 85 86 87 88 89 90 91 92

mn×mn

93

Remark 4. For the air traffic system model with finite buffer size, if b → 0, namely, the buffer size is zero, then the additional relation introduced in the original models will disappear eventually and it can transform into system model with none buffer, that is



x = A ⊗ x ⊕ B  ⊗ u y = C ⊗ x

¯ ⊗ x¯ (k − 1) ⊕ B¯ ⊗ u¯ (k) x¯ (k) = A ¯y = C¯ ⊗ x¯ (k)

94 95 96 97 98

(14)

Moreover, if we denote time instant at which the kth aircraft arrives at the entrance of sub-segment l1 and the earliest time instant at which the sub-segment li (i = 1, 2, . . . , n) starts to provide service for the kth aircraft by u (k) and xi (k), respectively. In addition, we denote time instant at which the kth aircraft flows out of the whole jet route and time period during which the kth aircraft travels through the sub-segment li (i = 1, 2, . . . , n) by y (k) and t i (k), respectively. The original system model with none buffer can be transformed into:



71

74

⎡

B 2m

69

73

⊕ ( A ⊗ A ⊗ A ⊗ · · · ⊗ A)   

.. ⎥ , . ⎦

68

72

A = E ⊕ A ⊕ ( A ⊗ A) ⊕ ( A ⊗ A ⊗ A) ⊕ · · ·

B 21 B 22 ⎥ ⎥

ε

(13)

∗

⎤

67

70

x(k) = A ∗ F x(k − 1) ⊕ A ∗ Bu (k) y (k) = C A ∗ F x(k − 1) ⊕ C A ∗ Bu (k)

where

⎥ ⎥ ⎥ ⎥, ⎥ F bm ⎦ 

38 40



A m,m−1 A mm

37 39

.

Remark 3. The analytical form of state vector x and output vector y in the proposed air traffic system max-plus model with finite buffer size can be expressed by:

⎥ ⎥ ⎥ ⎥ ⎥ ⎥, ⎥ ⎥ ⎥ ⎥ ⎦

.. ⎥ , .⎦

ε ··· ε B 11

..

.

ε ··· ε

⎢  = ⎣ ...

19 21

⎤

A 11 A 21

5

99 100 101 102 103 104 105 106 107 108 109 110 111

(15)

according to the system service principles as indicated above. As a consequent, we can also obtain the feedback control equation u (k) = K y (k − 1) in which K is feedback constant.

112 113 114 115 116 117

3. Modeling of the multiple jet routes with finite buffer size

118 119

Multiple jet routes with intersections air traffic flow max-plus model can be handled similarly compared with single jet route. In the following section, we will focus on the formulation of cross jet routes traffic flow max-plus model with finite buffer size and discuss the most typical scenario depicted in Fig. 2. It is easy to know that the flight conflict point lies in each split ¯ O B¯ and point near the intersection O for the cross jet routes A ¯  O B¯  . We now turn to the critical connection location for the two A cross jet routes. By regarding the intersection as a single-server node, then we can apply the FCFS queuing discipline to the different air traffic flow passing through O . Furthermore, we introduce a re-sequencing mechanism at the merge point O in order to obtain the constraint according to FCFS discipline. For the aircrafts

120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

AID:3948 /FLA

[m5G; v1.208; Prn:16/03/2017; 10:43] P.6 (1-12)

Y.-x. Han et al. / Aerospace Science and Technology ••• (••••) •••–•••

6

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44

instance, they can denote an average time or money cost of the delay for all the departures. In this case, it minimizes an average cost of total departure delays for the set of flights considered. Another application of the model coefficients is to use them as control parameters. The possibility to vary the parameters makes the model more realistic and more flexible in providing alternative solutions.

xh, O = u h, O ⊕ x g , O ⊗ t g , O

5. Worked example and discussion

(16)

It implies that the buffer size of O is zero. However, supposing that the buffer size of intersection O is b O , we can then re-sequence all the aircrafts on the basis of their initial arrival time at the merge point O following the equations (5)–(8) described above. That is, the subsequent aircrafts can only arrival at the re-sequencing node when all aircrafts that arrived at the node before have been served completely according to the rule of FCFS. Accordingly, we can obtain the analogical constrains as demonstrated before at O . Note that a basic airspace unit (such as a terminal control area) may involve a wide variety of flight conflicts. It is straightforward to construct air traffic system model adapted to the control of air traffic flow in the whole airspace network using the basic system models as described above. 4. Modeling of operation optimization cost function

47 48 49 50 51 52 53 54 55 56 57 58 59

As is known to all, one of the effective measures for air traffic system management is the total aircraft flight delay, which is calculated as a sum of delay times of all flights considered. The amount of delay mainly depends on how well the available air traffic system resource is utilized, especially during the congested periods. Therefore, total flight delay has been chosen as the optimality criterion. It is straightforward that the adjustment of aircraft arriving time at the jet route entrance corresponds to the parameter u inherently and the variation of aircraft flight speed corresponds to the parameters A, B and C . Supposing the number of aircraft involved in the conflict is m, we can present the following multi-aircraft operation optimization cost function:



60 61 62 63 64 65 66

87

¯ O B, ¯ all the sub-segments of the A¯  O B¯  are flight on the segment A considered as virtual resources and vice versa. We can also obtain a network that is analogous to Fig. 1 and the system model falls into the class of constraints introduced above. Taking two aircrafts ¯ O B¯ and called g and h as an example. Aircraft g is located on A ¯  O B¯  . If aircraft g arrives O firstly, this gives the h is located on A following constraint at the re-sequencing node O :

45 46

86

Fig. 2. Division of multiple climbing/descending jet routes.

21

min

α

n +m k=n+1

u (k) + β

n +m

  y (k)

(17)

k=n+1

α and β denote two weights for the system regulation methods mentioned above, α + β = 1 (α ≥ 0, β ≥ 0). Note that the parameters in the objective function (17) can have various meanings. For

88 89 90 91 92 93 94 95 96 97

A case study is presented where two airports system operation configurations for Shanghai terminal control area are modeled, analyzed and compared using the air traffic max-plus optimization model presented above. We consider the terminal airspace network surrounding Shanghai airports known as ZSPD and ZSSS in simulation as representative scenarios. Firstly, we consider the scheduled flights departing from ZSPD and ZSSS. The airspace networks are presented in Figs. 3–4 separately along with possible departure air routes configurations. For the two airspace networks, each node represents an air route characteristic location and the main location in the terminal control area is the corridor exit. The terminal control area starts with the aircraft taking off from the specified runway and then moving on each departure air route. All the aircrafts in the terminal control area follow the same flight procedure mentioned above but differ in departure air routes. In addition, we assume that the air traffic system networks previously described admit no routing and no internal overtaking. The initial flight information is presented in Table 3 over the 2-h period. During this period, five departure fixes are used for the outgoing air traffic flows, namely ‘ALDAP’, ‘PIKAS’, ‘SX’, ‘AND’ and ‘RUXIL’. In general, each aircraft can be placed into one of the fixes as indicated before. As shown in Table 3, all the flight conflicts can be classified into two categories: converging conflicts at ‘OLGAP’ and ‘NINAS’, along with following conflicts at ‘ALDAP’, ‘PIKAS’, ‘SX’, ‘AND’ and ‘RUXIL’. Specifically, aircraft flow between ‘SHA’ and ‘SX’ competes with ‘PUD’ to ‘SX’. Traffic flow between ‘SHA’ and ‘PIKAS’ competes with ‘PUD’ to ‘PIKAS’. In addition, the total-energy model and base of aircraft data are both applied for aircraft trajectory simulation. The total-energy model equates the rate of work done by forces acting on the aircraft to the rate of increase in potential and kinetic energy [14], that is:

( T H R − D ) V T A S = mg 0

dh dt

+ mV T A S

dV T A S dt

98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131

(18)

132

JID:AESCTE AID:3948 /FLA

[m5G; v1.208; Prn:16/03/2017; 10:43] P.7 (1-12)

Y.-x. Han et al. / Aerospace Science and Technology ••• (••••) •••–•••

7

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94

29

95

30

96

31

97

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40

106

41

107

42

108

43

109

Fig. 3. Simplified instrument departure chart of ZSPD.

44 45 46

110 111

Table 3 Initial flight information for ZSPD and ZSSS.

112

47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

113

Aircraft number

Aircraft type

Departure time

Departure airport

Destination airport

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20

B733 B738 B733 B763 A321 B737 B763 A319 B763 B763 B738 A320 A320 A321 A321 B738 A320 B763 B738 B738

1401 1402 1404 1405 1406 1408 1409 1410 1411 1411 1412 1414 1414 1416 1417 1419 1419 1421 1421 1422

ZSPD ZSPD ZSSS ZSPD ZSPD ZSSS ZSPD ZSPD ZSPD ZSSS ZSPD ZSSS ZSPD ZSPD ZSSS ZSSS ZSPD ZSPD ZSSS ZSSS

ZYDD ZHHH ZPPP WSSS ZGSZ ZBTJ KSFO ZYTL RJAA ZBTJ RKTU ZGHA ZUUU RKPC ZBAA ZGGG ZYHB RKSI ZWWW ZSFZ

Departure fix ALDAP PIKAS SX AND RUXIL PIKAS ALDAP PIKAS ALDAP PIKAS ALDAP SX PIKAS ALDAP PIKAS SX ALDAP ALDAP PIKAS SX (continued on next page)

114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

AID:3948 /FLA

[m5G; v1.208; Prn:16/03/2017; 10:43] P.8 (1-12)

Y.-x. Han et al. / Aerospace Science and Technology ••• (••••) •••–•••

8

1

67

2

68

3

69

4

70

5

71

6

72

7

73

8

74

9

75

10

76

11

77

12

78

13

79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94

29

95

30

96

31

97

32

98

33

99

34

100

35

101

Fig. 4. Simplified instrument departure chart of ZSSS.

36

102

37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

103 104

Table 3 (Continued) Aircraft number

Aircraft type

Departure time

Departure airport

Destination airport

21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50

A332 B738 B738 A321 A319 A320 B737 A320 B738 B733 B738 B738 B733 A321 B738 B733 B738 B772 B738 B738 B763 B763 B737 B738 B737 B738 A320 A320 A320 A346

1423 1424 1425 1427 1427 1429 1429 1431 1433 1434 1435 1437 1440 1443 1447 1448 1500 1504 1504 1505 1506 1507 1508 1509 1510 1511 1512 1514 1515 1517

ZSPD ZSSS ZSPD ZSSS ZSPD ZSPD ZSSS ZSPD ZSSS ZSPD ZSPD ZSSS ZSSS ZSPD ZSSS ZSSS ZSPD ZSSS ZSPD ZSSS ZSPD ZSSS ZSPD ZSPD ZSSS ZSSS ZSPD ZSPD ZSPD ZSPD

VMMC ZSQD ZLXY ZBSJ ZYTX RJGG ZGKL VTBS ZSFZ ZYCC RKPK ZSQD ZSWH VHHH ZBTJ ZSNB ZYTX ZBAA ZUUU ZGNN RJAA ZGGG RKSI VHHH ZGHA ZSJN ZYTL ZLXY VHHH KLAX

Departure fix AND PIKAS PIKAS PIKAS PIKAS ALDAP SX SX SX ALDAP ALDAP PIKAS PIKAS RUXIL PIKAS AND PIKAS PIKAS PIKAS SX ALDAP SX ALDAP RUXIL SX PIKAS PIKAS PIKAS RUXIL ALDAP (continued on next page)

105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:3948 /FLA

[m5G; v1.208; Prn:16/03/2017; 10:43] P.9 (1-12)

Y.-x. Han et al. / Aerospace Science and Technology ••• (••••) •••–•••

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23

9

67

Table 3 (Continued) Aircraft number

Aircraft type

Departure time

Departure airport

Destination airport

Departure fix

51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74

B763 A320 B763 B752 A320 A321 A320 A320 A320 A320 A319 A320 A321 B737 B738 A320 B763 B763 A320 B752 B752 A332 B738 B737

1518 1517 1519 1520 1521 1523 1527 1528 1529 1531 1531 1535 1535 1539 1542 1543 1545 1548 1551 1552 1552 1556 1557 1558

ZSPD ZSSS ZSSS ZSPD ZSPD ZSPD ZSSS ZSPD ZSPD ZSSS ZSPD ZSSS ZSPD ZSSS ZSSS ZSPD ZSSS ZSPD ZSPD ZSPD ZSSS ZSSS ZSPD ZSSS

RJAA ZGGG ZGSZ ZYHB ZGSD VHHH ZGGG ZYHB RKPC ZSWH ZUUU ZGHA ZYTL ZGSZ ZSAM ZYTX ZBAA RJAA ZYTX VMMC ZBAA ZGSZ ZHHH ZPPP

ALDAP SX SX ALDAP RUXIL RUXIL SX ALDAP ALDAP PIKAS PIKAS SX PIKAS SX SX PIKAS PIKAS ALDAP PIKAS AND PIKAS SX PIKAS SX

24

27

Symbols

Description

Units

29

THR

thrust acting parallel to the aircraft velocity vector aerodynamic drag aircraft mass geodetic altitude gravitational acceleration true airspeed time derivative

Newton

30 32 33 34 35

D m h g0 V T AS dV T A S /dt

Newton kilograms m m/s2 m/s s−1

36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66

70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 91

Table 4 Definition of variables in the total-energy model.

28

31

69

90

25 26

68

A summary of the parameters specified by the total-energy model is supplied in Table 4. Besides, the base of aircraft data provides a set of ASCII files containing performance and operating procedure coefficients for different aircraft types. The data can be used to calculate thrust, drag and fuel flow focusing on different flight profiles. The sample aircraft performance parameters are given in Table 5. Detailed information on how these parameters have been used for aircraft performance calculation is provided in [14]. Thus, we can use the ‘total energy model’ to generate aircraft trajectories consisting of horizontal position and vertical altitude over future time intervals. Trajectories are generated by integrating predefined waypoints, defined by the current aircraft position, filed flight plan route, and airspace description data. To simplify the optimization process, the horizontal and vertical components of the trajectory are decoupled. The horizontal trajectory is used as the basis for computing the vertical profile. Moreover, the candidate aircrafts are treated by similar statistical behavior to reduce the set of variables. This set defines what we call the “service activity” mentioned in section 2. We divide the whole flight trajectory into six sub-segments based on the aircraft altitude profile. Then we can obtain the system process time of six sub-segments. In the following section, we will mainly discuss the air traffic flow passing through ‘PIKAS’ and this can be extended easily to the other system fixes or the whole terminal control area. As is well known, the terminal control system comprises departure fixes and a runway system. There are two separate sets of departure fixes located in the near-terminal airspace area so that the departure fixes serve only departure flow. Furthermore, we can obtain the system resource occupy time or release time for each

aircraft in each departure fix so that we can detect the potential flight conflict of the terminal control system in advance. The flight operation at each characteristic location of ‘PIKAS’ and the corresponding required time (unit: s) for each aircraft are given in Table 6. The symbol ‘/’ implies that aircraft can’t receive the six serial services entirely because of the limited length of departure route. If the minimum requirement of safety separation time is 1 min, three scenarios are pursued and the system optimization results can be obtained. The specifics of the three cases are described below. Additionally, we denote the temperature deviation from international standard atmosphere by symbol ‘ISA+t’ in the following sections. (1) Case 1: The aircraft trajectories are optimized by modifying departure time instant with different temperature deviation such as ISA+0, ISA+10 and ISA-10. We denote the three temperature deviation scenarios by S11 , S12 and S13 respectively. (2) Case 2: The aircraft trajectories are optimized by modifying departure time instant with different climbing mode such as the reduced climbing power mode is employed or not. We denote the two climbing power modes by S21 and S22 respectively. (3) Case 3: The aircraft trajectories are optimized by modifying aircraft departure time instant and departure speed schedule. We denote the two adjustment patterns by symbols ‘DT’ and ‘DS’ respectively. By solving the system optimization model (17), we analyzed six typical cases without considering buffer and the corresponding results are shown in Table 7. It can be seen from Table 7 that the optimization results differing for different optimization models. Clearly, optimization model (6) performs best compared to the other models, but this increases the control complexity for air traffic controllers. For the typical cases with buffers, since the airspace network configurations are similar, detailed analysis is required to compare and choose the best among them. Using equations (5)–(8) and the planned take off time of each aircraft, the exact starting process time of every characteristic location in the airspace network can also be obtained. Then for the given system processing times, following the procedure in section 2, the max-plus system optimization results with buffers size of one or two are given in Table 8 assuming that all the buffers have the same size. The in-

92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

AID:3948 /FLA

[m5G; v1.208; Prn:16/03/2017; 10:43] P.10 (1-12)

Y.-x. Han et al. / Aerospace Science and Technology ••• (••••) •••–•••

10

1 2 3 4 5 6

67

Table 5 Summary of simplified performance parameters of A320.

68

Symbols

Description

Value

Symbols

Description

Value

69

C Tc,1 C Tc,2 C Tc,3

1st max-climb coefficient 2nd max-climb coefficient 3rd max-climb coefficient

1.41E+05 4.89E+04 6.50E−11

C f cr hMO hmax

0.95423E+00 3.90E+04 3.40E+04

70

C Tc,4 S C D0,CR C D2,CR C f1

1st temperature coefficient reference wing surface area parasitic drag coefficient induced drag coefficient 1st fuel coefficient

9.98E+00 1.23E+02 2.51E−02 3.61E−02 6.33E−01

Gt Gw mmax V stall_to C f2

cruise fuel coefficient max-operation altitude max-altitude at MTOW and ISA temperature gradient weight gradient maximum mass take-off stall speed 2nd fuel coefficient

7 8 9 10 11

−7.30E+01 2.65E−01 7.70E+01 1.14E+02 8.59E+02

12 13 14 15 16

19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40

43

46 47 48

51 52 53 54 55 56

Characteristic location Cl1

Cl2

Cl3

Cl4

Cl5

Cl6

59 60 61 62 63 64 65 66

76 77

2 6 8 10 13 15 19 22 23 24 25 32 33 35 37 38 39 46 47 48 60 61 63 66 67 69 71 73

43 41 43 45 45 44 43 43 43 44 43 43 41 43 43 52 43 43 45 45 45 43 44 45 45 45 45 43

55 52 54 58 58 57 55 55 55 57 54 55 53 55 55 64 55 55 58 58 58 54 57 58 58 58 58 55

49 48 49 50 50 50 49 49 49 50 49 49 48 49 49 53 49 49 50 50 50 49 50 50 50 50 50 49

38 36 37 39 39 39 38 38 38 39 37 38 36 38 38 41 38 38 39 39 39 37 39 39 39 39 39 38

35 42 50 28 39 34 35 35 35 34 50 35 34 35 35 38 35 35 39 39 39 50 34 39 28 39 20 35

124 / 131 / 138 / / / 124 / 131 / / / 124 / 124 / 138 138 / 131 139 138 / 138 / 124

82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108

Table 7 Trajectory optimization types and results without buffers.

109 110

Optimization model

Optimization result/s

Optimization model

Optimization result/s

(1) M1:DT ∨ S11 ∨ S21 (2) M2:DT ∨ S11 ∨ S22 (3) M3:DT ∨ S12 ∨ S21

600 660 720

(4) M4:DT ∨ S12 ∨ S22 (5) M5:DT ∨ S13 ∨ S21 (6) M6:DT ∨ S11 ∨ S21 ∨ DS

720 660 420

111 112 113 114 115

Table 8 Trajectory optimization types and results with one/two buffers. Optimization model

Optimization result/s

Optimization model

Optimization result/s

(1) M1 (2) M2 (3) M3

540/480 600/540 660/600

(4) M4 (5) M5 (6) M6

600/540 540/540 360/300

57 58

75

81

Aircraft number

49 50

74

80

44 45

73

79

Table 6 Flight operation and required time for each characteristic location in ‘PIKAS’.

41 42

72

78

17 18

71

troduction of buffer size can improve optimization results to some extent for the air traffic system. Moreover, the total fuel consumptions of all the flight using different optimization models are given in Fig. 5. Consequently, the aircraft fuel consumption can be easily calculated for different operation configurations, which can provide effective reference for air traffic system optimization regulation. This is very useful, especially in the airspace system design stages when the exact flight

time for the given air traffic flow on a specified characteristic location is unknown and the system designers want to know the effect of variation in flight time on the fuel efficiency of the air route. In a similar way, we can also obtain the fuel consumption of aircraft with distinct buffer size configurations considered above. Table 9 shows the amount of overall fuel consumption of all the flight corresponding to ‘PIKAS’ and it can be seen from Table 9 that the amount increases with considering buffer size. Furthermore, the same conclusions can be achieved when increasing the buffer size to two. As depicted in Table 9, the fuel consumption increases with increasing the size of buffers. Another useful application of Table 9 is in evaluating the effect of changing the flight times on a given performance measure. Additionally, it can also be used during the air traffic system operation phase to assess the merits and trade-off of installing new navigation equipment or increasing new runways which will improve the fuel efficiency of air traffic system. In this

116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE AID:3948 /FLA

[m5G; v1.208; Prn:16/03/2017; 10:43] P.11 (1-12)

Y.-x. Han et al. / Aerospace Science and Technology ••• (••••) •••–•••

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

Fig. 5. Plots illustrating the fuel consumption of each aircraft. (For interpretation of the colors in this figure, the reader is referred to the web version of this article.)

22 23 24

Table 9 Total fuel consumption of aircraft without & with buffers.

25

Model type

Without buffer/kg

With one buffer/kg

26

M1 M2 M3 M4 M5 M6

2.2860e+04 2.2164e+04 2.3091e+04 2.2354e+04 2.2728e+04 2.3537e+04

2.3587e+04 2.2783e+04 2.3692e+04 2.2983e+04 2.3462e+04 2.4158e+04

27 28 29 30 31 32 33 34

way, the airspace designers can easily determine the system optimal processing time that minimizes the total fuel cost.

35 36

6. Conclusions

37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55

Described in this paper is the new approach to model deterministic air traffic systems taking into consideration finite buffers. This approach extends the class of systems that can be analyzed via max-plus theory. The method is based on the observation that the jet route can be decomposed into different sub-segments with the same features. The proposed max-plus linear model of air traffic system also has explicit buffer control mechanism and the underlying max-plus model describing air traffic system can provide means to calculate performance characteristics for the system. The modeling approach can be used in both air traffic system design and operation phases. Moreover, the approach is not limited to deterministic air traffic systems and it can be extended to stochastic scenarios with minor modifications. However, the air traffic system max-plus linear model put forward in this paper is not straightforward for the system to enter into uniform periodic steady state which can enhance the balance of system operation and we have to prove the way in which we can improve the uniformity of air traffic system combining the elementary analysis results in future.

56 57

Conflict of interest statement

58 59

None declared.

60 61

References

62 63 64 65 66

[1] A. Alonso-Ayuso, L.F. Escudero, P. Olaso, C. Pizarro, Conflict avoidance: 0-1 linear models for conflict detection & resolution, Top 21 (3) (2013) 1–20. [2] A. Alonso-Ayuso, L.F. Escudero, F.J. Martín-Campo, Multiobjective optimization for aircraft conflict resolution. A metaheuristic approach, Eur. J. Oper. Res. 248 (2) (2016) 691–702.

11

[3] A. Seleim, H. Elmaraghy, Generating max-plus equations for efficient analysis of manufacturing flow lines, J. Manuf. Syst. 37 (2014) 426–436. [4] A.E. Vela, S. Solak, J.B. Clarke, W.E. Singhose, E.R. Barnes, E.L. Johnson, Near real-time fuel-optimal en route conflict resolution, IEEE Trans. Intell. Transp. Syst. 11 (4) (2010) 826–837. [5] A.K. Agogino, K. Tumer, A multiagent approach to managing air traffic flow, Auton. Agents Multi-Agent Syst. 24 (1) (2012) 1–25. [6] B. Heidergott, A characterisation of (max, +)-linear queueing systems, Queueing Syst. 35 (1) (2000) 237–262. [7] C. Peyronne, A.R. Conn, M. Mongeau, D. Delahaye, Solving air traffic conflict problems via local continuous optimization, Eur. J. Oper. Res. 241 (2) (2015) 502–512. [8] C. Tomlin, I. Mitchell, R. Ghosh, Safety verification of conflict resolution maneuvers, IEEE Trans. Intell. Transp. Syst. 2 (2) (2001) 110–120. [9] C.A. Nicolas, Barnier, 4D-trajectory de-confliction through departure time adjustment, in: 8th USA/Europe Air Traffic Management Research and Development Seminar, 2009. [10] C.A. Nicolas Barnier, Combining flight level allocation with ground holding to optimize 4D-deconfliction, in: 9th USA/Europe Air Traffic Management Research and Development Seminar, 2011. [11] C.J. Tomlin, J. Lygeros, S.S. Sastry, A game theoretic approach to controller design for hybrid systems, Proc. IEEE 88 (7) (2000) 949–970. [12] D. Sislak, P. Volf, M. Pechoucek, N. Suri, Automated conflict resolution utilizing probability collectives optimizer, IEEE Trans. Syst. Man Cybern., Part C, Appl. Rev. 41 (3) (2011) 365–375. [13] D.Z. Zheng, Q.C. Zhao, Discrete Event Dynamic System, Tsinghua University Press, Beijing, 2000, pp. 96–173. [14] E.E.C. Eurocontrol Experimental Centre, User Manual for the Base of Aircraft Data (BADA), Revision 3.8, Eurocontrol, Montreal, 2010. [15] G. Cohen, S. Eacute, P. Gaubert, J.P. Quadrat, Max-plus algebra and system theory: where we are and where to go now, Annu. Rev. Control 23 (1) (1999) 207–219. [16] G.J. Olsder, Max algebra approach to discrete event systems, in: IEE Colloquium on Discrete Event Systems: A New Challenge for Intelligent Control Systems, 1993. [17] H. Liang, H. Yang, B. Yang, Medium term conflict detection model based on approximate optimal solution, J. Comput. Inf. Syst. 8 (18) (2012) 7815–7821. [18] Í.R.D. Oliveira, F.S. Carvalho, J.B.C. Junior, L.M. Sato, Multi-agent tools for air traffic management, IEEE Int. Conf. Comput. Sci. Eng. Workshops (2008). [19] J. Omer, A space-discretized mixed-integer linear model for air-conflict resolution with speed and heading maneuvers, Comput. Oper. Res. 58 (2015) 75–86. [20] J. Tang, M.A. Piera, T. Guasch, Coloured Petri net-based traffic collision avoidance system encounter model for the analysis of potential induced collisions, Transp. Res., Part C, Emerg. Technol. 67 (2016) 357–377. [21] J.A. Besada, G. Frontera, J. Crespo, E. Casado, J. Lopez-Leones, Automated aircraft trajectory prediction based on formal intent-related language processing, IEEE Trans. Intell. Transp. Syst. 14 (3) (2013) 1067–1082. [22] J.H. Hu, M. Prandini, S. Sastry, Optimal coordinated maneuvers for threedimensional aircraft conflict resolution, J. Guid. Control Dyn. 25 (5) (2015) 888–900. [23] J.K. Kuchar, L.C. Yang, A review of conflict detection and resolution modeling methods, IEEE Trans. Intell. Transp. Syst. 1 (4) (2000) 179–189. [24] K. Tumer, A. Agogino, Improving air traffic management with a learning multiagent system, IEEE Intell. Syst. 24 (1) (2009) 18–21. [25] M. Nguyen-Duc, J.P. Briot, A. Drogoul, V. Duong, An application of multi-agent coordination techniques in air traffic management, in: IEEE/WIC International Conference on Intelligent Agent Technology, 2003. [26] N. Dougui, D. Delahaye, S. Puechmorel, M. Mongeau, A light-propagation model for aircraft trajectory planning, J. Glob. Optim. 56 (3) (2013) 873–895. [27] N. Durand, J.M. Alliot, J. Noailles, A.E. Belin, A.E.B.R. Camichel, Automatic aircraft conflict resolution using genetic algorithms, in: Proceedings of the ACM Symposium on Applied Computing, ACM, Association for Computing Machinery, 1996. [28] P. Menon, G. Sweriduk, K. Bilimoria, New approach for modeling, analysis, and control of air traffic flow, J. Guid. Control Dyn. 27 (5) (2004) 737–744. [29] R. Goverde, P.H.L. Bovy, G.J. Olsder, The max-plus algebra approach to transportation problems, in: World Transport Research: Selected Proceedings of the 8th World Conference on Transport Research, 1999. [30] S. Cafieri, N. Durand, Aircraft deconfliction with speed regulation: new models from mixed-integer optimization, J. Glob. Optim. 58 (4) (2014) 613–629. [31] S. Harry, B. Richard, L. Michael, Next Generation Air Transportation System (NGATS) Air Traffic Management (ATM)-Airspace Project, NASA, Washington, 2006. [32] S. Ruiz, M.A. Piera, J. Nosedal, A. Ranieri, Strategic de-confliction in the presence of a large number of 4D trajectories using a causal modeling approach, Transp. Res., Part C, Emerg. Technol. 39 (2014) 129–147. [33] S. Ruiz, M.A. Piera, I.D. Pozo, A medium term conflict detection and resolution system for terminal maneuvering area based on spatial data structures and 4D trajectories, Transp. Res., Part C, Emerg. Technol. 26 (1) (2013) 396–417.

67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132

JID:AESCTE

12

1 2 3 4 5 6 7 8 9 10 11 12 13

AID:3948 /FLA

[m5G; v1.208; Prn:16/03/2017; 10:43] P.12 (1-12)

Y.-x. Han et al. / Aerospace Science and Technology ••• (••••) •••–•••

[34] S. Wollkind, J. Valasek, T. Ioerger, Automated conflict resolution for air traffic management using cooperative multiagent negotiation, in: AIAA Guidance, Navigation, & Control Conference, 2013. [35] T. Tarnopolskaya, N. Fulton, H. Maurer, Synthesis of optimal bang–bang control for cooperative collision avoidance for aircraft (ships) with unequal linear speeds, J. Optim. Theory Appl. 155 (1) (2012) 115–144. [36] V.P. Jilkov, X.R. Li, J.H. Ledet, An efficient algorithm for aircraft conflict detection and resolution using List Viterbi Algorithm, in: International Conference on Information Fusion, 2015. [37] W. Chao, G. Jing, X. Xu, Analysis of air traffic flow control through agent-based modeling and simulation, in: International Conference on Computer Modeling and Simulation, 2009. [38] W. Schuster, W. Ochieng, Performance requirements of future trajectory prediction and conflict detection and resolution tools within SESAR and next-gen: framework for the derivation and discussion, J. Air Transp. Manag. 35 (2014) 92–101.

[39] X.B. Hu, S.F. Wu, J. Jiang, On-line free-flight path optimization based on improved genetic algorithms, Eng. Appl. Artif. Intell. 17 (8) (2004) 897–907.

67

[40] X.J. Zhang, X.M. Guan, H. Inseok, K.Q. Cai, A hybrid distributed-centralized conflict resolution approach for multi-aircraft based on cooperative coevolutionary, Sci. China Inf. Sci. 56 (12) (2013) 1–16.

69

[41] X.M. Tang, P. Chen, B. Li, Optimal air route flight conflict resolution based on receding horizon control, Aerosp. Sci. Technol. 50 (2016) 77–87. [42] Y. Matsuno, T. Tsuchiya, J. Wei, I. Hwang, N. Matayoshi, Stochastic optimal control for aircraft conflict resolution under wind uncertainty, Aerosp. Sci. Technol. 43 (2015) 77–88. [43] Y. Xu, H.H. Zhang, Z.H. Liao, L. Yang, A dynamic air traffic model for analyzing relationship patterns of traffic flow parameters in terminal airspace, Aerosp. Sci. Technol. 55 (2016) 10–23.

68 70 71 72 73 74 75 76 77 78 79

14

80

15

81

16

82

17

83

18

84

19

85

20

86

21

87

22

88

23

89

24

90

25

91

26

92

27

93

28

94

29

95

30

96

31

97

32

98

33

99

34

100

35

101

36

102

37

103

38

104

39

105

40

106

41

107

42

108

43

109

44

110

45

111

46

112

47

113

48

114

49

115

50

116

51

117

52

118

53

119

54

120

55

121

56

122

57

123

58

124

59

125

60

126

61

127

62

128

63

129

64

130

65

131

66

132