A rolling horizon optimization approach for dynamic airspace sectorization

IFAC Journal of Systems and Control 11 (2020) 100076

Contents lists available at ScienceDirect

IFAC Journal of Systems and Control journal homepage: www.elsevier.com/locate/ifacsc

A rolling horizon optimization approach for dynamic airspace sectorization ∗

C.S.Y. Wong a , , S. Suresh b , N. Sundararajan a a b

SCSE, Nanyang Technological University, Singapore 639798, Singapore Department of Aerospace Engineering, Indian Institute of Science, Bangalore, India

article

info

Article history: Received 1 January 2019 Received in revised form 12 May 2019 Accepted 31 January 2020 Available online 7 February 2020 Keywords: Dynamic sectorization Rolling horizon optimization Airspace Sector similarity

a b s t r a c t Dynamic airspace sectorization is the concept of dynamically changing the airspaces that the air traffic controllers are managing to handle sudden changes in the air traffic patterns and also the nonavailability of certain airspaces due to weather or other considerations. For the dynamically changing sector shapes to be operationally feasible, air traffic controllers need to be familiar with the sector shapes i.e. the new sector shapes need to have a high level of similarity with the previous sector shapes. This paper presents a rolling horizon optimization approach to plan and handle the above problem. Normally in a conventional dynamic airspace sectorization optimization approach, one changes the sector shapes for the next hour interval based on the current sector shapes and the planned next hour traffic pattern. However, in this paper, we propose a rolling horizon optimization approach where the shapes are changed based on a finite horizon of the traffic in future time intervals. Detailed performance evaluation of the proposed method is presented based on Singapore flight information region traffic and historical data. The results clearly show that the proposed approach is better able to provide feasible gradually changing sector shapes that are more suitable for operations when compared to the conventional approach. © 2020 Elsevier Ltd. All rights reserved.

1. Introduction Recent forecast of air traffic growth for the next twenty years have indicated an increasing trend (Federal Aviation Adminstration (F.A.A.), 2017) and this growth necessitates the need for new approaches in managing the air traffic safely and efficiently (Wang, Song, Wen, & Zhao, 2016; Xu, Zhang, Liao, & Yang, 2016). One of the major problems in the management of airspace is on the sectorization of the given airspace so controllers can manage them efficiently (Delahaye & Puechmorel, 2008; Kicinger & Yousefi, 2009; Sergeeva, Delahaye, & Mancel, 2015; Tang, Alam, Lokan, & Abbass, 2012a; Wong, Venugopalan, & Suresh, 2016; Xue, 2010). Given a set of traffic flows, a set of optimal sectors would ensure the Air Traffic Controllers (ATCs) workloads such that the ATCs are able to handle them comfortably and that there is no unnecessary coordination required between ATCs in handling different sectors. However, given the varying characteristics of the traffic flows, there may not be a single set of optimal sectors that will fulfil the desired objectives over all the time. Therefore, the idea of dynamic sectorization (Chen, Bi, Zhang, & Song, 2013; Schultz, Gerdes, Standfuß, & Temme, 2017) – where ∗ Corresponding author. E-mail addresses: [email protected] (C.S.Y. Wong), [email protected] (S. Suresh), [email protected] (N. Sundararajan). https://doi.org/10.1016/j.ifacsc.2020.100076 2468-6018/© 2020 Elsevier Ltd. All rights reserved.

the shapes and even the number of sectors are changed over time according to the traffic flows, is relevant. However, many argue that the idea of dynamic sectorization is not operationally feasible as it takes time for a traffic controller to adapt to the new shapes of the sectors. Following this, Tang, Alam, Lokan, and Abbass (2012b) proposed a sector similarity metric as an objective to be optimized in a multi-objective optimization framework to prevent drastic changes in sector shapes. Alternatively, Sergeeva (Sergeeva, 2017; Sergeeva et al., 2015) introduced the concept of functional airspace blocks, which could be separated or combined according to the changing traffic flows. This method trade-offs the flexibility of the sector shapes for a higher operational feasibility by building in a structure to the sector shapes. Recently, Gerdes, Temme, and Schultz (2018) addressed the issue of a smooth transition between two sectorizations for different air traffic demand. They introduced a method of changing the vertices of sector shapes with an evolutionary algorithm and interim diagrams to help with gradual sector shape changes such that the they can be operationally feasible. The need for similar sector shapes in Dynamic Airspace Sectorization (DAS) requires us to closely look at the problem with future traffic flows during the selection of current sector shapes. A Cognitive decision-making architecture for Dynamic Airspace Sectorization (CDAS) (Wong, Suresh, & Sundararajan, 2019) has been proposed for the user to consider the next hour traffic flows

2

C.S.Y. Wong, S. Suresh and N. Sundararajan / IFAC Journal of Systems and Control 11 (2020) 100076

when selecting the sector shapes for the current hour. However, CDAS is unable to guarantee solutions with similar sector shapes in the following time intervals. Hence, in this paper, we consider the future traffic flows over a finite number of time intervals for selecting the next hour sector shapes in the DAS problem. This is referred to the rolling horizon approach for determining sector shapes. The rolling horizon framework is commonly used for decisionmaking in dynamic stochastic problems (Sethi & Sorger, 1991). One such example would be the distributed control of robotic networks (Bullo, Cortés, & Martínez, 2009) where the position of robots are always changing. It is also used in path planning and optimizing trajectories (Sharma & York, 2018; Tang, Chen, & Li, 2016; Zhang, Li, & Wang, 2018) in aerospace applications, where control and state trajectories for a dynamic system are determined over a period of time based on the minimization of a performance index. In these problems, the trajectories can be mathematically represented by a system model derived from the physics considerations and optimal solutions are obtained by conventional methods such as dynamic programming and convex optimization. However, the dynamic airspace sectorization problem dealt in this paper cannot be fully represented by a mathematical model and hence optimization over finite time horizon is carried out with the help of a genetic algorithm, a solver meant for more general problems (Michalewicz, Krawczyk, Kazemi, & Janikow, 1990). The DAS problem is casted into a rolling horizon optimization framework. The current state in the framework is represented by current traffic flows and future states by predicted traffic flows. The control action is the choice of sector shapes for the next time interval. The optimal sector shapes for the next time interval is the one that minimizes the standard deviation of dynamic density-based (Laudeman & Center., 1998) monitoring workloads among sectors for a given time horizon. Moreover, the sector shapes has to fulfil all the required spatial, safety and sector similarity (Tang et al., 2012b) constraints. The dynamic densitybased monitoring workload measurement takes in consideration the traffic complexity in the airspace, such as changes in altitude of aircraft and the presence of potential conflicts. Details of the dynamic density measurement metric are given later. The spatial constraint restricts the generation of Voronoi sites, which are used to generate sector shapes, to the airspace of interest. The safety constraints consist of a coordination workload limit and that any crossing points between flights do not occur near the sector boundaries. This ensures that the ATC would be given sufficient time to execute required actions and will not be overwhelmed. The sector similarity constraint is an inter-time period constraint that compares the current sector shapes with the sector shapes in the previous time period for similarity in shapes. With the sector similarity constraint, one can ensure that sector shapes are operationally feasible i.e. does not vary too much from the previous sector shapes. The optimization is carried out with the use of a genetic algorithm. The performance of the proposed rolling horizon optimization approach has been evaluated using historical flight data in the Singapore Flight Information Region (FIR) and has been compared with the conventional approach of a single stage optimization results. The results clearly indicate that the rolling horizon optimization approach is better suited to handle future sudden traffic variations, meeting all the constraints and providing gradual sector shape changes for operational feasibility. The paper is organized as follows. First, Section 2 introduces the basic airspace sectorization concepts used in the proposed DAS model. Section 3 presents the DAS problem in a rolling horizon optimization framework, with details on performance index, constraints and chromosome representation in the genetic algorithm. Section 4 shows the results of the performance evaluation of the proposed method. Section 5 summarizes the conclusions from this study.

2. Airspace sectorization concepts Before presenting the rolling horizon framework, some basic elements of Air Traffic Control namely airspace sectorization concepts, definition of workload and similarity metrics are described in this section. 2.1. Airspace sector generation In this paper, the methodology used for the generation of 2D sector shapes is by Voronoi Diagrams (Fortune, 0000). Voronoi diagram is a partitioning technique of subdividing a given region into sub-regions. A set of points known as sites are used to generate these partitions satisfying the condition that, any point within each partition is closest to its corresponding site. Such diagrams have an interesting property that the resulting shape of the partitions are strictly convex (internal angles are lower than 180 deg). As convex sector shapes are preferred by air traffic controllers, this is a widely used technique in the automatic sectorization literature (Sergeeva et al., 2015; Tang et al., 2012b; Xue, 2010). The set of sites used to generate the sector shapes are denoted as V = v1 , . . . , vK where K is the total number of 2D sector shapes. vi =

( λ) vi φ

vi

,

i = 1, . . . , K

(1)

Superscript λ and φ denote the longitude and latitude of the site’s location respectively. The sector boundaries are determined through the Voronoi decomposition of V . The use of a Voronoi diagram is effective when the given airspace is of a cuboid shape. However, in real world, the largest regular divisions of the airspace — Flight Information Regions (FIRs) are usually determined territorially by different countries covering that airspace. The airspace owned by these countries are mostly made up of non-convex constrained spaces. The issue of non-convex airspace boundaries has been addressed by Gerdes et al. (2018) with an appropriate data structure called Doubly Connected Edge List (DCEL) (Berg, Cheong, Kreveld, & Overmars, 2008). The proposed approach in Gerdes et al. (2018) identifies the true boundaries of every sector, which can be used to find out the boundaries of each sector in the final solution given a nonconvex airspace. For the purpose of our experiment, the Voronoi sites are to be generated in the constrained airspaces to ensure that every sector shape is within the FIR and hence the provided solutions are feasible. To obtain feasible 3D sectors, the shapes of an airspace sector should fulfil the condition of a right prism for the ease of handling by controllers. Voronoi diagrams are first used to generate the shapes of airspace sectors in 2D (longitude and latitude). This is followed by the generation of 3D sectors through a vertical splitting of the given area at a certain altitude when required. 2.2. Definitions of airspace workload and similarity metrics A brief explanation of the following definitions for the metrics used in monitoring workloads, coordination workloads and sector similarities in the air traffic control area are given below. 2.2.1. Definition of dynamic density used for monitoring workloads Dynamic density is an air traffic management metric that was proposed to measure both the traffic density and complexity together (Laudeman & Center., 1998). The dynamic density, (ρ ) for a given sector is defined as:

ρ=

9 ∑ i=1

Wi (θi ),

(2)

C.S.Y. Wong, S. Suresh and N. Sundararajan / IFAC Journal of Systems and Control 11 (2020) 100076

where θi represents the traffic complexity factors and Wi represents the weights assigned to each factor. Here, the W is set as [2.40, 2.45, 2.94, 2.45, 1.83, 4, 3, 2.11, 1]. These weight values are relative to the traffic density (θTD ), which is given a weight of 1. The details of the study that was conducted to determine the weights can be found in Laudeman and Center. (1998). The definitions of traffic complexity factors θi and traffic density θTD are given by: 1. Heading Change (θHC ) - The number of aircraft that made a heading change of more than 15 degrees during a sample interval of two minutes 2. Speed Change (θSC ) - The number of aircraft that had a computed airspeed change of greater than 10 knots or 0.02 Mach during a sample interval of two minutes 3. Altitude Change (θAC ) - The number of aircraft that made an altitude change of greater than 750 ft during a sample interval of two minutes. 4. Minimum Distance 0–5 n. mi (θMD5 ) - The number of aircraft that had a Euclidean distance of 0–5 n. mi. to the nearest aircraft at the end of each two minute sample interval, excluding conflicting aircraft. 5. Minimum Distance 5–10 n. mi (θMD10 ) - The number of aircraft that had a Euclidean distance of 5–10 n. mi. to the nearest aircraft at the end of each two minute sample interval, excluding conflicting aircraft. 6. Conflict Predicted 0–25 n. mi. (θCP25 ) - The number of aircraft predicted to be in conflict with another aircraft whose lateral distance at the end of each two minute sample interval was 0–25 n mi. (with a vertical separation of less than 1000 ft) 7. Conflict Predicted 25–40 n. mi. (θCP40 ) - The number of aircraft predicted to be in conflict with another aircraft whose lateral distance at the end of each two minute sample interval was 25–40 n mi. (with a vertical separation of less than 1000 ft) 8. Conflict Predicted 40–70 n. mi. (θCP70 ) - The number of aircraft predicted to be in conflict with another aircraft whose lateral distance at the end of each two minute sample interval was 40–70 n mi. (with a vertical separation of less than 1000 ft) 9. Traffic Density (θTD ) - The sum of the number of aircraft in a particular area over time 2.2.2. Definition of coordination workloads The coordination workload is the total number of flights leaving a sector during the given time period and is represented by cj = #aircrafts leav ing sector j

(3)

2.2.3. Definition of sector similarity metric Sector similarity is a metric used to measure the changes in sector shapes (Tang et al., 2012b). It measures the overlapping areas of two sector shapes from two different airspace sectorization. First, the maximum overlapping area ak of Sk,d+1 with all the sectors from previous sector shapes is computed K

ak = max (Sk,d+1 ∩ Sj,d ) j=1

k ∈ [1, K ],

(4)

where Sj,d is one of the sectors from the previous airspace sectorization (d) and Sk,d+1 is a sector from the current airspace sectorization (d + 1). This overlap is computed using the Weiler– Atherton Clipping algorithm (Weiler & Atherton, 1977). The maximum overlap area does not reflect the similarity effectively due to varying sector sizes. Hence, a normalized sector similarity (rk ) of

3

sector Sk is computed as the ratio between the maximum overlap area (ak ) and the area of previous sector (Sj,d ). rk =

ak Sj,d

,

k = 1, . . . , K

(5)

where Sj,d is the original sector which has the biggest overlap with the current sector Sk . A higher value of rk means the sectors are more similar. For sector similarity to be an effective metric, a one-to-one parent–child relationship between the old and new sectors has to be established. This is being assigned using an integer linear-programming model with the constraint that each old sector is only mapped to one new sector, and vice-versa. More details on this integer linear-programming model are given in Venugopalan, Wong, Suresh, and Sundararajan (2018). 3. Dynamic airspace sectorization in a rolling horizon framework A rolling horizon approach allows optimization of the systems at the current time while taking into account the planned traffic flows in future based on a finite horizon, but only implementing the solution for the current (first) interval at the current time. This process is repeated at the next time instant with planned traffic flows in the next time horizon. Hence, the rolling horizon approach has the ability to take appropriate control actions now while keeping the future in mind. In the DAS problem, the states are the traffic flows in the airspace and the control actions are the determination of proper sector shapes for the airspace. The objective function for the optimization is to keep the monitoring workloads for the sectors balanced. The ultimate objective of the airspace sectorization problem is to provide an airspace sectorization for each time period such that the monitoring workload between the controller are balanced when the traffic flows change, meeting the required safety and coordination workload constraints and also ensuring minimum changes in the shapes of the sectors between these time periods. The use of rolling horizon in DAS allows one to ensure that the shapes of the sectors meet a user-defined sector similarity in the given finite-time horizon, which otherwise would not be captured when optimizing only for a single time period at a time. Being a non-linear optimization problem which cannot be fully represented in mathematical models, the actual optimization tool used in this work is a genetic algorithm. 3.1. Overall structure Fig. 1 gives an overview of the rolling horizon framework for the DAS problem. At the current time period t, the sectors used in t + 1 to t + p are being considered (where the finite horizon is p). The sectors in t + 1 to t + p periods are being optimized with the objective of balancing the ATC monitoring workload M over all the periods. The optimal sectors in each time period is subjected to the spatial constraints of the given airspace, safety constraints and the maximum coordination workloads of the sectors C = [c1 , c2 , . . . , cN ] where N is the total number of sectors. The number of sectors is being kept unchanged throughout such that the shapes of sectors can be similar and comparable. Between the current and next time period, the optimization of sectors is subjected to a lower bound on the sector similarity metric. A normalized standard deviation of monitoring workload (M t +1 , M t +2 , . . .) is calculated for each period based on the sectors in the solution set. With that, the mean M over p + 1 periods is utilized as the performance index for minimization.

4

C.S.Y. Wong, S. Suresh and N. Sundararajan / IFAC Journal of Systems and Control 11 (2020) 100076

Fig. 1. Rolling horizon framework for dynamic airspace sectorization.

3.2. Details on the optimization model formulation 3.2.1. Objective: Minimizing the standard deviation of monitoring workload across different time periods For a given traffic pattern in sector j, the normalized monitor√ ρj2 ∑N

ing workload (Mj =

j=1

ρj2

) is computed using the dynamic

density measure (ρ ) (Laudeman & Center., 1998) given in Eq. (2). This normalization allows for a fair comparison of monitoring workloads across different time periods. The objective for each time period, t, is to balance the monitoring workloads across the controllers in different sectors. This is mathematically expressed as

   t M = min(√ where µm =

N ∑

1 N −1 1 N

(Mj − µm )2 ),

(6)

j=1

∑N

j=1 Mj .

The main objective of the optimization is to minimize M, the mean standard deviation of monitoring workloads over p time periods and is formally defined as

∑t +p M=

t =t +1

Mt

p

(7)

Next, a brief description of the constraints are given. 3.2.2. Intra time period constraints 1: spatial constraints The Voronoi sites vi of all K 2D sector shapes must exist within the constrained airspace S. More details on sector generation using Voronoi sites can be found in Section 2.1. This can be mathematically presented as C 1 : vi ∈ S

∀i = 1, . . . , K

(8)

3.2.3. Intra time period constraint 2 (safety constraint): minimum distance between the crossing points (within a sector) to the boundaries of that sector For the safe operation of aircraft and avoidance of potential conflicts between two aircraft in a sector, ATC requires a minimum distance between the sector boundary and the crossing traffic flow. The crossing point (CP) is identified as an area prone to conflicts due to overlapping flows. CP is defined based on a wider range of points where aircraft are within lateral distances of 0–10 nautical miles and vertical separations of lower than 1000 ft. The minimum Euclidean distance (σj ) between the crossing points (CP) and the sector boundaries is defined as

σj = min ∥CP − Bjk ∥2 , k=1,...,r

j = 1, 2, . . . , N ,

(9)

where (CP) is a vector representing the crossing points and Bjk stores vertex k of the boundaries of sector j. In this study, the minimum distance has been set as 10 nautical miles (NM). Based on this, the following safety constraints are imposed: C 2 : σj > 10NM ,

∀j = 1, 2, . . . , N

(10)

3.2.4. Intra time period constraint 3: maximum sector coordination workload This constraint assigns a limit (clim ) on the maximum coordination workload of each sector j, ensuring that there is no overload of coordination workloads for controllers. More details on the coordination workload metric can be found in Section 2.2.2. C 3 : cj ≤ clim ,

∀j = 1, 2, . . . , N

(11)

C.S.Y. Wong, S. Suresh and N. Sundararajan / IFAC Journal of Systems and Control 11 (2020) 100076

5

Fig. 2. Input chromosome for N sectors and K areas for a single time interval.

3.2.5. Inter time period constraint: sector similarity For operational purposes, the shapes of previous sectors seen by the controllers should be similar to the current shapes of the sectors. Details on the sector similarity can be found in Section 2.2.3. The sector similarity SS t for hour t is the average sector similarity (measure) of all k sector shapes and is given by

s2 = q2 − kc ∗ βc ∗ xc ∗ (q2 − q1 )

(15)

Note that the sector similarity is measured in a 2D plane. To satisfy the right prism constraints, the sector similarity of 3D sectors are calculated based on the 2D mapping of the respective 3D sector shapes. The sector similarity SS t must be greater than a user-defined limit, SSlim . The constraint is provided as follows

s1 and s2 are the new child chromosomes produced by the parent chromosomes q1 and q2 in the population. The crossover operation is ran 0.5 ∗ Np times, giving a new population of Np . kc is a vector of binary values, that reflects 1 if a randomly generated value exceeds the crossover fraction of 2/nVars where nVars represents the number of variables. βc is a user-defined value that is typically set between [0,2], solutions tend to differ more when set above 1. xc is a randomly generated vector with values between 0 and 1. The Gaussian mutation process which follows after the intermediate crossover process is given below. First, the variable that controls the amount of mutation, scalem , is calculated.

C 4 : SS t ≥ SSlim ,

scalem = scalem − shrinkm ∗

SS t =

K 1 ∑

K

ri

(12)

i=1

∀t = t + 1 , . . . , t + p

(13)

3.2.6. Summary of rolling horizon model of DAS problem In a nutshell, the 3-Dimensional Dynamic Airspace Resectorization (3D-DARS) model in the rolling horizon framework can be summarized as follows. Input: Voronoi sites and altitude split (continuous-valued variables) Output: 3D Sector shapes Objective: 1. Minimize the standard deviation of monitoring workload between sectors across time (M) Constraints: C 1 − C 3: Implemented within and during each time period C 4: Implemented using the sector shapes in the previous time period (if available) 1. Voronoi sites must exist within the constrained airspace (C 1) 2. Distances between the crossing points and the sector boundaries > 10NM(C 2) 3. Maximum Sector Coordination Workloads are kept within user-defined limits (C 3) 4. Sector Similarity between different time periods are kept within user-defined limits (C 4) 3.3. Optimization solver The solver used for this implementation is an evolutionary algorithm known as the basic Genetic Algorithm (GA) (Holland, 1975). GA is a heuristic method that is developed from the idea of biological evolution. Solutions are modelled in chromosomes and they undergo changes through operations such as crossover, recombination, mutation, and selection. 3.3.1. Genetic algorithm In this experiment, the GA (Lin, 2011) undergoes selection, crossover and mutation operations to form a new population. GA follows a binary tournament selection process, where two random chromosomes are picked and the one with better fitness value will be kept. This selection process is performed as many times as the size of the population (Np ). Then, the population undergoes an intermediate crossover process described below. s1 = q1 + kc ∗ βc ∗ xc ∗ (q2 − q1 )

(14)

Gi Gmax

(16)

where scalem and shrinkm are user-defined parameters. Gi and Gmax represents the ith iteration and the maximum number of iterations respectively. scalem generally controls the amount of mutation and shrinkm allows one to decrease the amount of mutation towards the end of optimization. For the purpose of a local search, the range of scalem and shrinkm typically varies between [0.1, 0.5] and [0.5, 1.0] respectively. The following shows how each variable in the chromosome is mutated, given that a randomly generated value exceeds the mutation fraction of 2/nVars. s(i) = q(i) + scalem ∗ x(N(0, 1))

(17)

where s(i) and q(i) are the i variable of the child and parent chromosome respectively and x(N(0, 1)) is a normally distributed random number. After going through the binary tournament selection, intermediate crossover and Gaussian mutation operations, a new population is formed. This new population is then merged with the initial population, where the best Np solutions are kept, forming the final population for the current generation. This cycle then repeats for a user-defined number of generations. 3.3.2. Chromosome structure With the use of Voronoi diagrams for sector generation, the input variables for the GA algorithm are the site locations and the vertical splits. Fig. 2 provides a typical structure of a chromosome for the airspace sectorization problem in the GA solver for a single time period. The total number of 2D sectors is represented by K and the total number of vertical splits is indicated by N − K , giving a total of N sectors. For a total of N = 7 sectors with K = 4 areas (2D sectors), a total of 2 ∗ N + (N − K ) = 17 input φ variables are required. viλ and vi for i = 1, . . . , K represents the φ Voronoi sites used for generation of 2D sector shapes. viλ and vi for i = K + 1, . . . , N are points used to decide which sector should the vertical split occur in. In other words, the sector with Voronoi φ site that is closest to (vKλ +1 , vK +1 ) will experience a vertical split φ φ at alt1 . In the case where vKλ +1 , vK +1 and vKλ +2 , vK +2 points to the same sector, there would be 2 vertical splits in the selected sector at alt1 and alt2 . Given that the total number of time periods considered is 3, the total input variables for N = 7 sectors with K = 4 areas (2D sectors) is 17 ∗ 3 = 51.

6

C.S.Y. Wong, S. Suresh and N. Sundararajan / IFAC Journal of Systems and Control 11 (2020) 100076

Table 1 Parameter settings. Optimization settings

Value/Methodology

NGPM settings Population size (Np ) Maximum generations Selection Crossover Mutation

50/150 (Rolling horizon) 2000 Binary tournament Intermediate crossover (βc = 1.2) Gaussian (scalem = 0.1, shrinkm = 0.5)

Model settings Number of sectors (N) Number of areas (K ) Number of periods (p) (Rolling horizon) Max. coordination workload (clim ) Min. sector similarity (SSlim )

7 4 3 15 80%

4. Performance evaluation of the rolling horizon optimization approach In this section, an experimental study on the rolling horizon approach in comparison with the single interval optimization approach, i.e. not considering any future time periods during the optimization, is presented. Both approaches follow the same formulation of the problem reflected in Section 3.2.6. The rolling horizon approach considers M t +1 , M t +2 and M t +3 at period t during the optimization while the single interval approach only considers M t +1 . 4.1. Experimental details In this section, the details of the experimental study conducted using the Singapore FIR and the historical flight data are presented. The study has been conducted using a GA with a modified version of the software implementation, NGPM v1.4 (NSGA-II Program for MATLAB) (Lin, 2011). Table 1 shows the parameter settings used in this study. The population size is 50 and 150 for the single interval and rolling horizon approach respectively, varied according to the number of input variables. The maximum number of generations is set to ensure that the solutions found by the genetic algorithm converge. The ratio for immediate crossover is set to 1.2 (above one) for the diversification of solutions. The scale and shrink variables for mutation are set to act as a local search to improve the potentially sub-optimal solutions and allows for convergence towards the end of a run. The number of areas (K ) and sectors (N) were previously found to be favourable as 7 and 4 in Wong et al. (2019). A time interval of 1 h is analysed for this study. Typically, the traffic forecast up to 3 h is available. Therefore, p is chosen to be 3 in our experiment. In the initialization, the input variables are set to satisfy the spatial constraint (C 1). In the first time period of the rolling horizon approach, the input variables for each time period (1, 2, 3) are initialized with the same values, i.e. the same sector shapes are given for each time period at the start of the optimization. Sector similarity is applied in all circumstances except in the first time period. In the rolling horizon approach, the chromosome values of the previous time period sectorization are used for the measurement of sector similarity in the future time intervals. Further, the selected solution of the future time intervals (t + 2, t + 3) will be carried over to the optimization of the next time interval as an individual of the population.

Fig. 3. Air traffic flows used for experimental study.

4.2. Flight data description In this study, a set of historical flight data on 1st June 2015 obtained from FlightAware (2015) has been used to validate the rolling horizon approach. The flight data is made up of ADS-B data with 2–4 positions per minute depending on the coverage of the area and ANSP radar data (from governments) that are about 1 position per minute. However, the tracks provided only cover the location up to 108 deg in longitude. Hence, the flight data from the longitude of 108 to 116.3 degs was extrapolated from the given flight plans to reflect a better representative picture of the Singapore FIR. Furthermore, the missing information on

C.S.Y. Wong, S. Suresh and N. Sundararajan / IFAC Journal of Systems and Control 11 (2020) 100076 Table 2 Comparison of objective between the single interval and rolling horizon approach. Current period (t)

Single interval M t +1

0 1 2 3 4

Rolling horizon (p = 3) SS t +1

0.037 – 0.013 81% 0.087 80% 0.006 81% Unable to meet SSlim

M t +1

SS t +1

0.031 0.017 0.001 0.005 0.014

– 80% 81% 80% 83%

the altitude of certain flights was corrected using the historical information and flight plans. Only flight data above the flight level FL180 is considered for this study. Even though the Singapore FIR is a non-convex airspace, the Voronoi shapes are implemented on a cuboid shape of the latitude of [−0.5, 10.3], longitude of [102.10, 116.30] and altitude of [18 000, 45 000] ft for the sake of simplicity in this experiment. In the future, we can explore the use of non-convex shapes, recently proposed in Gerdes et al. (2018). A total of 5 h intervals of air traffic is being illustrated in Fig. 3. The 3D trajectories and 2D projections of the trajectories (on latitudes and longitudes) are provided on the left and right respectively. It can be observed that the traffic density increased rather significantly from Hour 0 to Hour 1. From the 2D projections of trajectories, the spread of traffic flow from Hour 1 to 2 can be seen to increase, with traffic flows from Hour 2 covering over a wider area. From Hour 2 to Hour 3, the traffic flow again changes, a greater number of intersecting flights can be observed. The traffic density in Hour 3 and 4 is quite similar, but the traffic flows are quite different. The traffic complexity components of the dynamic density metric are calculated using Eq. (2). The changes in heading, speed, and altitude (θHC , θSC , θAC ) are stored as locations where these changes occurred. The other factors (θMD5 , θMD10 , θCP25 , θCP40 , θCP70 ) are pegged to the flight and points at which they occurred. 4.3. Study results The results from the study is presented in two parts. First, the rolling horizon approach is presented along with the single interval optimization approach in Section 4.3.1. Next, a detailed analysis of the rolling horizon approach is given in Section 4.3.2. Lastly, the results are discussed in Section 4.3.3. 4.3.1. Comparison of the rolling horizon and single interval approaches The implemented solutions extracted from the rolling horizon approach are compared with the results obtained from the single interval approach in Table 2. In period t = 0, 1, 3, the M t +1 values for the single state and rolling horizon approach are comparable to each other. However, in the period t = 2, the solution obtained for the single interval approach is worse off than the rolling horizon approach. Fig. 4 provides the proposed sector shapes for period t = 0, 1, 2, to allow us better understand this observation. The sector shapes found in Hour 0 for both approaches are quite different, even though both solutions are able to obtain a sufficiently low standard deviation of monitoring workload (as reflected in Table 2). As the traffic density increases from Hour 0 to Hour 1, both approaches seem to be able to handle the traffic with very different sector shapes. This is because the vertical splits in 3D sectors occur in different areas in both approaches. In the rolling horizon approach, the 3 splits occur at 22,260 ft in the blue,

7

Table 3 Overall results of rolling horizon approach. Current interval (t)

M

M t +1

M t +2

M t +3

SS t +1

0 1 2 3 4 5 6 7 8 9

0.016 0.008 0.007 0.013 0.014 0.012 0.017 0.041 0.086 0.055

0.031 0.017 0.001 0.005 0.014 0.022 0.004 0.007 0.118 0.073

0.017 0.001 0.005 0.014 0.022 0.006 0.007 0.039 0.089 0.051

0.001 0.005 0.014 0.022 0.006 0.007 0.039 0.078 0.051 0.042

– 81% 80% 81% 83% 80% 84% 85% 86% 80%

29,630 ft in the yellow and 28,990 ft in the green areas, while the 3 splits occur at 21,260 ft in the blue, 28,950 ft in the yellow and 29,290 ft in the red areas instead. However, when the spread of traffic flows changes from Hour 1 to Hour 2, the single optimization approach (M t +1 = 0.087) starts to struggle in finding a good solution while meeting the sector similarity constraints. Eventually at t = 4, the single interval approach is unable to find a solution that meets the sector similarity constraint. 4.3.2. Detailed analysis of the rolling horizon approach The rolling horizon approach was ran for an interval of 10 h. Table 3 presents an overview of the results from the rolling horizon approach. Overall, feasible solutions meeting the given constraints were obtained in every interval. The mean standard deviation of the monitoring workload implemented at t, M t +1 is consistently below 0.05 until t = 8. Further, it can be observed that the standard deviation of monitoring workload, M 0+2 = M 1+1 and M 0+3 = M 1+2 . This similar occurrence is also observed in the next time intervals with the exception of M 6+1 , M 8+1 , M 8+2 and M 9+1 . This means that the solutions found in the earlier time interval are generally being selected and implemented in the future time intervals. The sector similarity SS t +1 of the solution obtained for each interval meets the sector similarity constraint SSlim of 80%, varying from 80% to 86%. The sector similarity SS t +1 tends to hover around 80% because the model is built to satisfy the minimum requirement of 80% rather than maximizing the sector similarity. At t = 8, there was a significant spike in the standard deviation of the monitoring workload in the proposed solution (M 8+1=9 = 0.118). In order to gain further insights, the proposed sector shapes at t = 7 and t = 8 are illustrated in Fig. 5. The proposed solution for M 9 at t = 7 and t = 8 are different, particularly in the yellow and blue regions. The solution for M 11 at t = 8 reveals that the previously proposed sector shapes in for M 9 at t = 7 might not be able to generate a feasible solution for M 11 . 4.3.3. Discussion This section provides some insights on the results presented above. The first observation is on the occurrence of vertical splits. They tend to be in regions with crossing traffic or near the bottom left corner of the airspace. This can be related to the dynamic density metric where traffic complexities such as crowded areas, potential conflicts and altitude changes are penalized with higher workloads. Thus, these vertical splits for the balance of workloads would naturally occur near the Singapore Changi airport with many ascending/descending flights near the bottom left corner of airspace and over regions with crossing traffic. Further, one can verify this by referring to the traffic flows found in Fig. 3. It can be seen from Table 2 that both approaches are able to meet the sector similarity limit constraint SSlim of 80% until

8

C.S.Y. Wong, S. Suresh and N. Sundararajan / IFAC Journal of Systems and Control 11 (2020) 100076

Fig. 4. Comparison of proposed sectorization between rolling horizon and single interval approaches in Hour 0, 1, 2 [Flight tracks are indicated by blue lines]. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

t = 3. It should be noted that the sector similarity measure is the average of the sector similarity of the 2D sector shapes given in Eq. (12). Hence, for the solution in Hour 1 of the rolling horizon approach, the high similarity for the blue (100%) and green (96%) areas are able to overcome the low similarity of the red area (52%). This might in turn allow a greater change in sector shapes than expected. In order to overcome this, one might consider using a higher SSlim . For a better range of solutions for a decision maker to choose from, the use of SSt +1 as an objective to maximize instead of a constraint would be explored in the future. The inability of the single interval approach to find a solution that satisfies the sector similarity constraints indicates that there could be many solutions suitable for the given traffic at the current time interval t. However, not all solutions are suitable for the traffic in future time intervals. The consideration of future time intervals helps to narrow down the pool of suitable solutions. Hence, it is effective in terms of planning in the DAS problem. It

allows one to continuously obtain feasible solutions in the future time horizons. The study also shows that there could be trade-offs in earlier time intervals, to implement lesser optimal solutions (in this case for t = 1) so that there are feasible solutions for future time intervals (t = 4). The detailed analysis continues to exemplify the significance of the rolling horizon approach in the planning of DAS, enabling one to obtain feasible solutions in future time intervals with gradual changes in sector shapes. Furthermore, the red, blue and green area consist of an unchanged overlapping airspace despite changes in sector shapes, giving the solutions a certain structure. This helps in the operational feasibility of solutions as air traffic controllers are usually trained and are more familiar with the traffic flows in a particular area. Therefore, it can be said that the consideration of future time intervals in the rolling horizon method helps to obtain better and feasible solutions over

C.S.Y. Wong, S. Suresh and N. Sundararajan / IFAC Journal of Systems and Control 11 (2020) 100076

9

Fig. 5. Proposed sector shapes with traffic flows at t = 7, 8. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

a time horizon and that the use of rolling horizon approach is appropriate for the DAS problem. 5. Conclusions The paper presents a rolling horizon framework for the DAS problem. Given the changing states of the air traffic flows, the rolling horizon model finds a sequence of gradually changing sector shapes that takes into account future time interval variations during the optimization. It also satisfies the given spatial, safety and sector similarity constraints. The performance of the rolling horizon approach was evaluated in comparison to the single time interval optimization approach using the historical flight data in the Singapore FIR over a finite time horizon. The single interval optimization faces problems in finding feasible sector shapes in the later time intervals. However, the rolling horizon model is able to provide feasible sector shapes with the trade-off of implementing lesser optimal solutions in the earlier time interval. Moreover, the sector shapes found in the rolling horizon approach are more similar between time intervals as when compared to the single interval method. This might help in the operational feasibility of the proposed solutions because air traffic controllers are familiar with the traffic flows in the given sector shapes. Hence, one can conclude that the rolling horizon model is suitable for the DAS problem to find feasible (gradually changing) sector shapes over a finite time horizon. Acknowledgement The authors wish to extend their thanks to the ATMRI:2014R8, Singapore, for providing financial support to conduct this study. References Berg, M. d., Cheong, O., Kreveld, M. v., & Overmars, M. (2008). Computational geometry: Algorithms and applications (3rd ed.). Santa Clara, CA, USA: Springer-Verlag TELOS.

Bullo, F., Cortés, J., & Martínez, S. (2009). Distributed control of robotic networks. In Applied Mathematics Series, Princeton University Press, Electronically available at http://coordinationbook.info. Chen, Y., Bi, H., Zhang, D., & Song, Z. (2013). Dynamic airspace sectorization via improved genetic algorithm. Journal of Modern Transportation, 21(2), 117–124. http://dx.doi.org/10.1007/s40534-013-0010-2. Delahaye, D., & Puechmorel, S. (2008). 3D airspace design by evolutionary computation. In 2008 IEEE/AIAA 27th digital avionics systems conference (pp. 3.B.6–1–3.B.6–13). http://dx.doi.org/10.1109/DASC.2008.4702803. Federal Aviation Adminstration (F. A. A. ) (2017). FAA aerospace forecast fiscal years 2017-2037. URL https://www.faa.gov/data_research/aviation/ aerospace_forecasts/media/FY2017-37_FAA_Aerospace_Forecast.pdf. FlightAware (2015). ADS-B flight data in Singapore regional airspace in june 2015. Fortune, S. Voronoi diagrams and delaunay triangulations. In Computing in euclidean geometry (pp. 225–265) arXiv:https://www.worldscientific.com/doi/ pdf/10.1142/9789812831699_0007. Gerdes, I., Temme, A., & Schultz, M. (2018). Dynamic airspace sectorisation for flight-centric operations. Transportation Research Part C: Emerging Technologies, 95, 460–480. http://dx.doi.org/10.1016/j.trc.2018.07.032, URL http: //www.sciencedirect.com/science/article/pii/S0968090X18310520. Holland, J. H. (1975). Adaptation in natural and artificial systems. an introductory analysis with applications to biology, control and artificial intelligence. Kicinger, R., & Yousefi, A. (2009). Heuristic method for 3D airspace partitioning: genetic algorithm and agent-based approach. In Proceedings of 9th AIAA aviation technology, integration and operations (ATIO) conference, Hilton Head, South Carolina (pp. 1–15). Laudeman, I. V., & Center., A. R. (1998). Dynamic density [microform] : an air traffic management metric (p. 1 v.). National Aeronautics and Space Administration, Ames Research Center ; National Technical Information Service, distributor Moffett Field, Calif. : Springfield, Va. Lin, S. (2011). NGPM0 – a NSGA-II program in matlab v1.4. http: //www.mathworks.com/matlabcentral/fileexchange/31166-ngpm-a-nsgaii-program-in-matlab-v1-4. URL http://www.mathworks.com/matlabcentral/ fileexchange/31166-ngpm-a-nsga-ii-program-in-matlab-v1-4. Michalewicz, Z., Krawczyk, J. B., Kazemi, M., & Janikow, C. Z. (1990). Genetic algorithms and optimal control problems. In 29th IEEE conference on decision and control, Vol. 3 (pp. 1664–1666). http://dx.doi.org/10.1109/CDC.1990.203904. Schultz, M., Gerdes, I., Standfuß, T., & Temme, A. (2017). Future airspace design by dynamic sectorization. In EIWAC. URL http://elib.dlr.de/115051/. Sergeeva, M. (2017). Automated airspace sectorization by genetic algorithm. optimization and control [math.OC] (Ph.D. thesis), Université Paul Sabatier (Toulouse 3). Sergeeva, M., Delahaye, D., & Mancel, C. (2015). 3D airspace sector design by genetic algorithm. In 2015 international conference on models and technologies for intelligent transportation systems (MT-ITS) (pp. 499–506). http://dx.doi.org/ 10.1109/MTITS.2015.7223300.

10

C.S.Y. Wong, S. Suresh and N. Sundararajan / IFAC Journal of Systems and Control 11 (2020) 100076

Sethi, S., & Sorger, G. (1991). A theory of rolling horizon decision making. Annals of Operations Research, 29(1–4), 387–416. http://dx.doi.org/10.1007/ BF02283607. Sharma, R., & York, G. W. (2018). Near optimal finite-time terminal controllers for space trajectories via sdre-based approach using dynamic programming. Aerospace Science and Technology, 75, 128–138. http://dx.doi.org/10. 1016/j.ast.2017.12.022, URL http://www.sciencedirect.com/science/article/pii/ S1270963817312026. Tang, J., Alam, S., Lokan, C., & Abbass, H. A. (2012a). A multi-objective approach for dynamic airspace sectorization using agent based and geometric models. Transportation Research Part C: Emerging Technologies, 21(1), 89–121. http://dx.doi.org/10.1016/j.trc.2011.08.008, URL http://www. sciencedirect.com/science/article/pii/S0968090X11001185. Tang, J., Alam, S., Lokan, C., & Abbass, H. A. (2012b). A multi-objective evolutionary method for dynamic airspace re-sectorization using sectors clipping and similarities. In 2012 IEEE congress on evolutionary computation (pp. 1–8). http://dx.doi.org/10.1109/CEC.2012.6253008. Tang, X.-M., Chen, P., & Li, B. (2016). Optimal air route flight conflict resolution based on receding horizon control. Aerospace Science and Technology, 50, 77–87. http://dx.doi.org/10.1016/j.ast.2015.12.024, URL http://www. sciencedirect.com/science/article/pii/S1270963815300262. Venugopalan, T. K., Wong, C. S. Y., Suresh, S., & Sundararajan, N. (2018). Simultaneous optimization of airway and sector design for air traffic management. Journal of Air Transportation, 26(1), 8–22. http://dx.doi.org/10.2514/1.D0090. Wang, H., Song, Z., Wen, R., & Zhao, Y. (2016). Study on evolution characteristics of air traffic situation complexity based on complex network theory. Aerospace Science and Technology, 58, 518–528. http://dx.doi.org/10. 1016/j.ast.2016.09.016, URL http://www.sciencedirect.com/science/article/pii/ S1270963816306897.

Weiler, K., & Atherton, P. (1977). Hidden surface removal using polygon area sorting. SIGGRAPH Computer Graphics, 11(2), 214–222. http://dx.doi.org/10. 1145/965141.563896. Wong, C. S. Y., Suresh, S., & Sundararajan, N. (2019). CDAS: A cognitive decision-making architecture for dynamic airspace sectorization for efficient operations. IEEE Transactions on Intelligent Transportation Systems, [ISSN: 1558-0016] 20(5), 1659–1668. http://dx.doi.org/10.1109/TITS.2018.2833151. Wong, C. S. Y., Venugopalan, T. K., & Suresh, S. (2016). A multi-objective approach for 3D airspace sectorization: A study on Singapore regional airspace. In 2016 IEEE symposium series on computational intelligence (SSCI) (pp. 1–8). http://dx.doi.org/10.1109/SSCI.2016.7850098. Xu, Y., Zhang, H., Liao, Z., & Yang, L. (2016). A dynamic air traffic model for analyzing relationship patterns of traffic flow parameters in terminal airspace. Aerospace Science and Technology, 55, 10–23. http://dx.doi.org/10. 1016/j.ast.2016.05.010, URL http://www.sciencedirect.com/science/article/pii/ S1270963816301808. Xue, M. (2010). Three dimensional sector design with optimal number of sectors. In Proceedings of the AIAA conference on guidance, navigation and control (GNC) and modeling and simulation technologies (MST). Zhang, Z., Li, J., & Wang, J. (2018). Sequential convex programming for nonlinear optimal control problems in uav path planning. Aerospace Science and Technology, 76, 280–290. http://dx.doi.org/10.1016/j.ast.2018.01.040, URL http://www.sciencedirect.com/science/article/pii/S1270963817304352.