# Sectorization and Configuration Transition in Airspace Design

Authors: Xiang Zou; Peng Cheng; Bang An; Jingyan Song
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Current airspace is sectorized according to some predefined rules that are not flexible. To facilitate utilizing the airspace more efficiently, methods to design sectors need to be promoted. In this paper, we propose an undirected graph cut-based approach that employs a memetic local search-embedded constrained evolution algorithm, NSGA-II, to generate nondominated airspace configurations. We also propose a new concave hull-based method to automatically depict sector boundaries. In addition, we also study the configuration transition problem. We define the similarity of the two different configurations and calculate their similarity with a bisection diagram and a minimum cost flow algorithm. We build a forward network to represent configuration transitions across several consecutive time periods and use multiobjective dynamic programming to determine a series of nondominated configuration links from the first period to the end. We test our approaches by simulation in high-altitude airspace controlled by Beijing Area Control Center. The results show that our sectorization method outperforms the current configuration in practice, providing a lower sector number, lower intersector flow, more balanced workload distribution among the different sectors, and no constraint violations, so that the proposed approach shows its significant potential as practical applications for dynamic airspace configuration.


---

## Body

## 1. Introduction and Literature Review

Airspace sectors are basic controlling units in Air Transportation Systems (Figure1). They were originally designed according to some predefined rules such as historical or geographic considerations or just according to experience. Sectors have essentially remained unchanged in terms of geometric shape and the total number of sectors inside a specific airspace. However, along with rapidly increasing air transportation, fixed sectors cannot accommodate varying traffic flows anymore; several problems have arisen, such as unbalanced workload distribution across different sectors, with overload in some sectors and very sparse flow density in others, and improper sector numbers, which means too many open sectors in off-peak time periods and too few sectors during busy times or too little flight time in a single sector for some flights.Figure 1
Configuration of Beijing Area Control.Original ideas to deal with the problem of fixed airspace structure is the “Merge and Divide” operation, meaning combining two or more adjacent sectors together when the traffic flow is low and splitting one sector into several during peak hours or choosing one airspace structure from a predefined experienced structure set [1–4]. However, this approach is not flexible enough because the boundaries of these sectors remain unchanged across different time periods. A more advanced concept, called Dynamic Airspace Configuration (DAC), was therefore proposed [5]. In DAC, both the boundaries of the sectors and the number of sectors are allowed to change according to varying traffic situations.One key issue in DAC is the sectorization problem, that is, how to divide an airspace into several sectors. The solution of a sectorization problem is always called the airspace configuration.To the best of our knowledge, the work by Delahaye et al. [6] may be one of the earliest studies to systematically research the sectorization problem, in which the author utilized a genetic algorithm to generate an optimal airspace configuration. Since then, many approaches have been developed. To summarize, relevant methods can be sorted into three categories [7]:(i)
Methods based on geometric computation.(ii)
Methods by cells (grids) growth (gathering) or by directly clustering trajectory points.(iii)
Methods based on undirected graph cuts.In the geometric computation category [8–13], approaches combining Voronoi diagrams with genetic algorithms were proposed in [8–10]. Tang et al. [11] used several kinds of geometric cuts such as bisection cuts and kd trees to split the airspace and compare different cutting methods.In the cell-growth category [11, 14–20], Brinton directly clustered trajectory points to form sectors [14]. Yousefi and Donohue divided an airspace into three layers with different altitude ranges [16, 17]. Each layer was discretized into hexagonal cells, with information about the controller workload. The hexagonal cells were then gathered into sectors. Based on [16], Drew utilized a boundary-smoothing method to eliminating jagged boundary segments [18]. Klein also divided an airspace into hexagonal cells [19], but, in his approach, sectors grew up from a set of seeding cells.The third category [20–24] is in fact another kind of clustering approach, but it is based on a weighted undirected graph and uses one subgraph to represent a sector. Li et al. constructed a weighted graph model that accurately represents the air-route network [21]. The sectorization problem was then formulated as a graph cut problem and solved by iterative spectral bisection. Martinez et al. proposed a method based on a weighted graph and a grid [22] and also utilized spectral bisection to cut the graph. Chen and Zhang proposed a spectral clustering-based approach to clustering vertices [23]. The spectral clustering solution was further refined by the ODLB algorithm and another heuristic algorithm to get better performance in terms of workload balancing. Trandac et al. proposed a method based on Constraint Programming [24].In [25], Zelinski gave a comprehensive comparison of different approaches. The results showed completely different sector shapes according to the different approaches. The performance of these approaches was evaluated, which revealed their strengths and weaknesses. To summarize, methods using geometric computation are simple and straightforward, but they are optimally inferior because the methods used to cut the plane or space are limited. Although the second category may be the best in terms of workload balancing, it can hardly handle other objectives or constraints in sectorization problems. Approaches based on undirected graph cuts show great potential in managing multiple objectives and constraints. However, these approaches have to face two critical issues. The first one is the validity of the graph cut method, which must ensure conformance to all constraints of the problem. The second one, which is much more challenging, is how to depict exact sector boundaries based on the generated subgraphs. In this paper, we focus on high-altitude airspace in en route areas. The top and bottom covers of this kind of airspace match completely, so that horizontal and vertical partitioning operations can be executed separately [7]. Because the operations for vertical partitioning are always much simpler than the horizontal ones [7, 10], in our work, we only carry out sectorization in the horizontal plane, and we call the sectorization problem solved here the 2D en route sectorization problem (2D ERSP).The 2D ERSP is a typical multiobjective optimization problem (MOOP) with several objectives, such as balancing the workload across different sectors and minimizing intersector flow, and several constraints, such as no overload in any sector and no reentry of any single sector for any flight.The first main contribution of this paper is that we propose an approach that comprehensively handles multiple objectives and ensures no violations against any constraint. This approach belongs to the undirected graph cut category, and it is based on the reasonable airspace model and the Constrained NSGA-II [26] accompanied by memetic local search [27]. In addition, thanks to the evolution procedure, we do not need to predefine desired sector numbers as many previous studies have had to do.Different from geometric computation-based and cell-growth-based approaches, undirected graph cut-based methods cannot naturally generate sector boundaries after vertex clustering. However, many studies of this type directly neglect the boundary depiction problem, only giving vertex clustering results. Among the very limited results in undirected graph cut methods that seriously consider the boundary depiction problem, Chen and Zhang used Voronoi polygons of vertices to form sector boundaries [23]. Li et al. proposed a method based on a shortest-path algorithm [21], and Trandac et al. usedconstrained triangulation [24] to form sector boundaries. These two studies both used manual selection. In fact, to the best of our knowledge, except for the methods employing Voronoi polygons, other boundary depiction methods serving undirected graph cut-based approaches generally require manual selection or manual modification.Hence, the second main contribution of this paper is its proposal of a new concave hull-based approach to automatically depict exact sector boundaries.In the DAC, although airspace configuration should remain fixed for a particular time period for practical reasons, such as the demand of air traffic controllers [2–4], it can and should be changed to match varying traffic situations. Generally, the change of airspace configuration should be carried out in conjunction with the change in air traffic controllers [2]. Hence, another important issue in DAC is how to avoid sharp changes in airspace configuration in order to maintain steady air traffic. This problem is called the configuration transition problem (CTP). Unfortunately, like the boundary depiction problem, the CTP is also seldom reported in the existing literatures. To our best knowledge, only the results in [4, 28] presented several rules to minimize the workload caused by the configuration transition.The third main contribution of this paper is that it represents a serious study of the CTP. We first define the similarity of the two configurations and calculate their similarity with a bisection diagram and a minimum cost flow algorithm [29]. Next, we build a forward network to represent configuration transitions across several consecutive time periods. Then, we take advantage of multiobjective dynamic programming to find a series of nondominated configuration links from the first period to the end.This paper is organized as follows. In Section2, we give a comprehensive description of the 2D ERSP, including its objectives and constraints. In Section 3, we discuss how to model the airspace network into an undirected graph and how to integrate traffic information into the model. In Section 4, approaches used to solve the 2D ERSP and to depict sector boundaries are presented. In Section 5, procedures to solve the CTP are stated in detail. The performance of the proposed approaches is tested and analyzed in Section 6. Conclusions and future works are discussed in the last section.


## 2. Description of the 2D ERSP

We describe the 2D ERSP by three objectives and five constraints that are commonly mentioned in the majority of the literatures working on the sectorization problem.The three objectives are as follows.Obj1. Balancing intrasector workload across different sectors.Obj2. Minimizing aggregate intersector workload (coordination workload).Obj3. Minimizing aggregate occurrences of the violations of the minimum staying time constraint.Obj1 andObj2 are general objectives that are nearly mentioned in all relevant literatures. Nevertheless,Obj1 always acts as the core objective, whileObj2 is sometimes not really considered when solving the sectorization problem in previous studies. Besides,Obj3 is sometimes considered as a constraint in literatures like [20]. However, when it acts as a constraint, it is always treated as “soft constraint”; that is, it may not be strictly satisfied. Therefore, we model it as an objective aiming at reducing its violations.Along with these three objectives, the 2D ERSR has five constraints as follows.Con1. Con1 is the connectivity constraint; that is, any sector must be connected, and no sector can consist of two or more separate parts. This is the basic constraint which can never be violated.Con2. Con2 is the convexity constraint; that is, no flight can enter one single sector more than once. Practically, this is also a strict constraint. However, it cannot be ensured in some cases.Con3. Con3 is the instantaneous flow constraints; that is, no sector can be overload, and the number of flights in one single sector at any moment cannot exceed a threshold. This constraint is often neglected in some literatures.Con4. Con4 is the minimum distance between fixes and boundaries (MDFB) constraint; that is, the distance between any fix and any segment of the sector boundary should not be below a minimum threshold. Few literatures have had proposed powerful methods forCon4.Con5. The sector boundary should be as compact as possible. Jagged boundary is unwanted.


## 3. Airspace Model

### 3.1. Initial Undirected Graph Construction

Basic elements in airspace are waypoints, navaids, ordinary intersections, and air routes connecting them. This natural network structure makes it very convenient to represent airspace by a graph.We model airspace into an undirected graphG(V,E), whereV={v1,v2,…,vn} is the set of vertices that a vertex stands for a waypoint, or a navaid, or an ordinary intersection, and E={eij∣eij=1,vi,vj∈V} is the set of edges. We set eij=1, if vi is connected with vj by a segment of air route or if the Voronoi polygons of vi and vj have a common side.It is an expansion to existing literatures considering the situation of common side of Voronoi polygons when constructing the graph, in which only air routes are treated as edges of the undirected graph. One of the major reasons of introducing these extra edges is to cover the shortcomings in the connectivity check in the existing graph cut-based methods. In Figure2(a), we assume a regional airspace consisting of three vertices and two route segments. If we want to cluster va and vb into one sector and split vc alone into another sector, the connectivity constraint in existing methods would be violated while it is satisfied in the expanded graph (see Figure 2(b)) proposed here, because va and vb are treated as connected to the common side of their Voronoi polygons. This situation is possible in practice, such as a regional hub-spoke structure in which all surrounding fixes are connected to the hub fix with air routes while no routes exist between these surrounding fixes. To relax the burden of the controller in charge of the airspace near the hub fix, one general choice is to divide the airspace into more parts so that some route-unconnected fixes may be included in one sector. This dividing action could lead connectivity violation based on existing type of graph, while it still functions well if our expanded graph is employed.Figure 2
Edge between Voronoi neighborhood.
(a) (b)


### 3.2. Graph Simplification and Modification

We set four rules to simplify and modify the initial undirected graph.Rule 1.
If any fix is only crossed by one air route and there is no difference between the inbound course and the outbound course, directly delete this fix. In this situation, the flow density on both sides of the fix is nearly the same so that it is not necessary to consider these two route segments separately.Rule 2.
If several fixes are less than twice of the MDFB to each other, merge them as a fix union. This rule is to avoid any generated sector boundaries crossing the gap between two fixes that are too near to each other. In other words, it is to satisfy the MDFB constraint. For any fix union, we represent it with a “union-point” at the geometric center of the union and define the radius of the “union-point” asr=maxdis+MDFB, where maxdis is the largest distance from the “union-point” to any member fixes in the fix union. Edges completely inside one fix union are deleted; edges connecting fixes belonging to different unions are represented by edges to connect the corresponding “union-points.” After Rule 2, we set each “union-point” as a vertex of the modified undirected graph.Rule 3.
If any edgeeij is less than the distance of MDFB from vertex vk that is not one of its endpoints, we split eij into two segments, eik and ekj.Rule 4.
If any edgeeij representing Voronoi neighborhood intersects with any other edge representing air route, delete eij. This action is to avoid generating unnecessary vertices.


### 3.3. Traffic Information Representation

Many previous studies only employ flow density information, such as the number of flights or the number of trajectory points, to describe the traffic situation. Although the flow density is indeed the most important factor when evaluating the workload of air traffic controllers, it is oversimplified in some extent. In fact, several comprehensive metrics [30–32] evaluating airspace complexity or air traffic controllers’ workload have been proposed. These metrics can be identified into four types: flow density, altitude change, speed change, and course change [2, 31, 32]. Accurately calculating each metric is very time-consuming, especially when solving the sectorization problem that needs iterative optimization. Hence, simplification or approximation must be used if we want to consider all these metrics simultaneously. There are examples such as representing several metrics with one integrated metric [2, 8], ignoring several metrics in the optimizing objectives [14], limiting the size of the target airspace and using much simpler airspace model [8], and precalculating the metrics in some basic airspace cells and directly using these values in the optimization process [15]. In the proposed approach, we introduce ten matrices to store these four types of information.(1) Flow Density Matrix F D N × N . This matrix stores the number of flights flying through each edge of G, where N is the number of vertices in G. If k flights fly through eij, FDij=k. Besides, if eij=1 and k=0, we still set FDij=ε, where ε is a small positive value, to remain the connectivity of vi and vj.(2) Symmetric Altitude Change Matrix S A C N × N . S A C stores the information about altitude changes. Suppose a set of trajectory points of flight f that are projected on eij, denoted as tpijf. We calculate the standard variance of the altitudes of tpijf, denoted as acijf. Define the set of flights flying through eij as Fij, and we have(1)SACij=SACji=∑f∈Fijtpijfacijf,where tpijf is the count of trajectory points in tpijf. Other elements in SAC with eij=0or eij=1, and FDij=εare all set to 0. As G is an undirected graph, SAC is symmetric.(3) Asymmetric Altitude Change Matrix A A C N × N . A A C also stores the information about altitude changes. We divide trajectory points that are projected on edge eij into two parts, tpij,if and tpij,jf. The trajectory points nearer to vi than vj are put into tpij,if, while the others are put into tpij,jf. Then, we calculate the standard variances of these two sets, denoted as acij,if and acij,jf. Hence, (2)AACij=∑f∈Fijtpij,ifacij,if,AACji=∑f∈Fijtpij,jfacij,jf.(4) Symmetric Speed Change Matrix S S C N × N. Similar to those notations in SAC, we define scijf and thus(3)SSCij=SSCji=∑f∈Fijtpijfscijf.Other elements inSSC with eij=0or eij=1 and FDij=εare all set to 0.(5) Asymmetric Speed Change Matrix A S C N × N. Similar to those notations in AAC, (4)ASCij=∑f∈Fijtpij,ifASCij,if,ASCji=∑f∈Fijtpij,jfASCij,jf.(6) Symmetric Trajectory Point Count Matrix S T P N × N. STP stores the count of trajectory points that are projected on each edge: (5)STPij=STPji=∑f∈Fijtpijf.(7) Asymmetric Trajectory Point Count Matrix A T P N × N. Consider (6)ATPij=∑f∈Fijtpij,if,ATPji=∑f∈Fijtpij,jf.(8) Flight Course Change Matrix F C C F × N. F is the number of flights considered in the sectorization problem and N is the number of vertices in G. Recall the rules in Section 3.2. In fact, any vertex in G represents a fix union. Hence, for any flight-vertex pair, such as f and vi, we set FCCfi as the aggregate course change times when flight f flies through vi.(9) Vertex Arrival Time Matrix V A T F × N. The element VATfi is defined as the earliest arriving time of flight f at any fix that belongs to vi.(10) Vertex Leaving Time Matrix V L T F × N. The element VLTfi is defined as the latest leaving time of flight f from any fix that belongs to vi.


## 3.1. Initial Undirected Graph Construction

Basic elements in airspace are waypoints, navaids, ordinary intersections, and air routes connecting them. This natural network structure makes it very convenient to represent airspace by a graph.We model airspace into an undirected graphG(V,E), whereV={v1,v2,…,vn} is the set of vertices that a vertex stands for a waypoint, or a navaid, or an ordinary intersection, and E={eij∣eij=1,vi,vj∈V} is the set of edges. We set eij=1, if vi is connected with vj by a segment of air route or if the Voronoi polygons of vi and vj have a common side.It is an expansion to existing literatures considering the situation of common side of Voronoi polygons when constructing the graph, in which only air routes are treated as edges of the undirected graph. One of the major reasons of introducing these extra edges is to cover the shortcomings in the connectivity check in the existing graph cut-based methods. In Figure2(a), we assume a regional airspace consisting of three vertices and two route segments. If we want to cluster va and vb into one sector and split vc alone into another sector, the connectivity constraint in existing methods would be violated while it is satisfied in the expanded graph (see Figure 2(b)) proposed here, because va and vb are treated as connected to the common side of their Voronoi polygons. This situation is possible in practice, such as a regional hub-spoke structure in which all surrounding fixes are connected to the hub fix with air routes while no routes exist between these surrounding fixes. To relax the burden of the controller in charge of the airspace near the hub fix, one general choice is to divide the airspace into more parts so that some route-unconnected fixes may be included in one sector. This dividing action could lead connectivity violation based on existing type of graph, while it still functions well if our expanded graph is employed.Figure 2
Edge between Voronoi neighborhood.
(a) (b)


## 3.2. Graph Simplification and Modification

We set four rules to simplify and modify the initial undirected graph.Rule 1.
If any fix is only crossed by one air route and there is no difference between the inbound course and the outbound course, directly delete this fix. In this situation, the flow density on both sides of the fix is nearly the same so that it is not necessary to consider these two route segments separately.Rule 2.
If several fixes are less than twice of the MDFB to each other, merge them as a fix union. This rule is to avoid any generated sector boundaries crossing the gap between two fixes that are too near to each other. In other words, it is to satisfy the MDFB constraint. For any fix union, we represent it with a “union-point” at the geometric center of the union and define the radius of the “union-point” asr=maxdis+MDFB, where maxdis is the largest distance from the “union-point” to any member fixes in the fix union. Edges completely inside one fix union are deleted; edges connecting fixes belonging to different unions are represented by edges to connect the corresponding “union-points.” After Rule 2, we set each “union-point” as a vertex of the modified undirected graph.Rule 3.
If any edgeeij is less than the distance of MDFB from vertex vk that is not one of its endpoints, we split eij into two segments, eik and ekj.Rule 4.
If any edgeeij representing Voronoi neighborhood intersects with any other edge representing air route, delete eij. This action is to avoid generating unnecessary vertices.


## 3.3. Traffic Information Representation

Many previous studies only employ flow density information, such as the number of flights or the number of trajectory points, to describe the traffic situation. Although the flow density is indeed the most important factor when evaluating the workload of air traffic controllers, it is oversimplified in some extent. In fact, several comprehensive metrics [30–32] evaluating airspace complexity or air traffic controllers’ workload have been proposed. These metrics can be identified into four types: flow density, altitude change, speed change, and course change [2, 31, 32]. Accurately calculating each metric is very time-consuming, especially when solving the sectorization problem that needs iterative optimization. Hence, simplification or approximation must be used if we want to consider all these metrics simultaneously. There are examples such as representing several metrics with one integrated metric [2, 8], ignoring several metrics in the optimizing objectives [14], limiting the size of the target airspace and using much simpler airspace model [8], and precalculating the metrics in some basic airspace cells and directly using these values in the optimization process [15]. In the proposed approach, we introduce ten matrices to store these four types of information.(1) Flow Density Matrix F D N × N . This matrix stores the number of flights flying through each edge of G, where N is the number of vertices in G. If k flights fly through eij, FDij=k. Besides, if eij=1 and k=0, we still set FDij=ε, where ε is a small positive value, to remain the connectivity of vi and vj.(2) Symmetric Altitude Change Matrix S A C N × N . S A C stores the information about altitude changes. Suppose a set of trajectory points of flight f that are projected on eij, denoted as tpijf. We calculate the standard variance of the altitudes of tpijf, denoted as acijf. Define the set of flights flying through eij as Fij, and we have(1)SACij=SACji=∑f∈Fijtpijfacijf,where tpijf is the count of trajectory points in tpijf. Other elements in SAC with eij=0or eij=1, and FDij=εare all set to 0. As G is an undirected graph, SAC is symmetric.(3) Asymmetric Altitude Change Matrix A A C N × N . A A C also stores the information about altitude changes. We divide trajectory points that are projected on edge eij into two parts, tpij,if and tpij,jf. The trajectory points nearer to vi than vj are put into tpij,if, while the others are put into tpij,jf. Then, we calculate the standard variances of these two sets, denoted as acij,if and acij,jf. Hence, (2)AACij=∑f∈Fijtpij,ifacij,if,AACji=∑f∈Fijtpij,jfacij,jf.(4) Symmetric Speed Change Matrix S S C N × N. Similar to those notations in SAC, we define scijf and thus(3)SSCij=SSCji=∑f∈Fijtpijfscijf.Other elements inSSC with eij=0or eij=1 and FDij=εare all set to 0.(5) Asymmetric Speed Change Matrix A S C N × N. Similar to those notations in AAC, (4)ASCij=∑f∈Fijtpij,ifASCij,if,ASCji=∑f∈Fijtpij,jfASCij,jf.(6) Symmetric Trajectory Point Count Matrix S T P N × N. STP stores the count of trajectory points that are projected on each edge: (5)STPij=STPji=∑f∈Fijtpijf.(7) Asymmetric Trajectory Point Count Matrix A T P N × N. Consider (6)ATPij=∑f∈Fijtpij,if,ATPji=∑f∈Fijtpij,jf.(8) Flight Course Change Matrix F C C F × N. F is the number of flights considered in the sectorization problem and N is the number of vertices in G. Recall the rules in Section 3.2. In fact, any vertex in G represents a fix union. Hence, for any flight-vertex pair, such as f and vi, we set FCCfi as the aggregate course change times when flight f flies through vi.(9) Vertex Arrival Time Matrix V A T F × N. The element VATfi is defined as the earliest arriving time of flight f at any fix that belongs to vi.(10) Vertex Leaving Time Matrix V L T F × N. The element VLTfi is defined as the latest leaving time of flight f from any fix that belongs to vi.


## 4. Solve the 2D ERSP

### 4.1. Mathematical Expression

As stated before, the airspace is modeled as an undirected graph,G. Solving the 2D ERSP is to cut G into k subgraphs, GsVs,Es, where s∈S=1,2,…,k, (7)G=⋃s=1kGs,V=⋃s=1kVs,Vi∩Vj=∅,i,j∈S,i≠j.Each subgraph stands for one sector.


#### 4.1.1. Objectives

(1) Balancing the Intrasector Workload. Consider (8) min ⁡ 1 k ∑ s = 1 k W L G s - W L - 2 . W L G s is the intraworkload of sector s and WL- is the mean value of WLGs, s∈S. WLGs is a weighted and scaled summation of the values of four metrics of sector s, including information about the flow density, the altitude change, the speed change, and the course change: (9)WLGs=wFDFD~s+wacAC~s+wscSC~s+wccCC~s,where wFD, wac, wsc, and wcc are weights of these four parts.We defineEs′ as the set of intrasector edges of sector s, where Es′=eij∣eij=1,vi,vj∈Vs and Es′′ are the set of intersector edges between s and other sectors, where Es′′=eij∣eij=1,vi∈Vs,vj∉Vs; then(10)FDs=∑eij∈Es′STPij+∑eij∈Es′′ATPij,ACs=∑eij∈Es′SACij+∑eij∈Es′′AACijFDs,SCs=∑eij∈Es′′SSCij+∑eij∈Es′′ASCijFDs,CCs=NTsNFs,where NTs is the total course change times of flights when flying inside sector s and NFs is the total number of flights flying through s. Hence, CCs stands for average course change times of single flight in sector s. Suppose the flight set considered in the sectorization problem as F, and we have(11)NTs=∑f∈Fvi∈VsFCCfi.Because of different ranges of the four parts in (9), and a more meaningful understanding of the summation, we normalize them to the range of [0,1]. We set four standard values, FDsd, ACsd, SCsd, and CCsd, that are large enough and do not change during the optimization. Then, we define(12)FD~s=FDsFDsd,AC~s=ACsACsd,SC~s=SCsSCsd,CC~s=CCsCCsd.(2) Minimizing the Aggregate Coordination Workload. We use the density of intersector flow to stand for the coordination workload. This objective can be expressed by(13)min⁡∑s≠rs,r∈SWLcoGs,Gr,where WLco(Gs,Gr) is the intersector flow density between Gs and Gr. It is the number of flights flying across Gs and Gr:(14)WLcoGs,Gr=∑vp∈Vsvq∈VrFDpq.(3) Minimizing the Occurrences of the Violations of the Minimum Staying Time Constraint. If any flight stays in the sector less than a minimum time threshold after it enters a sector, the minimum staying time constraint is violated.Obj3 is to minimize the occurrences of these situations.As we do not depict exact sector boundaries during the optimization process, it is impossible to accurately calculate the staying time of any flight. We approximate the accurate staying time as follows.Still suppose the flight path of flightf as fpf=vf,1,vf,2,…,vf,k, vf,1,vf,2,…,vf,k∈V, where vf,i is the ith vertex along the path. Then, we calculate the staying time of flight f in sector sas follows: (15)dtfs=12VATf,a-VLTf,a-1+12VATf,e+1-VLTf,e+VLTf,e-VATf,e+∑vf,m∈Vsvf,m+1∈VsVATf,m+1-VATf,m,where va and ve are the first and the last vertices inside sector s along fpf. The first term of (15) stands for the time spent to reach va after entering the sector. The second term stands for the time spent to leave the sector after leaving ve. The third term stands for the time spent in ve. The fourth term stands for other times spent inside sector s.We define the set of all considered flights asF and set mdtvfs=1, if any flight f∈F violates the minimum staying time constraint in sector s. Hence,Obj3 can be expressed as(16)min⁡∑f∈F∑s∈Smdtvfs.


#### 4.1.2. Constraints

(1) Connectivity Constraint. Each subgraph representing one sector must be self-connected.(2) Flight-Route-Based Convexity Constraint. For any flight f, fpf=vf,1,vf,2,…,vf,k, vf,1,vf,2,…,vf,k∈V. If vf,i,vf,j∈Vs, i≤j, for ∀vf,m, i≤m≤j, vf,m∈Vs.(3) Instantaneous Flow Constraint. Similar to the expressions above, we suppose va and ve as the first and last vertex inside sector s along fpf. Hence, the entry time of sector s can be denoted as mid(VLTa-1,VATa), which is the midtime of the time leaving va-1 and the time reaching va. Similarly, we denote the departure time of sector s as mid(VLTe,VATe+1). With the entry time and the departure time for every flight, we can calculate the number of flights in any sector at any moment. Suppose the maximum instantaneous flight number in sector s as Instflows and the instantaneous flow threshold as Threinsflow, andCon3 can be expressed as(17)Instflowf,s≤Threinsflow,∀s.Con4 andCon5 are not considered in the optimization phase. But they are ensured by the airspace modeling in Section 3 and the boundary depiction method in Section 4.3.


### 4.2. Constrained NSGA-II Based Evolution

Some objectives and constraints of the 2D ERSP cannot be written in explicit mathematical formulas [2, 24]. For example, consider the constraint that no flights can enter any single sector more than once. It is affected by the trajectory of flights and the boundary of a sector, which is too complex to express this constraint by an explicit inequality. Hence, it is impossible to construct the 2D ERSP into a typical form of mathematical programming below:(18)Min⁡fxAx≤bx∈X.As evolution algorithms can still work even if the constraints can only be expressed by text descriptions, they are good choices to solve the 2D ERSP. The Constrained Nondominated Sorting Genetic Algorithm II (CNSGA-II) [26] is a fast and elitist multiobjective genetic algorithm that is very appropriate for the 2D ERSP, considering its characteristics of multiobjectives and the strong capability in handling nonclassic constraints. We utilize the CNSGA-II accompanied with the memetic local search [27] to solve the 2D ERSP. The main flow diagram of the evolution algorithm is shown in Figure 3.Figure 3
Diagram of the evolution algorithm.


#### 4.2.1. Chromosome

Every chromosome hasN genes representing N vertices in graph G. The value of each gene is the index of the sector to which the related vertex belongs. For example, in Figure 4, the third vertex is assigned to sector 2 and the 50th vertex is assigned to sector 5.Figure 4
Chromosome.We set a common largest possible index of sectors for every individual; that is, the possible maximum number of sectors is fixed. On one hand, this number is large enough so that it has no impact on the results. On the other hand, because we only set the possible maximum number, indexes of sectors in different individuals are totally independent, which means that two sectors with the same index in different individuals may refer to different sectors.


#### 4.2.2. Reproduction

Tournament selection is used to select individuals in parent generation.


#### 4.2.3. Crossover

Two individuals can crossover only when they have large enough hamming distance. In single crossover operation, two randomly chosen genes are exchanged. To avoid generating individuals that violates connectivity constraint as much as possible, the crossover can be granted only when at least one of the vertices (genes) is moved to an adjacent sector or stays in original sector after the crossover.


#### 4.2.4. Mutation

In single mutation operation, one gene is randomly selected. The value of the gene is randomly changed and it is bounded by the value of the maximum sector number. However, similar to the crossover operation, only when the related vertex is changed to an adjacent sector, the mutation operation can be granted.


#### 4.2.5. Memetic Local Search

The memetic algorithm [27] is a local search method embedded in one global search iteration of evolution algorithms that can help improving the individuals in one generation. In this paper, two heuristic polices of memetic local search are applied.Policy A. For any sector, if FDs>FDup, where FDup is an upper bound, we split the sector into two or more parts. For the subgraph Gs corresponding to sector s, we build its similarity matrix WNs×NsGs, where Ns is the number of vertices that belong to Gs. For any element wijGs, 1≤i≤Ns, 1≤j≤Ns, suppose the indices of its two related vertices in graph G as p and q, and then wijGs=FDpq. We add up all elements of each row of WNs×NsGs to build a diagonal matrix DGs, where diiGs=∑j=1NswijGs. We define theLaplacian matrix of Gs, LGs=DGs-WGs. Besides, we introduce three lemmas in the Spectral Graph Theory [33].Lemma 1.
For a weighed undirected graph with nonnegative edge weights, itsLaplacian matrix has at least one eigenvalue with value of 0.Lemma 2.
The number of isolated subgraphs of a weighed undirected graph with nonnegative edge weights equals the number of 0 eigenvalues of itsLaplacian matrix.Lemma 3.
For a weighed undirected graph with nonnegative edge weights, denoted asG, suppose the number of 0 eigenvalues of itsLaplacian matrix as k. Then, the dimension of the eigenspace of 0 eigenvalues is k. This linear subspace is spanned by such a group of base vectors: 1s1,1s2,…,1sk, where 1si(1≤i≤k) is the indicator vector of ith isolated subgraph, denoted as Gi. If a vertex is in Gi, the corresponding element in 1si is 1; otherwise, it is 0.With these three lemmas, we calculate the eigenvalues ofLGs and the corresponding eigenvectors. If the multiplicity of 0 eigenvalues, k, is larger than 1, we just divide Gs into k parts according to the base vectors. If k=1, we get the eigenvector corresponding to the second smallest eigenvalue, calledFiedler Vector [34]. Then, the vertices of Gs are clustered into two groups. One is with the positive elements ofFiedler Vector and the other is with the negative ones. Vertices corresponding to the elements with value of 0 are randomly assigned to one of these two groups. Hence, when k=1, we split Gs into two parts.Policy B. Policy B is operated after the division by Policy A. For any individual, if FDs<FDlow for sector s, where FDlow is a lower bound, we treat the flow density in this sector as too low. Thus, we merge s into one of its adjacent sectors, assumed as s′. Values of genes (vertices) belonging to s are changed to the index of s′.To avoid loss of diversity and being trapped in regional optima, the memetic local search is operated with gradually changing intensity and frequency [35]. The frequency refers to the gap of generations between two callings of the local search, and the intensity refers to the proportion of individuals in one single generation that are applied with the local search. At the beginning of the evolution, the qualities of the individuals are still low, and, in general, we set low frequency and intensity. Along with the evolution, the searching space approaches the optimal solutions, and thus we need more accurate search. Therefore, the frequency and intensity of the local search both increase.


#### 4.2.6. Fitness Functions

Fitness functions correspond to the expressions in formulas (8), (13), and (16).


#### 4.2.7. Constraint Check

Connectivity. Similar to those stated in “Policy A” of the memetic local search, for each subgraph Gs (sectors), we calculate itsLaplacian matrix LGs and check the multiplicity of 0 eigenvalues. If the number is larger than 1, we record one more time of violation of the connectivity constraint for current individual.Flight-Route-Based-Convexity. Still suppose the flight path of flight f as vf,1,vf,2,…,vf,k. If vf,p∈Gs and vf,p+1∉Gs, we check if vf,q∈Gs, q>p+1. Once the condition vf,q∈Gs, q>p+1 is met, we record one more time of the violation of the convexity constraint for current individual.Instantaneous Flow Constraint. We check the validation of inequality (17). If inequality (17) is unsatisfied, we record one more time of the violation of the instantaneous flow constraint for current individual.


#### 4.2.8. Sort and Elitism

In the CNSGA-II, the elitism operation is executed for every generation. The dominating relationship between two individuals (solutions),i and j, is defined below:(1)
i≺j, if i is feasible (no constraint violation) while j is not.(2)
When solutionsi and j are both feasible, i≺j if the values of all fitness functions of i are not more than that of j and at least one of such values is smaller than that of j.(3)
When solutionsi and j are both infeasible, we compare their times of constraint violations in the following order of decreasing priority: connectivity > convexity > instantaneous flow. That is, i≺j, if the connectivity constraint violating times of solution i is less than that of solution j. Or, if the connectivity constraint violating times are equal, we continue to compare the violating times of the convexity constraint and so on.(4)
When solutioni and j are both infeasible and they are equivalent in terms of the violating times of all three constraints, i≺j if the values of all fitness functions of i are not more than that of j and at least one such value is smaller than that of j.


#### 4.2.9. Initial Generation

Most previous studies using evolutionary algorithms always concern about how to refine the “genetic operations” (such as the crossover or mutation rules) or try to propose heuristics to improve the evolution procedure. However, for those problems that the gross optima or a boundary of the gross optima cannot be preliminarily determined (such as the 2D ERSP), no matter how excellent a heuristic or an improved operation seems to be, we are unable to determine whether the gross optima can be found with the help of these heuristics. In these problems, we cannot guarantee reaching the gross optima from any type of initial generations. Hence, we have to realize that different initial generations have different impact on final solutions. However, discussions about the initial generation of evolution algorithms and its impact to the final solutions are absent in relevant literatures.In this paper, we propose four types of initial generations and compare them with the corresponding nondominated solution sets produced by the CNSGA-II.The first type is the most general choice that just randomly produces the initial generation. In following text, we call this caseRandomly Initialization, or RI.The second type takes advantage of the current configuration in practice. We represent the current configuration with a chromosome and randomly mutate half of the genes several times to produce the initial generation. In following text, we call this caseCurrent Initialization, or CI.Thirdly, we employ a method based on the normalized spectral clustering algorithm (NSCA) [36]. The NSCA is a data clustering algorithm that shows its talent in the graph N-cut problem if each data point is represented by one vertex in a weighted graph and the weights of edges stand for similarities between two data points [37]. We choose the NSCA because of two objectives of the 2D ERSP, balancing intrasector workload and minimizing gross intersector workload, which are very similar to the goals of the undirected graph N-cut problem. We first change the values of elements in matrix FD with value of ε to 0 and call the modified density matrix FD′. Then, FD′ is defined as the similarity matrix of G; that is, WG=FD′ and we construct a diagonal matrix DG, where diiG=∑j=1NwijG. Following steps of the method based on the NSCA are listed in Procedure 1. Note that diiG may equal to 0 if no aircraft flies through any routes connecting to vertex vi. This may lead singularity of DG and cause an error in the 4th step in Procedure 1. Hence, we must prune WG and DG firstly, deleting rows and columns where diiG=0. The method in Procedure 1 outputs one configuration of the airspace. We also use a chromosome to represent this solution and produce the initial generation by randomly mutating half genes of this individual several times. In this method, we need to set desired cluster number for the NSCA. However, this number is not the final counts of sectors because of possible change in the following evolution procedure. Hence, we directly set it as the count of sectors in current configuration. In following text, we call this caseNSCA Initialization, or NI.Procedure 1:Procedures of the method based on the NSCA.
(1)   Input DG, WG and desired cluster number, k. (2)   Prune DG and WG by deleting rows and columns where diiG=0. Define the pruned matrices asPDG and PWG. Define the set of deleted vertices as VD where VD=r and the set of remained vertices asRD. Define the dimension of PDG and PWG as PN=N-r. (3)   Calculate the Laplacian matrix, PLG=PDG-PWG. (4)   Compute the normalized Laplacian matrix, PLnor=IPN×PN-PDG-1/2PWGPDG-1/2. (5)   Compute the eigenvector corresponding to the smallest k eigenvalues of PLnor. (6)   Construct a matrix EPN×k, each of its column corresponds to the eigenvector of one of the smallestk eigenvalue. (7)   Normalize every row of EPN×k, meaning normalizing EPN×k into E~PN×k, E~i=Ei/∑keij21/2, where E~i and Ei stand for the ith row of E~ and E respectively, eij is the element at theith row and jth column of E. (8)   For each row of E~, set a corresponding vector yi  (i=1,2,…,PN). (9)   Cluster vectors y1,y2,…,yPN into k clusters using k-means cluster. The index of the cluster for each row is the index of the sector for theith vertex of G. (10) For any vertex in VD, assign it to the sector to which its nearest vertex in RD belongs.In the fourth type, half individuals in the initial generation are produced by CI and the other half is produced by NI. In following text, we call this caseMix Initialization, or MI.


### 4.3. Sector Boundary Depiction

As stated in Section1, the boundary depiction is much more challenging in graph cut-based approaches than that in geometric based and “cell gathering” based approaches. The methods directly using Voronoi polygons of vertices to form the boundaries are simple and straightforward. But they may not be very suitable in the sectorization problem for two reasons. Firstly, the shape of sectors is not strictly convex so that Voronoi polygons may not fit the outline of concave vertex cluster very well for its characteristic of convexity. Secondly, the MDFB constraint cannot be ensured by Voronoi polygons.We propose a concave hull-based method to depict the sector boundaries in this section. We firstly construct the concave hull for each vertex cluster (representing one sector) in Procedure2 with necessary explanations below.Procedure 2:Procedures to construct Concave(Bs).
(1) Initialize, input G=⋃s∈SGs, V=⋃s∈SVs, MP=⋃s∈SMPs, and OP=⋃sOPs (2) Fors∈S (3)   Set Bs=∅ (4)   Ifs is adjacent to the outer-boundary of the target airspace (5)    Obtain BMPs (See explain in Section 4.3.1), and set NBMs=∅ (6)    For each point pt∈BMPs (7)    Find point pto, that pto∈OP, and pto is nearest to pt (8)    NBMs=NBMs∪pto (9)    End (10)  Obtain OPs (See explain in Section 4.3.1) (11)  Bs=OPs∪MPs (12) End (13) Else (14)  Bs=MPs (15) End (16) Construct Concave(Bs) using α-shape algorithm (See explain in Section 4.3.1) (17) End


#### 4.3.1. Concave Hull Generation

Input Data. Input information of the concave hull generation method is in three categories.The first is the cut ofG: (19)G=⋃s=1kGs,V=⋃s=1kVs.The second is the set of midpoints of edges connecting two neighboring sectors,MP=mpi,j∣eij=1,vi∈Vsm,vj∈Vsn,sm≠sn. For any sector s, define the set of midpoints of edges connecting it and other sectors as MPs=mpi,j∣eij=1,vi∈Vs,vj∉Vs. Hence, we get that MP=⋃s∈SMPs and MPs1∩MPs2≠∅ if s1 and s2 are neighboring sectors.The third is a set of points located on the outerboundary of the target airspace. We suppose points on the outerboundary of the target airspace are distributed clockwise for the purpose of unification and define the point set asOP=opi∣1≤i≤n, OP=⋃sOPs, s∈S, and OPs=∅ if s is not next to the outerboundary.Obtain B M P s . For any sector s (subgraph Gs) adjacent to the outerboundary, we find out one or more couples of points in MPs, and denote the subset of MPs as BMPs. For any point mpi,j∈BMPs with corresponding edge eij, the following conditions must be satisfied:(i)
Both endpoints ofeij, denoted as vi and vj, are next to the outerboundary.(ii)
The Voronoi polygons ofvi and vj have one common side and this side intersects with the outerboundary.Note that, as the shape of any sector is not strictly convex, there may be more than one couple of points inBMPs. In Figure 5(a), we assume that V2 and V3 belong to sector s and V1 and V4 belong to the other two different sectors, respectively. In this situation, BMPs is made up of MP1 and MP2. In Figure 5(b), we assume that sector s consists of V2, V3, and V4, which surround V1. The shape of s is concave, but s still conforms to the convexity constraint because no air routes connect V3V1 or V4V1. In this situation, MP1, MP2, MP3, and MP4 are all elements of BMPs. BMPs with more than two couples of points is also possible which is not illustrated here. Note that BMPs1∩BMPs2≠∅, if s1 and s2 are neighboring sectors.Figure 5
Midpoints adjacent to the outerboundary.
(a) (b)Obtain O P s . For any sector s adjacent to the outerboundary, sort NBMs clockwise (as points in NBMs all belong to OP, sorting clockwise means sorting with increasing index in OP) as NBMs=op1s,op2s,…,op2k-1s,op2ks, where k is the number of couples nBMPs. Then, we define (20)OPs=⋃i=1kop2i-1s,…,op2is.We are aware that two neighboring points in NBMs may not be consecutive in OP. We define(21)op2i-1s,…,op2is=opind1,opind1+1,…,opind2-1,opind2,where ind1 is the index of op2i-1s in OP and ind2 is the index of op2is in OP.Construct Concave Hull of B s. We use the concave hull of Bs, denoted as ConcaveBs, to fit the outline of the vertex cluster of sector s. α-shape algorithm [38] is applied here to generate the concave hull. The basic idea of α-shape algorithm is to iteratively search α-extreme points and α-neighbors and to connect the found α-neighbors sequentially to generate α-hull. As large amount of extra definitions and notations must be introduced to formulate α-shape algorithm clearly, we do not present the details here.


#### 4.3.2. Concave Hull Amendment

C o n c a v e B s is a promising candidate to be the boundary of sector s because all vertices of Gs are inside it and it fits the shape of s very well. However, two more considerations make it necessary to recheck ConcaveBs. The first is the compactness of the boundary; that is, fewer segments are preferred on the boundary. The second is about the MDFB constraint; that is, the distance between some segments of ConcaveBs and some vertices of G may be less than the MDFB. Hence, we propose the Concave Hull Amendment Method (CHAM) to deal with these two considerations (Procedure 3).Procedure 3:Procedures of the CHAM.
(1) Fors∈S (2)   Get path(s) on Concave(Bs), set as P1s,P2s,…,Pkss (See Definition 4 in Section 4.3.2) (3)   For each path Pis, 1≤i≤ks (4)    Define vertices along Pis as V1p,V2p,…,VNip and corresponding BI1p,BI2p,…,BINip (note thatBIjp=Vjp when corresponding cross-sector edge has not been designed) (5)    Set j=1, gap=1. (6)    Whilej<Ni (7)     If cross-sector edges referred by BIjp and BIj+1p are both designed (8)      If line segment BIjpBIj+1p is infeasible (See Definition 5 in Section 4.3.2) (9)        Execute the OABD to BIjpBIj+1p (See Explain in Section 4.3.2) (10)    End (11)    j=j+1, go to line (6) (12)   End (13)   While line segment BIjpBIj+gapp is feasible (See Definition 5 in Section 4.3.2) (14)    gap=gap+1 (15)    If cross-sector edges referred by BIj+gap-1p and BIj+gapp are both designed (16)     break (17)    End (18)   End (19)   Ifgap=1 (20)    Execute the OABD to BIjpBIj+1p (See Explain in Section 4.3.2). (21)    Set the cross-sector edge referred by BIj+1p as designed, remain BIj+1p=Vj+1p (22)    j=j+1, go to line (6) (23)   End (24)   Else (25)    gap=gap-1 (26)    Draw line segment BIjpBIj+gapp, (27)    Set cross-sector edges referred by BIj+1p,BIj+2p,…,BIj+gapp as designed (28)    Set BIj+1p,BIj+2p,…,BIj+gapp as the intersections of related cross-sector edges and line segmentBIjpBIj+gapp (29)    j=j+gap, reset gap=1, go to line (6) (30)   End (31)  End (32) End (33) EndThere are two ideas in theCHAM. The first is to replace the polygonal line connecting several vertices on the concave hull with one line segment from the first vertex to the last one. The second idea is to avoid obstacle circles centering at the vertices of G with radius of that vertex (refer to Section 3.1 for the definition of the radius of a vertex). We call modified concave hulls after theCHAM boundary polygons and present two definitions about the concave hull and the boundary polygon first.Definition 4 (path(s) on the concave hull).
One scans the vertices ofConcaveBs clockwise. If s is adjacent to the outerboundary of the target airspace, whenever we encounter a sequence of vertices V1p-V2p-⋯-Vr-1p,Vrp, where V1p, Vrp∈OP or others belonging to MP, we define such a sequence as a path on the concave hull. If s is not adjacent to the outerboundary, that is, it is surrounded by other sectors, we simply define a path on ConcaveBs as a clockwise closed walk of all vertices of ConcaveBs.Definition 5.
A line segment connecting any two vertices on the boundary polygon of one sector is defined as feasible if it does not intersect with any inner-edges of any sector and it does not go through any obstacle circles centering at the vertices ofG.Next, we present another definition to avoid redepiction of common boundary between two neighboring sectors.Definition 6.
In theCHAM, once a segment on the boundary polygon is depicted, set all cross-sector edges that intersect with this segment as designed. Define the intersection of this segment and a cross-sector edge, eij, as the boundary crossing point of eij, denoted as BIij. When one cross-sector edge has not been designed, its default boundary crossing point is set as its midpoint, that is, the intersection of this edge and the concave hull. To ensure unification, when one endpoint of a depicted segment lies on the outerboundary of the target airspace, we also define a virtual cross-sector edge corresponding to this endpoint as designed and set the boundary crossing point as the endpoint itself.We also give two properties about the concave hull.Property 1.
For any sectors, line segments on ConcaveBs do not intersect with any inner-edges of any sector. This property is obvious from the definition and the generating process of α-shape [38].Property 2.
The distance between any vertex on the concave hull, representing the midpoint of a cross-sector edge, and any vertex ofG is at least MDFB.Proof.
For any vertex on the concave hull, if the corresponding cross-sector edge represents a real air-route segment, the edge is at least MDFB away from any vertex ofG according to Rule 3 in Section 3.2 except its two endpoints. In addition, as all edges of G are at least twice as long as MDFB, their midpoints are at least MDFB away from the two endpoints. For one vertex on the concave hull that is the midpoint of a sector-cross edge only standing for Voronoi neighboring relationship, it locates on the Voronoi polygons of both endpoints of this edge, which means that the nearest vertices of G to this midpoint are just the two endpoints of the edge. Considering the possible distance between this midpoint and the two endpoints is at least MDFB, if any other vertex of G is less than MDFB from this midpoint, a contradictory occurs.These two properties show two pieces of important information when we construct the boundary polygons:(1)
For any two consecutive vertices on the concave hull of a sector, if we can find a line segment or a polygonal line that does not go through any obstacle circles, this line segment or polygonal line can be set as boundary segment(s) of that sector.(2)
The boundary polygons after theCHAM can cross any vertices on the concave hulls.We show the procedures of theCHAM in Table 1 with necessary explanations below.Table 1
Results of six sectorization methods.
Time period Sectorization method Obj1 Obj2 Obj3 Con1 Con2 Con3 Nos Density gap RT (s) 12:00–14:00 Current in practice 0.8145 674 238 0 1 0 16 4218 — Only NSCA 1.0518 471 133 0 6 1 16 6027 — RI 0.3158 404 84 1 0 0 9 2521 934 CI 0.3124 360 59 0 0 0 9 2232 940 NI 0.3404 375 58 0 0 0 8 2356 942 MI 0.4329 386 68 0 0 0 8 2512 964  16:00–18:00 Current in practice 1.1466 707 234 0 4 2 16 6705 — Only NSCA 1.6230 428 86 0 0 3 16 6732 — RI 0.3148 525 106 1 0 0 10 2409 926 CI 0.4508 608 130 0 0 0 10 3007 937 NI 0.4259 499 104 0 0 0 10 2521 924 MI 0.4502 511 91 0 0 0 10 2619 940  20:00–22:00 Current in practice 1.3666 710 247 0 3 2 16 7198 — Only NSCA 2.0712 403 72 0 1 5 16 9150 — RI 0.4514 684 206 2 0 0 12 2375 936 CI 0.3487 719 218 0 0 0 12 1902 934 NI 0.4142 592 134 0 0 0 12 2258 930 MI 0.3993 598 136 0 0 0 12 2017 915  2:00–24:00 Current in practice 1.3850 427 166 0 1 2 16 7604 — Only NSCA 1.6230 428 86 1 2 3 16 9058 — RI 0.4318 468 205 2 0 0 12 2635 944 CI 0.4100 335 79 0 0 0 12 2239 937 NI 0.3678 350 105 0 0 0 12 2159 924 MI 0.2982 351 118 0 0 0 12 2073 930Obstacle Avoidance Boundary Design (OABD). We illustrate theOABD with Figure 6. Firstly, we assume two consecutive vertices on the concave hull as Vs and Ve. VsVe is assumed to intrude one or more obstacle circles (Figure 6(a)). TheOABD examines the first obstacle (such as Ob1 in Figure 6) intruded by VsVe. Then, we draw two tangents from Vs and Ve to the obstacle circle, respectively. The intersection of these two tangents is denoted as V0. Thus, the boundary between Vs and Ve is now VsV0Ve (Figure 6(b)). Next, we examine if VsV0Ve still intrudes any other obstacle (such as Ob2 in Figure 6). If it does, we truncate V0Ve at the boundary of the second obstacle (V1 in Figure 6(c)). Then, we draw two tangents of Ob2 originating from V1 and Ve, respectively, and get their intersection (V2 in Figure 6(c)). The boundary between Vs and Ve is now changed to VsV0V1V2Ve. TheOABD continues until the curve VsV0Ve (1 iteration) or VsV0V1⋯Vk+1Ve (k iterations, k>1) is free of all obstacles (Figure 6(c)). If too many iterations are needed in theOABD, the generated sector boundaries may be jagged, which is prohibited. Throughout our simulation, such situations do not appear. Relevant discussions will be in Section 6.5.Figure 6
Procedures of theOABD.
(a) (b) (c)In theOABD, there are always two intersections when we draw tangents, such as I1 and I2, in Figure 7. As I1 and I2 always locate symmetrically on opposite sides of the obstacle circle, one of the generated polygonal lines definitely intersects with inner-edges of sectors while the other does not. The one without intersections with inner-edges is our choice (such as VsI2Ve in Figure 7).Figure 7
Selection when drawing tangents.


#### 4.3.3. Vacuum Space Handle

After theCHAM, there still remains several vacuum spaces at the intersecting areas of three or more sectors, such as quadrilateralABCD in Figure 8(a). These vacuum spaces contain no vertices or edges and are not assigned to any sector. In this section, we propose a heuristic method to assign vacuum spaces. Firstly, we represent a vacuum space with a polygon, calling the vacuum polygon. The vacuum polygons can be sorted into four categories according to types of their sides:(I)
Triangles with no side from theOABD.(II)
Polygons with more than three sides and no side that is from theOABD.(III)
Quadrangles with two sides from theOABD.(IV)
Other types.Figure 8
Vacuum polygon handle, Category (II).
(a) (b) (c)For vacuum polygons in Category (I), we propose two choices, one is directly erasing the opposite side of the biggest inner angle (Figure9(a)), and the other is to find the “most vertical” midline and erase those two sides on which the midpoint does not locate (Figure 9(b)).Figure 9
Vacuum polygon handle, Categories (I) and (III).
(a) (b) (c)For vacuum polygons in Category (II), we first give the following property.Property 3.
The vacuum polygon in Category (II) is convex.Proof.
Suppose there is a concave vacuum polygon in Category (II) (such as the one in Figure8(b)). We consider the vertex causing the concavity (point C in Figure 8(b)). As this vertex is on a cross-sector edge, one of the endpoints of this cross-sector edge locates inside the vacuum polygon because of the concavity, obviously causing a contradictory against the definition of the vacuum space. In another situation, it locates outside the polygon while the cross-sector edge must intersect with one extra side of the polygon. However, neither of the two endpoints of this extra side is the vertex that causes the concavity; that is, a cross-sector edge between two sectors intersects with the boundary of an irrelevant sector. It is also impossible. Hence, we show that the vacuum polygons in Category (II) must be convex.With Property3, we can calculate the geometric center of a polygon in Category (II), which must locate inside the polygon because of its convexity (point O in Figure 8(b)). Then, we connect the center to every vertices of the polygon and erase original sides of the polygon (Figure 8(c)).For polygons in Category (III), which stand for an intersecting area of three sectors, we find the vertex produced by the OABD (pointI in Figure 9(c)), draw the diagonal from it, and erase those two opposite sides of the vertex produced by the OABD.For polygons in Category (IV), situation of their sides may be very complex. We directly assign the entire vacuum space to one sector, whose boundary segments are sides of the polygon afterOABD.


## 4.1. Mathematical Expression

As stated before, the airspace is modeled as an undirected graph,G. Solving the 2D ERSP is to cut G into k subgraphs, GsVs,Es, where s∈S=1,2,…,k, (7)G=⋃s=1kGs,V=⋃s=1kVs,Vi∩Vj=∅,i,j∈S,i≠j.Each subgraph stands for one sector.


### 4.1.1. Objectives

(1) Balancing the Intrasector Workload. Consider (8) min ⁡ 1 k ∑ s = 1 k W L G s - W L - 2 . W L G s is the intraworkload of sector s and WL- is the mean value of WLGs, s∈S. WLGs is a weighted and scaled summation of the values of four metrics of sector s, including information about the flow density, the altitude change, the speed change, and the course change: (9)WLGs=wFDFD~s+wacAC~s+wscSC~s+wccCC~s,where wFD, wac, wsc, and wcc are weights of these four parts.We defineEs′ as the set of intrasector edges of sector s, where Es′=eij∣eij=1,vi,vj∈Vs and Es′′ are the set of intersector edges between s and other sectors, where Es′′=eij∣eij=1,vi∈Vs,vj∉Vs; then(10)FDs=∑eij∈Es′STPij+∑eij∈Es′′ATPij,ACs=∑eij∈Es′SACij+∑eij∈Es′′AACijFDs,SCs=∑eij∈Es′′SSCij+∑eij∈Es′′ASCijFDs,CCs=NTsNFs,where NTs is the total course change times of flights when flying inside sector s and NFs is the total number of flights flying through s. Hence, CCs stands for average course change times of single flight in sector s. Suppose the flight set considered in the sectorization problem as F, and we have(11)NTs=∑f∈Fvi∈VsFCCfi.Because of different ranges of the four parts in (9), and a more meaningful understanding of the summation, we normalize them to the range of [0,1]. We set four standard values, FDsd, ACsd, SCsd, and CCsd, that are large enough and do not change during the optimization. Then, we define(12)FD~s=FDsFDsd,AC~s=ACsACsd,SC~s=SCsSCsd,CC~s=CCsCCsd.(2) Minimizing the Aggregate Coordination Workload. We use the density of intersector flow to stand for the coordination workload. This objective can be expressed by(13)min⁡∑s≠rs,r∈SWLcoGs,Gr,where WLco(Gs,Gr) is the intersector flow density between Gs and Gr. It is the number of flights flying across Gs and Gr:(14)WLcoGs,Gr=∑vp∈Vsvq∈VrFDpq.(3) Minimizing the Occurrences of the Violations of the Minimum Staying Time Constraint. If any flight stays in the sector less than a minimum time threshold after it enters a sector, the minimum staying time constraint is violated.Obj3 is to minimize the occurrences of these situations.As we do not depict exact sector boundaries during the optimization process, it is impossible to accurately calculate the staying time of any flight. We approximate the accurate staying time as follows.Still suppose the flight path of flightf as fpf=vf,1,vf,2,…,vf,k, vf,1,vf,2,…,vf,k∈V, where vf,i is the ith vertex along the path. Then, we calculate the staying time of flight f in sector sas follows: (15)dtfs=12VATf,a-VLTf,a-1+12VATf,e+1-VLTf,e+VLTf,e-VATf,e+∑vf,m∈Vsvf,m+1∈VsVATf,m+1-VATf,m,where va and ve are the first and the last vertices inside sector s along fpf. The first term of (15) stands for the time spent to reach va after entering the sector. The second term stands for the time spent to leave the sector after leaving ve. The third term stands for the time spent in ve. The fourth term stands for other times spent inside sector s.We define the set of all considered flights asF and set mdtvfs=1, if any flight f∈F violates the minimum staying time constraint in sector s. Hence,Obj3 can be expressed as(16)min⁡∑f∈F∑s∈Smdtvfs.


### 4.1.2. Constraints

(1) Connectivity Constraint. Each subgraph representing one sector must be self-connected.(2) Flight-Route-Based Convexity Constraint. For any flight f, fpf=vf,1,vf,2,…,vf,k, vf,1,vf,2,…,vf,k∈V. If vf,i,vf,j∈Vs, i≤j, for ∀vf,m, i≤m≤j, vf,m∈Vs.(3) Instantaneous Flow Constraint. Similar to the expressions above, we suppose va and ve as the first and last vertex inside sector s along fpf. Hence, the entry time of sector s can be denoted as mid(VLTa-1,VATa), which is the midtime of the time leaving va-1 and the time reaching va. Similarly, we denote the departure time of sector s as mid(VLTe,VATe+1). With the entry time and the departure time for every flight, we can calculate the number of flights in any sector at any moment. Suppose the maximum instantaneous flight number in sector s as Instflows and the instantaneous flow threshold as Threinsflow, andCon3 can be expressed as(17)Instflowf,s≤Threinsflow,∀s.Con4 andCon5 are not considered in the optimization phase. But they are ensured by the airspace modeling in Section 3 and the boundary depiction method in Section 4.3.


## 4.1.1. Objectives

(1) Balancing the Intrasector Workload. Consider (8) min ⁡ 1 k ∑ s = 1 k W L G s - W L - 2 . W L G s is the intraworkload of sector s and WL- is the mean value of WLGs, s∈S. WLGs is a weighted and scaled summation of the values of four metrics of sector s, including information about the flow density, the altitude change, the speed change, and the course change: (9)WLGs=wFDFD~s+wacAC~s+wscSC~s+wccCC~s,where wFD, wac, wsc, and wcc are weights of these four parts.We defineEs′ as the set of intrasector edges of sector s, where Es′=eij∣eij=1,vi,vj∈Vs and Es′′ are the set of intersector edges between s and other sectors, where Es′′=eij∣eij=1,vi∈Vs,vj∉Vs; then(10)FDs=∑eij∈Es′STPij+∑eij∈Es′′ATPij,ACs=∑eij∈Es′SACij+∑eij∈Es′′AACijFDs,SCs=∑eij∈Es′′SSCij+∑eij∈Es′′ASCijFDs,CCs=NTsNFs,where NTs is the total course change times of flights when flying inside sector s and NFs is the total number of flights flying through s. Hence, CCs stands for average course change times of single flight in sector s. Suppose the flight set considered in the sectorization problem as F, and we have(11)NTs=∑f∈Fvi∈VsFCCfi.Because of different ranges of the four parts in (9), and a more meaningful understanding of the summation, we normalize them to the range of [0,1]. We set four standard values, FDsd, ACsd, SCsd, and CCsd, that are large enough and do not change during the optimization. Then, we define(12)FD~s=FDsFDsd,AC~s=ACsACsd,SC~s=SCsSCsd,CC~s=CCsCCsd.(2) Minimizing the Aggregate Coordination Workload. We use the density of intersector flow to stand for the coordination workload. This objective can be expressed by(13)min⁡∑s≠rs,r∈SWLcoGs,Gr,where WLco(Gs,Gr) is the intersector flow density between Gs and Gr. It is the number of flights flying across Gs and Gr:(14)WLcoGs,Gr=∑vp∈Vsvq∈VrFDpq.(3) Minimizing the Occurrences of the Violations of the Minimum Staying Time Constraint. If any flight stays in the sector less than a minimum time threshold after it enters a sector, the minimum staying time constraint is violated.Obj3 is to minimize the occurrences of these situations.As we do not depict exact sector boundaries during the optimization process, it is impossible to accurately calculate the staying time of any flight. We approximate the accurate staying time as follows.Still suppose the flight path of flightf as fpf=vf,1,vf,2,…,vf,k, vf,1,vf,2,…,vf,k∈V, where vf,i is the ith vertex along the path. Then, we calculate the staying time of flight f in sector sas follows: (15)dtfs=12VATf,a-VLTf,a-1+12VATf,e+1-VLTf,e+VLTf,e-VATf,e+∑vf,m∈Vsvf,m+1∈VsVATf,m+1-VATf,m,where va and ve are the first and the last vertices inside sector s along fpf. The first term of (15) stands for the time spent to reach va after entering the sector. The second term stands for the time spent to leave the sector after leaving ve. The third term stands for the time spent in ve. The fourth term stands for other times spent inside sector s.We define the set of all considered flights asF and set mdtvfs=1, if any flight f∈F violates the minimum staying time constraint in sector s. Hence,Obj3 can be expressed as(16)min⁡∑f∈F∑s∈Smdtvfs.


## 4.1.2. Constraints

(1) Connectivity Constraint. Each subgraph representing one sector must be self-connected.(2) Flight-Route-Based Convexity Constraint. For any flight f, fpf=vf,1,vf,2,…,vf,k, vf,1,vf,2,…,vf,k∈V. If vf,i,vf,j∈Vs, i≤j, for ∀vf,m, i≤m≤j, vf,m∈Vs.(3) Instantaneous Flow Constraint. Similar to the expressions above, we suppose va and ve as the first and last vertex inside sector s along fpf. Hence, the entry time of sector s can be denoted as mid(VLTa-1,VATa), which is the midtime of the time leaving va-1 and the time reaching va. Similarly, we denote the departure time of sector s as mid(VLTe,VATe+1). With the entry time and the departure time for every flight, we can calculate the number of flights in any sector at any moment. Suppose the maximum instantaneous flight number in sector s as Instflows and the instantaneous flow threshold as Threinsflow, andCon3 can be expressed as(17)Instflowf,s≤Threinsflow,∀s.Con4 andCon5 are not considered in the optimization phase. But they are ensured by the airspace modeling in Section 3 and the boundary depiction method in Section 4.3.


## 4.2. Constrained NSGA-II Based Evolution

Some objectives and constraints of the 2D ERSP cannot be written in explicit mathematical formulas [2, 24]. For example, consider the constraint that no flights can enter any single sector more than once. It is affected by the trajectory of flights and the boundary of a sector, which is too complex to express this constraint by an explicit inequality. Hence, it is impossible to construct the 2D ERSP into a typical form of mathematical programming below:(18)Min⁡fxAx≤bx∈X.As evolution algorithms can still work even if the constraints can only be expressed by text descriptions, they are good choices to solve the 2D ERSP. The Constrained Nondominated Sorting Genetic Algorithm II (CNSGA-II) [26] is a fast and elitist multiobjective genetic algorithm that is very appropriate for the 2D ERSP, considering its characteristics of multiobjectives and the strong capability in handling nonclassic constraints. We utilize the CNSGA-II accompanied with the memetic local search [27] to solve the 2D ERSP. The main flow diagram of the evolution algorithm is shown in Figure 3.Figure 3
Diagram of the evolution algorithm.


### 4.2.1. Chromosome

Every chromosome hasN genes representing N vertices in graph G. The value of each gene is the index of the sector to which the related vertex belongs. For example, in Figure 4, the third vertex is assigned to sector 2 and the 50th vertex is assigned to sector 5.Figure 4
Chromosome.We set a common largest possible index of sectors for every individual; that is, the possible maximum number of sectors is fixed. On one hand, this number is large enough so that it has no impact on the results. On the other hand, because we only set the possible maximum number, indexes of sectors in different individuals are totally independent, which means that two sectors with the same index in different individuals may refer to different sectors.


### 4.2.2. Reproduction

Tournament selection is used to select individuals in parent generation.


### 4.2.3. Crossover

Two individuals can crossover only when they have large enough hamming distance. In single crossover operation, two randomly chosen genes are exchanged. To avoid generating individuals that violates connectivity constraint as much as possible, the crossover can be granted only when at least one of the vertices (genes) is moved to an adjacent sector or stays in original sector after the crossover.


### 4.2.4. Mutation

In single mutation operation, one gene is randomly selected. The value of the gene is randomly changed and it is bounded by the value of the maximum sector number. However, similar to the crossover operation, only when the related vertex is changed to an adjacent sector, the mutation operation can be granted.


### 4.2.5. Memetic Local Search

The memetic algorithm [27] is a local search method embedded in one global search iteration of evolution algorithms that can help improving the individuals in one generation. In this paper, two heuristic polices of memetic local search are applied.Policy A. For any sector, if FDs>FDup, where FDup is an upper bound, we split the sector into two or more parts. For the subgraph Gs corresponding to sector s, we build its similarity matrix WNs×NsGs, where Ns is the number of vertices that belong to Gs. For any element wijGs, 1≤i≤Ns, 1≤j≤Ns, suppose the indices of its two related vertices in graph G as p and q, and then wijGs=FDpq. We add up all elements of each row of WNs×NsGs to build a diagonal matrix DGs, where diiGs=∑j=1NswijGs. We define theLaplacian matrix of Gs, LGs=DGs-WGs. Besides, we introduce three lemmas in the Spectral Graph Theory [33].Lemma 1.
For a weighed undirected graph with nonnegative edge weights, itsLaplacian matrix has at least one eigenvalue with value of 0.Lemma 2.
The number of isolated subgraphs of a weighed undirected graph with nonnegative edge weights equals the number of 0 eigenvalues of itsLaplacian matrix.Lemma 3.
For a weighed undirected graph with nonnegative edge weights, denoted asG, suppose the number of 0 eigenvalues of itsLaplacian matrix as k. Then, the dimension of the eigenspace of 0 eigenvalues is k. This linear subspace is spanned by such a group of base vectors: 1s1,1s2,…,1sk, where 1si(1≤i≤k) is the indicator vector of ith isolated subgraph, denoted as Gi. If a vertex is in Gi, the corresponding element in 1si is 1; otherwise, it is 0.With these three lemmas, we calculate the eigenvalues ofLGs and the corresponding eigenvectors. If the multiplicity of 0 eigenvalues, k, is larger than 1, we just divide Gs into k parts according to the base vectors. If k=1, we get the eigenvector corresponding to the second smallest eigenvalue, calledFiedler Vector [34]. Then, the vertices of Gs are clustered into two groups. One is with the positive elements ofFiedler Vector and the other is with the negative ones. Vertices corresponding to the elements with value of 0 are randomly assigned to one of these two groups. Hence, when k=1, we split Gs into two parts.Policy B. Policy B is operated after the division by Policy A. For any individual, if FDs<FDlow for sector s, where FDlow is a lower bound, we treat the flow density in this sector as too low. Thus, we merge s into one of its adjacent sectors, assumed as s′. Values of genes (vertices) belonging to s are changed to the index of s′.To avoid loss of diversity and being trapped in regional optima, the memetic local search is operated with gradually changing intensity and frequency [35]. The frequency refers to the gap of generations between two callings of the local search, and the intensity refers to the proportion of individuals in one single generation that are applied with the local search. At the beginning of the evolution, the qualities of the individuals are still low, and, in general, we set low frequency and intensity. Along with the evolution, the searching space approaches the optimal solutions, and thus we need more accurate search. Therefore, the frequency and intensity of the local search both increase.


### 4.2.6. Fitness Functions

Fitness functions correspond to the expressions in formulas (8), (13), and (16).


### 4.2.7. Constraint Check

Connectivity. Similar to those stated in “Policy A” of the memetic local search, for each subgraph Gs (sectors), we calculate itsLaplacian matrix LGs and check the multiplicity of 0 eigenvalues. If the number is larger than 1, we record one more time of violation of the connectivity constraint for current individual.Flight-Route-Based-Convexity. Still suppose the flight path of flight f as vf,1,vf,2,…,vf,k. If vf,p∈Gs and vf,p+1∉Gs, we check if vf,q∈Gs, q>p+1. Once the condition vf,q∈Gs, q>p+1 is met, we record one more time of the violation of the convexity constraint for current individual.Instantaneous Flow Constraint. We check the validation of inequality (17). If inequality (17) is unsatisfied, we record one more time of the violation of the instantaneous flow constraint for current individual.


### 4.2.8. Sort and Elitism

In the CNSGA-II, the elitism operation is executed for every generation. The dominating relationship between two individuals (solutions),i and j, is defined below:(1)
i≺j, if i is feasible (no constraint violation) while j is not.(2)
When solutionsi and j are both feasible, i≺j if the values of all fitness functions of i are not more than that of j and at least one of such values is smaller than that of j.(3)
When solutionsi and j are both infeasible, we compare their times of constraint violations in the following order of decreasing priority: connectivity > convexity > instantaneous flow. That is, i≺j, if the connectivity constraint violating times of solution i is less than that of solution j. Or, if the connectivity constraint violating times are equal, we continue to compare the violating times of the convexity constraint and so on.(4)
When solutioni and j are both infeasible and they are equivalent in terms of the violating times of all three constraints, i≺j if the values of all fitness functions of i are not more than that of j and at least one such value is smaller than that of j.


### 4.2.9. Initial Generation

Most previous studies using evolutionary algorithms always concern about how to refine the “genetic operations” (such as the crossover or mutation rules) or try to propose heuristics to improve the evolution procedure. However, for those problems that the gross optima or a boundary of the gross optima cannot be preliminarily determined (such as the 2D ERSP), no matter how excellent a heuristic or an improved operation seems to be, we are unable to determine whether the gross optima can be found with the help of these heuristics. In these problems, we cannot guarantee reaching the gross optima from any type of initial generations. Hence, we have to realize that different initial generations have different impact on final solutions. However, discussions about the initial generation of evolution algorithms and its impact to the final solutions are absent in relevant literatures.In this paper, we propose four types of initial generations and compare them with the corresponding nondominated solution sets produced by the CNSGA-II.The first type is the most general choice that just randomly produces the initial generation. In following text, we call this caseRandomly Initialization, or RI.The second type takes advantage of the current configuration in practice. We represent the current configuration with a chromosome and randomly mutate half of the genes several times to produce the initial generation. In following text, we call this caseCurrent Initialization, or CI.Thirdly, we employ a method based on the normalized spectral clustering algorithm (NSCA) [36]. The NSCA is a data clustering algorithm that shows its talent in the graph N-cut problem if each data point is represented by one vertex in a weighted graph and the weights of edges stand for similarities between two data points [37]. We choose the NSCA because of two objectives of the 2D ERSP, balancing intrasector workload and minimizing gross intersector workload, which are very similar to the goals of the undirected graph N-cut problem. We first change the values of elements in matrix FD with value of ε to 0 and call the modified density matrix FD′. Then, FD′ is defined as the similarity matrix of G; that is, WG=FD′ and we construct a diagonal matrix DG, where diiG=∑j=1NwijG. Following steps of the method based on the NSCA are listed in Procedure 1. Note that diiG may equal to 0 if no aircraft flies through any routes connecting to vertex vi. This may lead singularity of DG and cause an error in the 4th step in Procedure 1. Hence, we must prune WG and DG firstly, deleting rows and columns where diiG=0. The method in Procedure 1 outputs one configuration of the airspace. We also use a chromosome to represent this solution and produce the initial generation by randomly mutating half genes of this individual several times. In this method, we need to set desired cluster number for the NSCA. However, this number is not the final counts of sectors because of possible change in the following evolution procedure. Hence, we directly set it as the count of sectors in current configuration. In following text, we call this caseNSCA Initialization, or NI.Procedure 1:Procedures of the method based on the NSCA.
(1)   Input DG, WG and desired cluster number, k. (2)   Prune DG and WG by deleting rows and columns where diiG=0. Define the pruned matrices asPDG and PWG. Define the set of deleted vertices as VD where VD=r and the set of remained vertices asRD. Define the dimension of PDG and PWG as PN=N-r. (3)   Calculate the Laplacian matrix, PLG=PDG-PWG. (4)   Compute the normalized Laplacian matrix, PLnor=IPN×PN-PDG-1/2PWGPDG-1/2. (5)   Compute the eigenvector corresponding to the smallest k eigenvalues of PLnor. (6)   Construct a matrix EPN×k, each of its column corresponds to the eigenvector of one of the smallestk eigenvalue. (7)   Normalize every row of EPN×k, meaning normalizing EPN×k into E~PN×k, E~i=Ei/∑keij21/2, where E~i and Ei stand for the ith row of E~ and E respectively, eij is the element at theith row and jth column of E. (8)   For each row of E~, set a corresponding vector yi  (i=1,2,…,PN). (9)   Cluster vectors y1,y2,…,yPN into k clusters using k-means cluster. The index of the cluster for each row is the index of the sector for theith vertex of G. (10) For any vertex in VD, assign it to the sector to which its nearest vertex in RD belongs.In the fourth type, half individuals in the initial generation are produced by CI and the other half is produced by NI. In following text, we call this caseMix Initialization, or MI.


## 4.2.1. Chromosome

Every chromosome hasN genes representing N vertices in graph G. The value of each gene is the index of the sector to which the related vertex belongs. For example, in Figure 4, the third vertex is assigned to sector 2 and the 50th vertex is assigned to sector 5.Figure 4
Chromosome.We set a common largest possible index of sectors for every individual; that is, the possible maximum number of sectors is fixed. On one hand, this number is large enough so that it has no impact on the results. On the other hand, because we only set the possible maximum number, indexes of sectors in different individuals are totally independent, which means that two sectors with the same index in different individuals may refer to different sectors.


## 4.2.2. Reproduction

Tournament selection is used to select individuals in parent generation.


## 4.2.3. Crossover

Two individuals can crossover only when they have large enough hamming distance. In single crossover operation, two randomly chosen genes are exchanged. To avoid generating individuals that violates connectivity constraint as much as possible, the crossover can be granted only when at least one of the vertices (genes) is moved to an adjacent sector or stays in original sector after the crossover.


## 4.2.4. Mutation

In single mutation operation, one gene is randomly selected. The value of the gene is randomly changed and it is bounded by the value of the maximum sector number. However, similar to the crossover operation, only when the related vertex is changed to an adjacent sector, the mutation operation can be granted.


## 4.2.5. Memetic Local Search

The memetic algorithm [27] is a local search method embedded in one global search iteration of evolution algorithms that can help improving the individuals in one generation. In this paper, two heuristic polices of memetic local search are applied.Policy A. For any sector, if FDs>FDup, where FDup is an upper bound, we split the sector into two or more parts. For the subgraph Gs corresponding to sector s, we build its similarity matrix WNs×NsGs, where Ns is the number of vertices that belong to Gs. For any element wijGs, 1≤i≤Ns, 1≤j≤Ns, suppose the indices of its two related vertices in graph G as p and q, and then wijGs=FDpq. We add up all elements of each row of WNs×NsGs to build a diagonal matrix DGs, where diiGs=∑j=1NswijGs. We define theLaplacian matrix of Gs, LGs=DGs-WGs. Besides, we introduce three lemmas in the Spectral Graph Theory [33].Lemma 1.
For a weighed undirected graph with nonnegative edge weights, itsLaplacian matrix has at least one eigenvalue with value of 0.Lemma 2.
The number of isolated subgraphs of a weighed undirected graph with nonnegative edge weights equals the number of 0 eigenvalues of itsLaplacian matrix.Lemma 3.
For a weighed undirected graph with nonnegative edge weights, denoted asG, suppose the number of 0 eigenvalues of itsLaplacian matrix as k. Then, the dimension of the eigenspace of 0 eigenvalues is k. This linear subspace is spanned by such a group of base vectors: 1s1,1s2,…,1sk, where 1si(1≤i≤k) is the indicator vector of ith isolated subgraph, denoted as Gi. If a vertex is in Gi, the corresponding element in 1si is 1; otherwise, it is 0.With these three lemmas, we calculate the eigenvalues ofLGs and the corresponding eigenvectors. If the multiplicity of 0 eigenvalues, k, is larger than 1, we just divide Gs into k parts according to the base vectors. If k=1, we get the eigenvector corresponding to the second smallest eigenvalue, calledFiedler Vector [34]. Then, the vertices of Gs are clustered into two groups. One is with the positive elements ofFiedler Vector and the other is with the negative ones. Vertices corresponding to the elements with value of 0 are randomly assigned to one of these two groups. Hence, when k=1, we split Gs into two parts.Policy B. Policy B is operated after the division by Policy A. For any individual, if FDs<FDlow for sector s, where FDlow is a lower bound, we treat the flow density in this sector as too low. Thus, we merge s into one of its adjacent sectors, assumed as s′. Values of genes (vertices) belonging to s are changed to the index of s′.To avoid loss of diversity and being trapped in regional optima, the memetic local search is operated with gradually changing intensity and frequency [35]. The frequency refers to the gap of generations between two callings of the local search, and the intensity refers to the proportion of individuals in one single generation that are applied with the local search. At the beginning of the evolution, the qualities of the individuals are still low, and, in general, we set low frequency and intensity. Along with the evolution, the searching space approaches the optimal solutions, and thus we need more accurate search. Therefore, the frequency and intensity of the local search both increase.


## 4.2.6. Fitness Functions

Fitness functions correspond to the expressions in formulas (8), (13), and (16).


## 4.2.7. Constraint Check

Connectivity. Similar to those stated in “Policy A” of the memetic local search, for each subgraph Gs (sectors), we calculate itsLaplacian matrix LGs and check the multiplicity of 0 eigenvalues. If the number is larger than 1, we record one more time of violation of the connectivity constraint for current individual.Flight-Route-Based-Convexity. Still suppose the flight path of flight f as vf,1,vf,2,…,vf,k. If vf,p∈Gs and vf,p+1∉Gs, we check if vf,q∈Gs, q>p+1. Once the condition vf,q∈Gs, q>p+1 is met, we record one more time of the violation of the convexity constraint for current individual.Instantaneous Flow Constraint. We check the validation of inequality (17). If inequality (17) is unsatisfied, we record one more time of the violation of the instantaneous flow constraint for current individual.


## 4.2.8. Sort and Elitism

In the CNSGA-II, the elitism operation is executed for every generation. The dominating relationship between two individuals (solutions),i and j, is defined below:(1)
i≺j, if i is feasible (no constraint violation) while j is not.(2)
When solutionsi and j are both feasible, i≺j if the values of all fitness functions of i are not more than that of j and at least one of such values is smaller than that of j.(3)
When solutionsi and j are both infeasible, we compare their times of constraint violations in the following order of decreasing priority: connectivity > convexity > instantaneous flow. That is, i≺j, if the connectivity constraint violating times of solution i is less than that of solution j. Or, if the connectivity constraint violating times are equal, we continue to compare the violating times of the convexity constraint and so on.(4)
When solutioni and j are both infeasible and they are equivalent in terms of the violating times of all three constraints, i≺j if the values of all fitness functions of i are not more than that of j and at least one such value is smaller than that of j.


## 4.2.9. Initial Generation

Most previous studies using evolutionary algorithms always concern about how to refine the “genetic operations” (such as the crossover or mutation rules) or try to propose heuristics to improve the evolution procedure. However, for those problems that the gross optima or a boundary of the gross optima cannot be preliminarily determined (such as the 2D ERSP), no matter how excellent a heuristic or an improved operation seems to be, we are unable to determine whether the gross optima can be found with the help of these heuristics. In these problems, we cannot guarantee reaching the gross optima from any type of initial generations. Hence, we have to realize that different initial generations have different impact on final solutions. However, discussions about the initial generation of evolution algorithms and its impact to the final solutions are absent in relevant literatures.In this paper, we propose four types of initial generations and compare them with the corresponding nondominated solution sets produced by the CNSGA-II.The first type is the most general choice that just randomly produces the initial generation. In following text, we call this caseRandomly Initialization, or RI.The second type takes advantage of the current configuration in practice. We represent the current configuration with a chromosome and randomly mutate half of the genes several times to produce the initial generation. In following text, we call this caseCurrent Initialization, or CI.Thirdly, we employ a method based on the normalized spectral clustering algorithm (NSCA) [36]. The NSCA is a data clustering algorithm that shows its talent in the graph N-cut problem if each data point is represented by one vertex in a weighted graph and the weights of edges stand for similarities between two data points [37]. We choose the NSCA because of two objectives of the 2D ERSP, balancing intrasector workload and minimizing gross intersector workload, which are very similar to the goals of the undirected graph N-cut problem. We first change the values of elements in matrix FD with value of ε to 0 and call the modified density matrix FD′. Then, FD′ is defined as the similarity matrix of G; that is, WG=FD′ and we construct a diagonal matrix DG, where diiG=∑j=1NwijG. Following steps of the method based on the NSCA are listed in Procedure 1. Note that diiG may equal to 0 if no aircraft flies through any routes connecting to vertex vi. This may lead singularity of DG and cause an error in the 4th step in Procedure 1. Hence, we must prune WG and DG firstly, deleting rows and columns where diiG=0. The method in Procedure 1 outputs one configuration of the airspace. We also use a chromosome to represent this solution and produce the initial generation by randomly mutating half genes of this individual several times. In this method, we need to set desired cluster number for the NSCA. However, this number is not the final counts of sectors because of possible change in the following evolution procedure. Hence, we directly set it as the count of sectors in current configuration. In following text, we call this caseNSCA Initialization, or NI.Procedure 1:Procedures of the method based on the NSCA.
(1)   Input DG, WG and desired cluster number, k. (2)   Prune DG and WG by deleting rows and columns where diiG=0. Define the pruned matrices asPDG and PWG. Define the set of deleted vertices as VD where VD=r and the set of remained vertices asRD. Define the dimension of PDG and PWG as PN=N-r. (3)   Calculate the Laplacian matrix, PLG=PDG-PWG. (4)   Compute the normalized Laplacian matrix, PLnor=IPN×PN-PDG-1/2PWGPDG-1/2. (5)   Compute the eigenvector corresponding to the smallest k eigenvalues of PLnor. (6)   Construct a matrix EPN×k, each of its column corresponds to the eigenvector of one of the smallestk eigenvalue. (7)   Normalize every row of EPN×k, meaning normalizing EPN×k into E~PN×k, E~i=Ei/∑keij21/2, where E~i and Ei stand for the ith row of E~ and E respectively, eij is the element at theith row and jth column of E. (8)   For each row of E~, set a corresponding vector yi  (i=1,2,…,PN). (9)   Cluster vectors y1,y2,…,yPN into k clusters using k-means cluster. The index of the cluster for each row is the index of the sector for theith vertex of G. (10) For any vertex in VD, assign it to the sector to which its nearest vertex in RD belongs.In the fourth type, half individuals in the initial generation are produced by CI and the other half is produced by NI. In following text, we call this caseMix Initialization, or MI.


## 4.3. Sector Boundary Depiction

As stated in Section1, the boundary depiction is much more challenging in graph cut-based approaches than that in geometric based and “cell gathering” based approaches. The methods directly using Voronoi polygons of vertices to form the boundaries are simple and straightforward. But they may not be very suitable in the sectorization problem for two reasons. Firstly, the shape of sectors is not strictly convex so that Voronoi polygons may not fit the outline of concave vertex cluster very well for its characteristic of convexity. Secondly, the MDFB constraint cannot be ensured by Voronoi polygons.We propose a concave hull-based method to depict the sector boundaries in this section. We firstly construct the concave hull for each vertex cluster (representing one sector) in Procedure2 with necessary explanations below.Procedure 2:Procedures to construct Concave(Bs).
(1) Initialize, input G=⋃s∈SGs, V=⋃s∈SVs, MP=⋃s∈SMPs, and OP=⋃sOPs (2) Fors∈S (3)   Set Bs=∅ (4)   Ifs is adjacent to the outer-boundary of the target airspace (5)    Obtain BMPs (See explain in Section 4.3.1), and set NBMs=∅ (6)    For each point pt∈BMPs (7)    Find point pto, that pto∈OP, and pto is nearest to pt (8)    NBMs=NBMs∪pto (9)    End (10)  Obtain OPs (See explain in Section 4.3.1) (11)  Bs=OPs∪MPs (12) End (13) Else (14)  Bs=MPs (15) End (16) Construct Concave(Bs) using α-shape algorithm (See explain in Section 4.3.1) (17) End


### 4.3.1. Concave Hull Generation

Input Data. Input information of the concave hull generation method is in three categories.The first is the cut ofG: (19)G=⋃s=1kGs,V=⋃s=1kVs.The second is the set of midpoints of edges connecting two neighboring sectors,MP=mpi,j∣eij=1,vi∈Vsm,vj∈Vsn,sm≠sn. For any sector s, define the set of midpoints of edges connecting it and other sectors as MPs=mpi,j∣eij=1,vi∈Vs,vj∉Vs. Hence, we get that MP=⋃s∈SMPs and MPs1∩MPs2≠∅ if s1 and s2 are neighboring sectors.The third is a set of points located on the outerboundary of the target airspace. We suppose points on the outerboundary of the target airspace are distributed clockwise for the purpose of unification and define the point set asOP=opi∣1≤i≤n, OP=⋃sOPs, s∈S, and OPs=∅ if s is not next to the outerboundary.Obtain B M P s . For any sector s (subgraph Gs) adjacent to the outerboundary, we find out one or more couples of points in MPs, and denote the subset of MPs as BMPs. For any point mpi,j∈BMPs with corresponding edge eij, the following conditions must be satisfied:(i)
Both endpoints ofeij, denoted as vi and vj, are next to the outerboundary.(ii)
The Voronoi polygons ofvi and vj have one common side and this side intersects with the outerboundary.Note that, as the shape of any sector is not strictly convex, there may be more than one couple of points inBMPs. In Figure 5(a), we assume that V2 and V3 belong to sector s and V1 and V4 belong to the other two different sectors, respectively. In this situation, BMPs is made up of MP1 and MP2. In Figure 5(b), we assume that sector s consists of V2, V3, and V4, which surround V1. The shape of s is concave, but s still conforms to the convexity constraint because no air routes connect V3V1 or V4V1. In this situation, MP1, MP2, MP3, and MP4 are all elements of BMPs. BMPs with more than two couples of points is also possible which is not illustrated here. Note that BMPs1∩BMPs2≠∅, if s1 and s2 are neighboring sectors.Figure 5
Midpoints adjacent to the outerboundary.
(a) (b)Obtain O P s . For any sector s adjacent to the outerboundary, sort NBMs clockwise (as points in NBMs all belong to OP, sorting clockwise means sorting with increasing index in OP) as NBMs=op1s,op2s,…,op2k-1s,op2ks, where k is the number of couples nBMPs. Then, we define (20)OPs=⋃i=1kop2i-1s,…,op2is.We are aware that two neighboring points in NBMs may not be consecutive in OP. We define(21)op2i-1s,…,op2is=opind1,opind1+1,…,opind2-1,opind2,where ind1 is the index of op2i-1s in OP and ind2 is the index of op2is in OP.Construct Concave Hull of B s. We use the concave hull of Bs, denoted as ConcaveBs, to fit the outline of the vertex cluster of sector s. α-shape algorithm [38] is applied here to generate the concave hull. The basic idea of α-shape algorithm is to iteratively search α-extreme points and α-neighbors and to connect the found α-neighbors sequentially to generate α-hull. As large amount of extra definitions and notations must be introduced to formulate α-shape algorithm clearly, we do not present the details here.


### 4.3.2. Concave Hull Amendment

C o n c a v e B s is a promising candidate to be the boundary of sector s because all vertices of Gs are inside it and it fits the shape of s very well. However, two more considerations make it necessary to recheck ConcaveBs. The first is the compactness of the boundary; that is, fewer segments are preferred on the boundary. The second is about the MDFB constraint; that is, the distance between some segments of ConcaveBs and some vertices of G may be less than the MDFB. Hence, we propose the Concave Hull Amendment Method (CHAM) to deal with these two considerations (Procedure 3).Procedure 3:Procedures of the CHAM.
(1) Fors∈S (2)   Get path(s) on Concave(Bs), set as P1s,P2s,…,Pkss (See Definition 4 in Section 4.3.2) (3)   For each path Pis, 1≤i≤ks (4)    Define vertices along Pis as V1p,V2p,…,VNip and corresponding BI1p,BI2p,…,BINip (note thatBIjp=Vjp when corresponding cross-sector edge has not been designed) (5)    Set j=1, gap=1. (6)    Whilej<Ni (7)     If cross-sector edges referred by BIjp and BIj+1p are both designed (8)      If line segment BIjpBIj+1p is infeasible (See Definition 5 in Section 4.3.2) (9)        Execute the OABD to BIjpBIj+1p (See Explain in Section 4.3.2) (10)    End (11)    j=j+1, go to line (6) (12)   End (13)   While line segment BIjpBIj+gapp is feasible (See Definition 5 in Section 4.3.2) (14)    gap=gap+1 (15)    If cross-sector edges referred by BIj+gap-1p and BIj+gapp are both designed (16)     break (17)    End (18)   End (19)   Ifgap=1 (20)    Execute the OABD to BIjpBIj+1p (See Explain in Section 4.3.2). (21)    Set the cross-sector edge referred by BIj+1p as designed, remain BIj+1p=Vj+1p (22)    j=j+1, go to line (6) (23)   End (24)   Else (25)    gap=gap-1 (26)    Draw line segment BIjpBIj+gapp, (27)    Set cross-sector edges referred by BIj+1p,BIj+2p,…,BIj+gapp as designed (28)    Set BIj+1p,BIj+2p,…,BIj+gapp as the intersections of related cross-sector edges and line segmentBIjpBIj+gapp (29)    j=j+gap, reset gap=1, go to line (6) (30)   End (31)  End (32) End (33) EndThere are two ideas in theCHAM. The first is to replace the polygonal line connecting several vertices on the concave hull with one line segment from the first vertex to the last one. The second idea is to avoid obstacle circles centering at the vertices of G with radius of that vertex (refer to Section 3.1 for the definition of the radius of a vertex). We call modified concave hulls after theCHAM boundary polygons and present two definitions about the concave hull and the boundary polygon first.Definition 4 (path(s) on the concave hull).
One scans the vertices ofConcaveBs clockwise. If s is adjacent to the outerboundary of the target airspace, whenever we encounter a sequence of vertices V1p-V2p-⋯-Vr-1p,Vrp, where V1p, Vrp∈OP or others belonging to MP, we define such a sequence as a path on the concave hull. If s is not adjacent to the outerboundary, that is, it is surrounded by other sectors, we simply define a path on ConcaveBs as a clockwise closed walk of all vertices of ConcaveBs.Definition 5.
A line segment connecting any two vertices on the boundary polygon of one sector is defined as feasible if it does not intersect with any inner-edges of any sector and it does not go through any obstacle circles centering at the vertices ofG.Next, we present another definition to avoid redepiction of common boundary between two neighboring sectors.Definition 6.
In theCHAM, once a segment on the boundary polygon is depicted, set all cross-sector edges that intersect with this segment as designed. Define the intersection of this segment and a cross-sector edge, eij, as the boundary crossing point of eij, denoted as BIij. When one cross-sector edge has not been designed, its default boundary crossing point is set as its midpoint, that is, the intersection of this edge and the concave hull. To ensure unification, when one endpoint of a depicted segment lies on the outerboundary of the target airspace, we also define a virtual cross-sector edge corresponding to this endpoint as designed and set the boundary crossing point as the endpoint itself.We also give two properties about the concave hull.Property 1.
For any sectors, line segments on ConcaveBs do not intersect with any inner-edges of any sector. This property is obvious from the definition and the generating process of α-shape [38].Property 2.
The distance between any vertex on the concave hull, representing the midpoint of a cross-sector edge, and any vertex ofG is at least MDFB.Proof.
For any vertex on the concave hull, if the corresponding cross-sector edge represents a real air-route segment, the edge is at least MDFB away from any vertex ofG according to Rule 3 in Section 3.2 except its two endpoints. In addition, as all edges of G are at least twice as long as MDFB, their midpoints are at least MDFB away from the two endpoints. For one vertex on the concave hull that is the midpoint of a sector-cross edge only standing for Voronoi neighboring relationship, it locates on the Voronoi polygons of both endpoints of this edge, which means that the nearest vertices of G to this midpoint are just the two endpoints of the edge. Considering the possible distance between this midpoint and the two endpoints is at least MDFB, if any other vertex of G is less than MDFB from this midpoint, a contradictory occurs.These two properties show two pieces of important information when we construct the boundary polygons:(1)
For any two consecutive vertices on the concave hull of a sector, if we can find a line segment or a polygonal line that does not go through any obstacle circles, this line segment or polygonal line can be set as boundary segment(s) of that sector.(2)
The boundary polygons after theCHAM can cross any vertices on the concave hulls.We show the procedures of theCHAM in Table 1 with necessary explanations below.Table 1
Results of six sectorization methods.
Time period Sectorization method Obj1 Obj2 Obj3 Con1 Con2 Con3 Nos Density gap RT (s) 12:00–14:00 Current in practice 0.8145 674 238 0 1 0 16 4218 — Only NSCA 1.0518 471 133 0 6 1 16 6027 — RI 0.3158 404 84 1 0 0 9 2521 934 CI 0.3124 360 59 0 0 0 9 2232 940 NI 0.3404 375 58 0 0 0 8 2356 942 MI 0.4329 386 68 0 0 0 8 2512 964  16:00–18:00 Current in practice 1.1466 707 234 0 4 2 16 6705 — Only NSCA 1.6230 428 86 0 0 3 16 6732 — RI 0.3148 525 106 1 0 0 10 2409 926 CI 0.4508 608 130 0 0 0 10 3007 937 NI 0.4259 499 104 0 0 0 10 2521 924 MI 0.4502 511 91 0 0 0 10 2619 940  20:00–22:00 Current in practice 1.3666 710 247 0 3 2 16 7198 — Only NSCA 2.0712 403 72 0 1 5 16 9150 — RI 0.4514 684 206 2 0 0 12 2375 936 CI 0.3487 719 218 0 0 0 12 1902 934 NI 0.4142 592 134 0 0 0 12 2258 930 MI 0.3993 598 136 0 0 0 12 2017 915  2:00–24:00 Current in practice 1.3850 427 166 0 1 2 16 7604 — Only NSCA 1.6230 428 86 1 2 3 16 9058 — RI 0.4318 468 205 2 0 0 12 2635 944 CI 0.4100 335 79 0 0 0 12 2239 937 NI 0.3678 350 105 0 0 0 12 2159 924 MI 0.2982 351 118 0 0 0 12 2073 930Obstacle Avoidance Boundary Design (OABD). We illustrate theOABD with Figure 6. Firstly, we assume two consecutive vertices on the concave hull as Vs and Ve. VsVe is assumed to intrude one or more obstacle circles (Figure 6(a)). TheOABD examines the first obstacle (such as Ob1 in Figure 6) intruded by VsVe. Then, we draw two tangents from Vs and Ve to the obstacle circle, respectively. The intersection of these two tangents is denoted as V0. Thus, the boundary between Vs and Ve is now VsV0Ve (Figure 6(b)). Next, we examine if VsV0Ve still intrudes any other obstacle (such as Ob2 in Figure 6). If it does, we truncate V0Ve at the boundary of the second obstacle (V1 in Figure 6(c)). Then, we draw two tangents of Ob2 originating from V1 and Ve, respectively, and get their intersection (V2 in Figure 6(c)). The boundary between Vs and Ve is now changed to VsV0V1V2Ve. TheOABD continues until the curve VsV0Ve (1 iteration) or VsV0V1⋯Vk+1Ve (k iterations, k>1) is free of all obstacles (Figure 6(c)). If too many iterations are needed in theOABD, the generated sector boundaries may be jagged, which is prohibited. Throughout our simulation, such situations do not appear. Relevant discussions will be in Section 6.5.Figure 6
Procedures of theOABD.
(a) (b) (c)In theOABD, there are always two intersections when we draw tangents, such as I1 and I2, in Figure 7. As I1 and I2 always locate symmetrically on opposite sides of the obstacle circle, one of the generated polygonal lines definitely intersects with inner-edges of sectors while the other does not. The one without intersections with inner-edges is our choice (such as VsI2Ve in Figure 7).Figure 7
Selection when drawing tangents.


### 4.3.3. Vacuum Space Handle

After theCHAM, there still remains several vacuum spaces at the intersecting areas of three or more sectors, such as quadrilateralABCD in Figure 8(a). These vacuum spaces contain no vertices or edges and are not assigned to any sector. In this section, we propose a heuristic method to assign vacuum spaces. Firstly, we represent a vacuum space with a polygon, calling the vacuum polygon. The vacuum polygons can be sorted into four categories according to types of their sides:(I)
Triangles with no side from theOABD.(II)
Polygons with more than three sides and no side that is from theOABD.(III)
Quadrangles with two sides from theOABD.(IV)
Other types.Figure 8
Vacuum polygon handle, Category (II).
(a) (b) (c)For vacuum polygons in Category (I), we propose two choices, one is directly erasing the opposite side of the biggest inner angle (Figure9(a)), and the other is to find the “most vertical” midline and erase those two sides on which the midpoint does not locate (Figure 9(b)).Figure 9
Vacuum polygon handle, Categories (I) and (III).
(a) (b) (c)For vacuum polygons in Category (II), we first give the following property.Property 3.
The vacuum polygon in Category (II) is convex.Proof.
Suppose there is a concave vacuum polygon in Category (II) (such as the one in Figure8(b)). We consider the vertex causing the concavity (point C in Figure 8(b)). As this vertex is on a cross-sector edge, one of the endpoints of this cross-sector edge locates inside the vacuum polygon because of the concavity, obviously causing a contradictory against the definition of the vacuum space. In another situation, it locates outside the polygon while the cross-sector edge must intersect with one extra side of the polygon. However, neither of the two endpoints of this extra side is the vertex that causes the concavity; that is, a cross-sector edge between two sectors intersects with the boundary of an irrelevant sector. It is also impossible. Hence, we show that the vacuum polygons in Category (II) must be convex.With Property3, we can calculate the geometric center of a polygon in Category (II), which must locate inside the polygon because of its convexity (point O in Figure 8(b)). Then, we connect the center to every vertices of the polygon and erase original sides of the polygon (Figure 8(c)).For polygons in Category (III), which stand for an intersecting area of three sectors, we find the vertex produced by the OABD (pointI in Figure 9(c)), draw the diagonal from it, and erase those two opposite sides of the vertex produced by the OABD.For polygons in Category (IV), situation of their sides may be very complex. We directly assign the entire vacuum space to one sector, whose boundary segments are sides of the polygon afterOABD.


## 4.3.1. Concave Hull Generation

Input Data. Input information of the concave hull generation method is in three categories.The first is the cut ofG: (19)G=⋃s=1kGs,V=⋃s=1kVs.The second is the set of midpoints of edges connecting two neighboring sectors,MP=mpi,j∣eij=1,vi∈Vsm,vj∈Vsn,sm≠sn. For any sector s, define the set of midpoints of edges connecting it and other sectors as MPs=mpi,j∣eij=1,vi∈Vs,vj∉Vs. Hence, we get that MP=⋃s∈SMPs and MPs1∩MPs2≠∅ if s1 and s2 are neighboring sectors.The third is a set of points located on the outerboundary of the target airspace. We suppose points on the outerboundary of the target airspace are distributed clockwise for the purpose of unification and define the point set asOP=opi∣1≤i≤n, OP=⋃sOPs, s∈S, and OPs=∅ if s is not next to the outerboundary.Obtain B M P s . For any sector s (subgraph Gs) adjacent to the outerboundary, we find out one or more couples of points in MPs, and denote the subset of MPs as BMPs. For any point mpi,j∈BMPs with corresponding edge eij, the following conditions must be satisfied:(i)
Both endpoints ofeij, denoted as vi and vj, are next to the outerboundary.(ii)
The Voronoi polygons ofvi and vj have one common side and this side intersects with the outerboundary.Note that, as the shape of any sector is not strictly convex, there may be more than one couple of points inBMPs. In Figure 5(a), we assume that V2 and V3 belong to sector s and V1 and V4 belong to the other two different sectors, respectively. In this situation, BMPs is made up of MP1 and MP2. In Figure 5(b), we assume that sector s consists of V2, V3, and V4, which surround V1. The shape of s is concave, but s still conforms to the convexity constraint because no air routes connect V3V1 or V4V1. In this situation, MP1, MP2, MP3, and MP4 are all elements of BMPs. BMPs with more than two couples of points is also possible which is not illustrated here. Note that BMPs1∩BMPs2≠∅, if s1 and s2 are neighboring sectors.Figure 5
Midpoints adjacent to the outerboundary.
(a) (b)Obtain O P s . For any sector s adjacent to the outerboundary, sort NBMs clockwise (as points in NBMs all belong to OP, sorting clockwise means sorting with increasing index in OP) as NBMs=op1s,op2s,…,op2k-1s,op2ks, where k is the number of couples nBMPs. Then, we define (20)OPs=⋃i=1kop2i-1s,…,op2is.We are aware that two neighboring points in NBMs may not be consecutive in OP. We define(21)op2i-1s,…,op2is=opind1,opind1+1,…,opind2-1,opind2,where ind1 is the index of op2i-1s in OP and ind2 is the index of op2is in OP.Construct Concave Hull of B s. We use the concave hull of Bs, denoted as ConcaveBs, to fit the outline of the vertex cluster of sector s. α-shape algorithm [38] is applied here to generate the concave hull. The basic idea of α-shape algorithm is to iteratively search α-extreme points and α-neighbors and to connect the found α-neighbors sequentially to generate α-hull. As large amount of extra definitions and notations must be introduced to formulate α-shape algorithm clearly, we do not present the details here.


## 4.3.2. Concave Hull Amendment

C o n c a v e B s is a promising candidate to be the boundary of sector s because all vertices of Gs are inside it and it fits the shape of s very well. However, two more considerations make it necessary to recheck ConcaveBs. The first is the compactness of the boundary; that is, fewer segments are preferred on the boundary. The second is about the MDFB constraint; that is, the distance between some segments of ConcaveBs and some vertices of G may be less than the MDFB. Hence, we propose the Concave Hull Amendment Method (CHAM) to deal with these two considerations (Procedure 3).Procedure 3:Procedures of the CHAM.
(1) Fors∈S (2)   Get path(s) on Concave(Bs), set as P1s,P2s,…,Pkss (See Definition 4 in Section 4.3.2) (3)   For each path Pis, 1≤i≤ks (4)    Define vertices along Pis as V1p,V2p,…,VNip and corresponding BI1p,BI2p,…,BINip (note thatBIjp=Vjp when corresponding cross-sector edge has not been designed) (5)    Set j=1, gap=1. (6)    Whilej<Ni (7)     If cross-sector edges referred by BIjp and BIj+1p are both designed (8)      If line segment BIjpBIj+1p is infeasible (See Definition 5 in Section 4.3.2) (9)        Execute the OABD to BIjpBIj+1p (See Explain in Section 4.3.2) (10)    End (11)    j=j+1, go to line (6) (12)   End (13)   While line segment BIjpBIj+gapp is feasible (See Definition 5 in Section 4.3.2) (14)    gap=gap+1 (15)    If cross-sector edges referred by BIj+gap-1p and BIj+gapp are both designed (16)     break (17)    End (18)   End (19)   Ifgap=1 (20)    Execute the OABD to BIjpBIj+1p (See Explain in Section 4.3.2). (21)    Set the cross-sector edge referred by BIj+1p as designed, remain BIj+1p=Vj+1p (22)    j=j+1, go to line (6) (23)   End (24)   Else (25)    gap=gap-1 (26)    Draw line segment BIjpBIj+gapp, (27)    Set cross-sector edges referred by BIj+1p,BIj+2p,…,BIj+gapp as designed (28)    Set BIj+1p,BIj+2p,…,BIj+gapp as the intersections of related cross-sector edges and line segmentBIjpBIj+gapp (29)    j=j+gap, reset gap=1, go to line (6) (30)   End (31)  End (32) End (33) EndThere are two ideas in theCHAM. The first is to replace the polygonal line connecting several vertices on the concave hull with one line segment from the first vertex to the last one. The second idea is to avoid obstacle circles centering at the vertices of G with radius of that vertex (refer to Section 3.1 for the definition of the radius of a vertex). We call modified concave hulls after theCHAM boundary polygons and present two definitions about the concave hull and the boundary polygon first.Definition 4 (path(s) on the concave hull).
One scans the vertices ofConcaveBs clockwise. If s is adjacent to the outerboundary of the target airspace, whenever we encounter a sequence of vertices V1p-V2p-⋯-Vr-1p,Vrp, where V1p, Vrp∈OP or others belonging to MP, we define such a sequence as a path on the concave hull. If s is not adjacent to the outerboundary, that is, it is surrounded by other sectors, we simply define a path on ConcaveBs as a clockwise closed walk of all vertices of ConcaveBs.Definition 5.
A line segment connecting any two vertices on the boundary polygon of one sector is defined as feasible if it does not intersect with any inner-edges of any sector and it does not go through any obstacle circles centering at the vertices ofG.Next, we present another definition to avoid redepiction of common boundary between two neighboring sectors.Definition 6.
In theCHAM, once a segment on the boundary polygon is depicted, set all cross-sector edges that intersect with this segment as designed. Define the intersection of this segment and a cross-sector edge, eij, as the boundary crossing point of eij, denoted as BIij. When one cross-sector edge has not been designed, its default boundary crossing point is set as its midpoint, that is, the intersection of this edge and the concave hull. To ensure unification, when one endpoint of a depicted segment lies on the outerboundary of the target airspace, we also define a virtual cross-sector edge corresponding to this endpoint as designed and set the boundary crossing point as the endpoint itself.We also give two properties about the concave hull.Property 1.
For any sectors, line segments on ConcaveBs do not intersect with any inner-edges of any sector. This property is obvious from the definition and the generating process of α-shape [38].Property 2.
The distance between any vertex on the concave hull, representing the midpoint of a cross-sector edge, and any vertex ofG is at least MDFB.Proof.
For any vertex on the concave hull, if the corresponding cross-sector edge represents a real air-route segment, the edge is at least MDFB away from any vertex ofG according to Rule 3 in Section 3.2 except its two endpoints. In addition, as all edges of G are at least twice as long as MDFB, their midpoints are at least MDFB away from the two endpoints. For one vertex on the concave hull that is the midpoint of a sector-cross edge only standing for Voronoi neighboring relationship, it locates on the Voronoi polygons of both endpoints of this edge, which means that the nearest vertices of G to this midpoint are just the two endpoints of the edge. Considering the possible distance between this midpoint and the two endpoints is at least MDFB, if any other vertex of G is less than MDFB from this midpoint, a contradictory occurs.These two properties show two pieces of important information when we construct the boundary polygons:(1)
For any two consecutive vertices on the concave hull of a sector, if we can find a line segment or a polygonal line that does not go through any obstacle circles, this line segment or polygonal line can be set as boundary segment(s) of that sector.(2)
The boundary polygons after theCHAM can cross any vertices on the concave hulls.We show the procedures of theCHAM in Table 1 with necessary explanations below.Table 1
Results of six sectorization methods.
Time period Sectorization method Obj1 Obj2 Obj3 Con1 Con2 Con3 Nos Density gap RT (s) 12:00–14:00 Current in practice 0.8145 674 238 0 1 0 16 4218 — Only NSCA 1.0518 471 133 0 6 1 16 6027 — RI 0.3158 404 84 1 0 0 9 2521 934 CI 0.3124 360 59 0 0 0 9 2232 940 NI 0.3404 375 58 0 0 0 8 2356 942 MI 0.4329 386 68 0 0 0 8 2512 964  16:00–18:00 Current in practice 1.1466 707 234 0 4 2 16 6705 — Only NSCA 1.6230 428 86 0 0 3 16 6732 — RI 0.3148 525 106 1 0 0 10 2409 926 CI 0.4508 608 130 0 0 0 10 3007 937 NI 0.4259 499 104 0 0 0 10 2521 924 MI 0.4502 511 91 0 0 0 10 2619 940  20:00–22:00 Current in practice 1.3666 710 247 0 3 2 16 7198 — Only NSCA 2.0712 403 72 0 1 5 16 9150 — RI 0.4514 684 206 2 0 0 12 2375 936 CI 0.3487 719 218 0 0 0 12 1902 934 NI 0.4142 592 134 0 0 0 12 2258 930 MI 0.3993 598 136 0 0 0 12 2017 915  2:00–24:00 Current in practice 1.3850 427 166 0 1 2 16 7604 — Only NSCA 1.6230 428 86 1 2 3 16 9058 — RI 0.4318 468 205 2 0 0 12 2635 944 CI 0.4100 335 79 0 0 0 12 2239 937 NI 0.3678 350 105 0 0 0 12 2159 924 MI 0.2982 351 118 0 0 0 12 2073 930Obstacle Avoidance Boundary Design (OABD). We illustrate theOABD with Figure 6. Firstly, we assume two consecutive vertices on the concave hull as Vs and Ve. VsVe is assumed to intrude one or more obstacle circles (Figure 6(a)). TheOABD examines the first obstacle (such as Ob1 in Figure 6) intruded by VsVe. Then, we draw two tangents from Vs and Ve to the obstacle circle, respectively. The intersection of these two tangents is denoted as V0. Thus, the boundary between Vs and Ve is now VsV0Ve (Figure 6(b)). Next, we examine if VsV0Ve still intrudes any other obstacle (such as Ob2 in Figure 6). If it does, we truncate V0Ve at the boundary of the second obstacle (V1 in Figure 6(c)). Then, we draw two tangents of Ob2 originating from V1 and Ve, respectively, and get their intersection (V2 in Figure 6(c)). The boundary between Vs and Ve is now changed to VsV0V1V2Ve. TheOABD continues until the curve VsV0Ve (1 iteration) or VsV0V1⋯Vk+1Ve (k iterations, k>1) is free of all obstacles (Figure 6(c)). If too many iterations are needed in theOABD, the generated sector boundaries may be jagged, which is prohibited. Throughout our simulation, such situations do not appear. Relevant discussions will be in Section 6.5.Figure 6
Procedures of theOABD.
(a) (b) (c)In theOABD, there are always two intersections when we draw tangents, such as I1 and I2, in Figure 7. As I1 and I2 always locate symmetrically on opposite sides of the obstacle circle, one of the generated polygonal lines definitely intersects with inner-edges of sectors while the other does not. The one without intersections with inner-edges is our choice (such as VsI2Ve in Figure 7).Figure 7
Selection when drawing tangents.


## 4.3.3. Vacuum Space Handle

After theCHAM, there still remains several vacuum spaces at the intersecting areas of three or more sectors, such as quadrilateralABCD in Figure 8(a). These vacuum spaces contain no vertices or edges and are not assigned to any sector. In this section, we propose a heuristic method to assign vacuum spaces. Firstly, we represent a vacuum space with a polygon, calling the vacuum polygon. The vacuum polygons can be sorted into four categories according to types of their sides:(I)
Triangles with no side from theOABD.(II)
Polygons with more than three sides and no side that is from theOABD.(III)
Quadrangles with two sides from theOABD.(IV)
Other types.Figure 8
Vacuum polygon handle, Category (II).
(a) (b) (c)For vacuum polygons in Category (I), we propose two choices, one is directly erasing the opposite side of the biggest inner angle (Figure9(a)), and the other is to find the “most vertical” midline and erase those two sides on which the midpoint does not locate (Figure 9(b)).Figure 9
Vacuum polygon handle, Categories (I) and (III).
(a) (b) (c)For vacuum polygons in Category (II), we first give the following property.Property 3.
The vacuum polygon in Category (II) is convex.Proof.
Suppose there is a concave vacuum polygon in Category (II) (such as the one in Figure8(b)). We consider the vertex causing the concavity (point C in Figure 8(b)). As this vertex is on a cross-sector edge, one of the endpoints of this cross-sector edge locates inside the vacuum polygon because of the concavity, obviously causing a contradictory against the definition of the vacuum space. In another situation, it locates outside the polygon while the cross-sector edge must intersect with one extra side of the polygon. However, neither of the two endpoints of this extra side is the vertex that causes the concavity; that is, a cross-sector edge between two sectors intersects with the boundary of an irrelevant sector. It is also impossible. Hence, we show that the vacuum polygons in Category (II) must be convex.With Property3, we can calculate the geometric center of a polygon in Category (II), which must locate inside the polygon because of its convexity (point O in Figure 8(b)). Then, we connect the center to every vertices of the polygon and erase original sides of the polygon (Figure 8(c)).For polygons in Category (III), which stand for an intersecting area of three sectors, we find the vertex produced by the OABD (pointI in Figure 9(c)), draw the diagonal from it, and erase those two opposite sides of the vertex produced by the OABD.For polygons in Category (IV), situation of their sides may be very complex. We directly assign the entire vacuum space to one sector, whose boundary segments are sides of the polygon afterOABD.


## 5. Configuration Transition

The methods presented from Sections2–4 are for the sectorization problem in single time period. This section focuses on solving the CTP. First of all, we present the definition of the similarity of two different configurations followed by its calculation method.


### 5.1. Similarity of Two Configurations

We still use graph-based notations to denote configurations. We still model the whole target airspace asG, and define C1 and C2 as two different configurations:(22)C1:G=⋃s=1k1Gs1,V=⋃s=1k1Vs1,C2:G=⋃s=1k2Gs2,V=⋃s=1k2Vs2.We construct a bisection diagram (Figure10(a)), where nodes on each side stand for the sectors in each configuration. The weight of an edge connecting two nodes on opposite sides is the similarity of two sectors, which is defined as the number of same vertices of these two sectors:(23)SimGs11,Gs22=Vs11∩Vs22s1∈1,2,…,k1,s2∈1,2,…,k2.Figure 10
Configuration similarity and transition.
(a) (b) (c)To calculate the similarity ofC1 and C2, we need to get the maximum match of the bisection diagram, which ensures a monoprojection for every node on the side with less nodes while the aggregate weights of selected edges are maximized. Without loss of generality, we suppose k1≤k2 (it is similar when k2≤k1), the similarity of C1 and C2 is defined as:(24)SimC1,C2=max⁡∑s1=1k1∑s2=1k2xs2s1ωs2s1,∑s2=1k2xs2s1=1,∀s1∈1,2,…,k1,∑s1=1k1xs2s1≤1,∀s2∈1,2,…,k2,xs2s1∈0,1,where ωs2s1=SimGs11,Gs22.This is a typical 0-1 programming problem, which can be solved with network flow algorithms. We modify the bisection diagram (Figure10(b)), adding one virtual source node to the side with less nodes and connect the source node with each node on the same side. Besides, we also add a virtual sink node to the side with more nodes, connecting the sink node with each node on the same side. The weights on all newly added edges are 0 so that they have no impact on the result. The weights on original edges are also changed. For an edge connecting Gs11 and Gs22, its weight is changed to M-SimGs11,Gs22 from SimGs11,Gs22, where M is a constant. So the objective function is changed to the form of k1M-f(x). Because k1 is known, original maximization problem in (24) is now transformed to a minimization problem. The supply of the source node and the demand of the sink node are both equal to k1 and the supply/demand at other nodes are all 0. Finally, the capacities on all edges are set as 1. Based on these modifications, we build a typical network that applies to any generic minimum cost flow algorithms [29]. We do not list the procedures of the minimum cost flow algorithms in detail. Readers can easily find them in [29] or any other textbooks on network flow or Discrete Optimization.


### 5.2. Nondominated Configuration Links

If we want to make a configuration plan for next several time periods, we should consider the quality of selected configurations in terms of the values inObj1 toObj3 and the similarities between configurations of consecutive time periods. Hence, finding an optimal configuration plan is another multiobjective problem.We build a forward network to represent the CTP (Figure10(c)). Each column of nodes represents one of several nondominated configurations in one period. A node in the network is denoted as Cti, 1≤t≤T, where T is the number of time periods and 1≤i≤kt, where kt is the number of nondominated configurations in period t. Each edge in the network stands for a possible transition from the previous period to the next period. Edges standing for transitions with too low similarity are truncated.We define two kinds of parameters in the network. The first one is the weight of edgeet,t+1ij, denoted as wt,t+1ij, set as the inverse of the similarity of configuration Cti and Ct+1j. The second one is the values ofObj1 toObj3 of node Cti, denoted as fti=ft,1i,ft,2i,ft,3i.Definition 7.
One defines a route in the forward network from one node to any node in the last period as a configuration link. A link without its starting node is called a path.For nodeCti, one denotes the cost along a path pi,tp from it to any node in the last period as gpi,tp=g1pi,tp,g2pi,tp,g3pi,tp,g4pi,tp, where g1pi,tp is the sum of similarities along the path and the other three items correspond to the sum of values of three objective functions of nodes along the path. Hence, one can define the dominating relationship between two paths.Definition 8.
Assuming two paths form the same time periodt to the last period as pi,tp and pj,tq, pi,tp≺pj,tq, if gpi,tp≤gpj,tq and ∃m,1≤m≤4, gmpi,tp<gmpj,tq.From Definition8, we can define a nondominated set of paths from node Cti, denoted as NDti. Then, we assume another path from one node in previous period, Ct-1j, which goes across Cti as pj,t-1r, and the part after Cti along pj,t-1r is pi,tp. Then, we can get the following relationships:(25a)pj,t-1r=pi,tp∪et-1,tji,Cti,(25b)g1pj,t-1r=g1pi,tp+wt-1,tji,(25c)gmpj,t-1r=gmpi,tp+ft,m-1i,m=2,3,4.Hence, we obtain the following proposition.Proposition 9.
Assuming a nondominated set of paths originating fromCt-1j and paths in this set go across Cti, the part after Cti along such path belongs to NDti.Proof.
Assume two paths originating fromCt-1j and going across Cti as pj,t-1p and pj,t-1q. Also, assume their parts after Cti as pi,tr and pi,ts, respectively. Considering the calculations in (25b) and (25c), if pi,tr≺pi,ts, then pj,t-1p≺pj,t-1q.Proposition9 shows that the subpaths not being in NDti can be discarded when we calculate NDt-1j. We consider all nodes, Cti, 1≤i≤kt, in period t, which can be transitioned from Ct-1j. We denote the nondominated paths originating from Ct-1j and going across Cti as NDt-1,tji, and we get (26)ND⌣t-1j=⋃iNDt-1,tji.We elite out nondominated paths in ND⌣t-1j and get the nondominated set of paths originating from Ct-1j, denoted as NDt-1j.We have presented the procedure to obtainNDt-1j so far when NDti is known. This is similar to one of the steps in typical dynamic programming except multiple objectives. In fact, theprinciple of optimality in monoobjective DP also works in multiobjective cases [39] and Proposition 9 just acts the same as theprinciple of optimality.Hence, we apply a backward multiobjective DP algorithm starting at stageT-1, getting the nondominated path set for each node in period T-1. In stage n, 1≤n≤T-1, we deal with edges from period t to t+1 and the nodes in period t+1 to generate sets of nondominated paths from nodes in period t.After the calculation in stage 1, we generate nondominated path sets from each node in period 1, denoted asND11,ND12,…,ND1k1. However, values of objective functions of configurations in period 1 are not involved yet. Next, we integrate nodes in period 1 into the paths to get final nondominated links.For each configuration in period 1,C1i, 1≤i≤k1 and corresponding ND1i, suppose a path in ND1i as pi,1p and the link li,1p that integrates C1i with pi,1p; then(27a)li,1p=pi,1p∪C1i,(27b)gmli,1p=gmpi,1p+f1,m-1i,m=2,3,4.We take operations like (27a) and (27b) for all nondominated paths in ND1i, 1≤i≤k1. Then, we classify the generated links into a set FND⌣. Finally, we select nondominated individuals from FND⌣ to obtain final nondominated configuration link set FND.


## 5.1. Similarity of Two Configurations

We still use graph-based notations to denote configurations. We still model the whole target airspace asG, and define C1 and C2 as two different configurations:(22)C1:G=⋃s=1k1Gs1,V=⋃s=1k1Vs1,C2:G=⋃s=1k2Gs2,V=⋃s=1k2Vs2.We construct a bisection diagram (Figure10(a)), where nodes on each side stand for the sectors in each configuration. The weight of an edge connecting two nodes on opposite sides is the similarity of two sectors, which is defined as the number of same vertices of these two sectors:(23)SimGs11,Gs22=Vs11∩Vs22s1∈1,2,…,k1,s2∈1,2,…,k2.Figure 10
Configuration similarity and transition.
(a) (b) (c)To calculate the similarity ofC1 and C2, we need to get the maximum match of the bisection diagram, which ensures a monoprojection for every node on the side with less nodes while the aggregate weights of selected edges are maximized. Without loss of generality, we suppose k1≤k2 (it is similar when k2≤k1), the similarity of C1 and C2 is defined as:(24)SimC1,C2=max⁡∑s1=1k1∑s2=1k2xs2s1ωs2s1,∑s2=1k2xs2s1=1,∀s1∈1,2,…,k1,∑s1=1k1xs2s1≤1,∀s2∈1,2,…,k2,xs2s1∈0,1,where ωs2s1=SimGs11,Gs22.This is a typical 0-1 programming problem, which can be solved with network flow algorithms. We modify the bisection diagram (Figure10(b)), adding one virtual source node to the side with less nodes and connect the source node with each node on the same side. Besides, we also add a virtual sink node to the side with more nodes, connecting the sink node with each node on the same side. The weights on all newly added edges are 0 so that they have no impact on the result. The weights on original edges are also changed. For an edge connecting Gs11 and Gs22, its weight is changed to M-SimGs11,Gs22 from SimGs11,Gs22, where M is a constant. So the objective function is changed to the form of k1M-f(x). Because k1 is known, original maximization problem in (24) is now transformed to a minimization problem. The supply of the source node and the demand of the sink node are both equal to k1 and the supply/demand at other nodes are all 0. Finally, the capacities on all edges are set as 1. Based on these modifications, we build a typical network that applies to any generic minimum cost flow algorithms [29]. We do not list the procedures of the minimum cost flow algorithms in detail. Readers can easily find them in [29] or any other textbooks on network flow or Discrete Optimization.


## 5.2. Nondominated Configuration Links

If we want to make a configuration plan for next several time periods, we should consider the quality of selected configurations in terms of the values inObj1 toObj3 and the similarities between configurations of consecutive time periods. Hence, finding an optimal configuration plan is another multiobjective problem.We build a forward network to represent the CTP (Figure10(c)). Each column of nodes represents one of several nondominated configurations in one period. A node in the network is denoted as Cti, 1≤t≤T, where T is the number of time periods and 1≤i≤kt, where kt is the number of nondominated configurations in period t. Each edge in the network stands for a possible transition from the previous period to the next period. Edges standing for transitions with too low similarity are truncated.We define two kinds of parameters in the network. The first one is the weight of edgeet,t+1ij, denoted as wt,t+1ij, set as the inverse of the similarity of configuration Cti and Ct+1j. The second one is the values ofObj1 toObj3 of node Cti, denoted as fti=ft,1i,ft,2i,ft,3i.Definition 7.
One defines a route in the forward network from one node to any node in the last period as a configuration link. A link without its starting node is called a path.For nodeCti, one denotes the cost along a path pi,tp from it to any node in the last period as gpi,tp=g1pi,tp,g2pi,tp,g3pi,tp,g4pi,tp, where g1pi,tp is the sum of similarities along the path and the other three items correspond to the sum of values of three objective functions of nodes along the path. Hence, one can define the dominating relationship between two paths.Definition 8.
Assuming two paths form the same time periodt to the last period as pi,tp and pj,tq, pi,tp≺pj,tq, if gpi,tp≤gpj,tq and ∃m,1≤m≤4, gmpi,tp<gmpj,tq.From Definition8, we can define a nondominated set of paths from node Cti, denoted as NDti. Then, we assume another path from one node in previous period, Ct-1j, which goes across Cti as pj,t-1r, and the part after Cti along pj,t-1r is pi,tp. Then, we can get the following relationships:(25a)pj,t-1r=pi,tp∪et-1,tji,Cti,(25b)g1pj,t-1r=g1pi,tp+wt-1,tji,(25c)gmpj,t-1r=gmpi,tp+ft,m-1i,m=2,3,4.Hence, we obtain the following proposition.Proposition 9.
Assuming a nondominated set of paths originating fromCt-1j and paths in this set go across Cti, the part after Cti along such path belongs to NDti.Proof.
Assume two paths originating fromCt-1j and going across Cti as pj,t-1p and pj,t-1q. Also, assume their parts after Cti as pi,tr and pi,ts, respectively. Considering the calculations in (25b) and (25c), if pi,tr≺pi,ts, then pj,t-1p≺pj,t-1q.Proposition9 shows that the subpaths not being in NDti can be discarded when we calculate NDt-1j. We consider all nodes, Cti, 1≤i≤kt, in period t, which can be transitioned from Ct-1j. We denote the nondominated paths originating from Ct-1j and going across Cti as NDt-1,tji, and we get (26)ND⌣t-1j=⋃iNDt-1,tji.We elite out nondominated paths in ND⌣t-1j and get the nondominated set of paths originating from Ct-1j, denoted as NDt-1j.We have presented the procedure to obtainNDt-1j so far when NDti is known. This is similar to one of the steps in typical dynamic programming except multiple objectives. In fact, theprinciple of optimality in monoobjective DP also works in multiobjective cases [39] and Proposition 9 just acts the same as theprinciple of optimality.Hence, we apply a backward multiobjective DP algorithm starting at stageT-1, getting the nondominated path set for each node in period T-1. In stage n, 1≤n≤T-1, we deal with edges from period t to t+1 and the nodes in period t+1 to generate sets of nondominated paths from nodes in period t.After the calculation in stage 1, we generate nondominated path sets from each node in period 1, denoted asND11,ND12,…,ND1k1. However, values of objective functions of configurations in period 1 are not involved yet. Next, we integrate nodes in period 1 into the paths to get final nondominated links.For each configuration in period 1,C1i, 1≤i≤k1 and corresponding ND1i, suppose a path in ND1i as pi,1p and the link li,1p that integrates C1i with pi,1p; then(27a)li,1p=pi,1p∪C1i,(27b)gmli,1p=gmpi,1p+f1,m-1i,m=2,3,4.We take operations like (27a) and (27b) for all nondominated paths in ND1i, 1≤i≤k1. Then, we classify the generated links into a set FND⌣. Finally, we select nondominated individuals from FND⌣ to obtain final nondominated configuration link set FND.


## 6. Experiment and Result Analysis

### 6.1. Setting

We test our methods with practical navigation database, recorded radar track, and practical flight plan teletext on July 20th, 2014. The target airspace is the high-altitude airspace controlled by Beijing Area Control Center (as in Figure1). Graph G representing the airspace is made up of 117 vertices and 293 edges (239 edges represent air routes), which is shown in Figure 11.Figure 11
The structure of the undirected graph.The MDFB is set as 15 nautical miles and the minimum staying time in one sector is set as 7 minutes. The threshold of instantaneous flow,Threinsflow, is set as 15. Other two thresholds in the memetic local search are set as FDup=5400 and FDlow=2500. Weights in (9) are determined by interviews with air traffic controllers, wFD=5, wac=1, wsc=1, wcc=1. Besides, we divide one day into 10 time periods, 0 a.m. to 6 a.m. is set as the first period considering its sparse traffic. After that, every two hours are set as one period such as 10 a.m. to 12 a.m. In the evolution algorithm, we set 100 individuals in one generation and the maximum times of iteration are set as 500. All simulations are executed on a PC with Core-i3 2120 CPU and 4 GB RAM.


### 6.2. Configuration Examples

To save space, we only demonstrate examples in four periods, 12:00–14:00, 16:00–18:00, 20:00–22:00, and 22:00–24:00. The period between 20:00–22:00 is the busiest period according to the recorded data.We show the four nondominated configurations in the nondominated configuration link with the largest gross similarity (see Section6.7) in Figure 12. The airspace is split into 9, 10, and 12 sectors in 12:00–14:00, 16:00–18:00, and 20:00–22:00 along with the increasing flow density. Sector number in the last period is also 12. However the traffic flow in this period is not so dense. In fact, at the beginning of period 22:00–24:00, the traffic continues the situation in 20:00–22:00, which is very dense and 12 sectors are necessary. The density gradually decreases and becomes much sparser near the midnight. However, as the configuration in one time period is fixed and it has to cover the busiest situation in this period because of the instantaneous flow constraint, the sector number during period 22:00–24:00 is also 12. If we can identify the flow density change more accurately, the pattern of time period division could be better. This is another critical and important problem in the DAC but not the scope of this paper.Figure 12
Configuration examples by the proposed approach.
(a) (b) (c) (d)In almost all periods except the period 20:00–22:00 with huge high density, northeastern, northern, and northwestern regions of the airspace are almost served by at most two sectors. This corresponds to the fact that the economy in northeast and northwest China is not so developed and flights from these regions are limited. On the contrary, in the southeastern and southern regions of the airspace, configurations changes frequently between two adjacent periods and much more sectors are needed during busy periods. In some situations, small sectors that even contain only one vertex are produced around the main stem of China’s north-south air route, “A461.” The peak flow density along “A461” occurs because flights from Guangzhou, Changsha, Wuhan, Zhengzhou and other big cities to Beijing are all probable to fly along “A461.” Hence, we have to limit the control area of any single sector so as to limit the controllers’ workload during high flow density periods.


### 6.3. Comparisons

In Table1, we summarize the majority of the metrics of configurations using six different methods, respectively, including current configuration in practice, the configuration generated only by the NSCA-based method and four evolution-based methods, that is, RI, CI, NI, and MI. We compare these results in terms of values of the three objective functions (Obj1 toObj3) and violating times of the three constraints (Con1 toCon3). “Nos” in Table 1 stands for the number of sectors. Density gap stands for the gap of the maximum and the minimum count of trajectory points in different sectors. For last four evolution-based methods, the values from the third column to the values in the fifth column are all averaged values of nondominated individuals in the final generation; values from the sixth to the eighth column and values in the column “density gap” are the biggest related values in the nondominated set. Column “Nos” gives the smallest sector number of single configuration in the nondominated set. Column “RT” gives the running time of evolution-based methods.


#### 6.3.1. Comparison with Current Configuration in Practice

In Table1, we notice the reductions in values of all three objective functions as well as the reduction of the density gap from configuration in practice to those by our approach. In addition, there are several occurrences of constraint violations in current configuration while any kind of violations are absent in our approach. Furthermore, much less sectors are needed in the proposed method comparing to the configuration in practice.


#### 6.3.2. Approximation of the NSCA

The values in rows named with “Only NSCA” show that the NSCA is just an approximation of the graphN-cut. Even if we do not consider constraint violations, the configuration obtained by only applying the NSCA shows no superiority to the current configuration in practice, especially from the aspect of balancing intrasector workloads. However, “Only NSCA” often performs well in terms ofObj2 orObj3. Hence, we may expect “nice genes” from the NSCA-based method and this is why we execute NI and MI to produce the initial generation.


#### 6.3.3. Comparison of Different Initial Generations

For every period, we set four cases about the initial generation of the CNSGA-II and list related results in last four rows of that period in Table1. Solutions obtained in RI are the worst considering even no feasible individuals can be produced (We have tried much more iterations but no feasible individual can be produced yet.). All individuals in the final generation in CI, NI, and MI are feasible, but there is no obvious dominance considering the averaged values of three objective functions. Although in the period between 12:00–14:00, the averaged value in CI dominates the averaged results in NI and MI, and several individuals in NI and MI are still not dominated by individuals in CI (see Figure 13). Similar situations occur in other periods. That is, if we unite the nondominated sets by CI, NI, and MI together and select once again, the final approximate Pareto Front is always made up by solutions from all these three cases.Figure 13
Nondominated sets of the period 12:00–14:00.


### 6.4. Final Nondominated Set

From the discussion in Section6.3, we propose the idea of parallel computation with CI, NI, and MI when executing the evolution algorithms. Then we merge their final generation together and select nondominated individuals from this merged set so as to get the final solutions of the 2D ERSP in one time period.


### 6.5. Applicability of the Automatic Boundary Depiction Method

The most important factor that affects the applicability of the automatic boundary depiction method is the conciseness of the depicted boundaries. If the boundaries are composed of too much segments or even jagged, the method is undesirable. We randomly select 10 optimal configuration links (that correspond to 100 configurations and 8500 sectors in total) and draw the boundaries using the proposed boundary depiction method. 59284 segments are generated in total and the boundary of single sector is made up by 6.97 segments in average. Among the 59284 segments, only 34 are produced by theOABD. TheOABD is called 17 times, which means that one segment in the concave hull is almost replaced by a two-part polygonal line in theOABD. This phenomenon may attribute to the advanced graph structure proposed in this paper that involves edges representing “Vonoroi neighborhood.” With these newly added edges, the distance between two adjacent vertices on the concave hull is limited so that the line segment connecting them seldom walks through too much “obstacle circles.” In all, the limited average number of segments consisting of the boundary of single sector and the fact that no jagged boundaries are produced by theOABD show the applicability of our automatic boundary depiction method.


### 6.6. Running Time

Because the sectorization problem is not a real-time decision problem and almost all solving methods need iterative optimization, few previous studies have reported the running time of their methods. However, in our mind, the running time is also important because short running time is helpful in utilizing latest traffic information.


#### 6.6.1. Time for Data Preparation

In Section3.3, we propose multiple data forms to store the traffic information. In our simulation, time spent to obtain traffic information for single time period ranges from 169 seconds to 324 seconds. The heavier the traffic, the more the time needed to calculate the metrics.


#### 6.6.2. Time for Optimization

The values in the last column of Table1 show that the evolution algorithm-based optimization is around 15 minutes for any single time period.


#### 6.6.3. Time for Boundary Depiction

The time spent to draw the sector boundaries for single configuration is about 1.5 minutes. The most time-consuming part is judging whether segments on boundary polygons intersect with inner-edges of sectors or whether they intrude obstacle circles. Besides, as we do not need to depict the boundaries during the optimization process and the depiction is only needed once after the evolution, the 1.5 min running time is acceptable and it is expected to have advantages over methods with manual operations.In all, in our simulation, it takes about 20 minutes from the data preparation to the sector boundary depiction for single configuration during a single time period. This implies if we want to get a configuration for next 2 hours, we can start just about 20 minutes before the starting moment of the period, which will benefit our approach with latest air traffic information so that to minimize the impact of uncertainty. In addition, as the approaches are all executed on a PC with Core-i3 CPU, we can expect much shorter running time if we use a computer with more computing power.


### 6.7. Results of the Nondominated Configuration Links

We randomly select 20 nondominated solutions (configurations) for each time period; that is, the number of links of configurations in the forward network is about 2010. Among them, we get 5284 nondominated links. One of the links with the largest gross similarity is C113→C222→C37→C497→C556→C6143→C7215→C897→C9112→C109. In fact, we have demonstrated the results of C556, C7215, C9112, and C109 in Figure 12.We must point out that these planned links can serve as candidates and they however may not be employed definitely. Air traffic management officers can choose one such link at the beginning of a day in order to start some necessary preparations for subsequent periods, such as determining the schedule of controllers and previewing the sectors. Although exact configurations applied may change because of uncertainties in the air traffic system, such changes should not be too serious if no emergency occurs.


## 6.1. Setting

We test our methods with practical navigation database, recorded radar track, and practical flight plan teletext on July 20th, 2014. The target airspace is the high-altitude airspace controlled by Beijing Area Control Center (as in Figure1). Graph G representing the airspace is made up of 117 vertices and 293 edges (239 edges represent air routes), which is shown in Figure 11.Figure 11
The structure of the undirected graph.The MDFB is set as 15 nautical miles and the minimum staying time in one sector is set as 7 minutes. The threshold of instantaneous flow,Threinsflow, is set as 15. Other two thresholds in the memetic local search are set as FDup=5400 and FDlow=2500. Weights in (9) are determined by interviews with air traffic controllers, wFD=5, wac=1, wsc=1, wcc=1. Besides, we divide one day into 10 time periods, 0 a.m. to 6 a.m. is set as the first period considering its sparse traffic. After that, every two hours are set as one period such as 10 a.m. to 12 a.m. In the evolution algorithm, we set 100 individuals in one generation and the maximum times of iteration are set as 500. All simulations are executed on a PC with Core-i3 2120 CPU and 4 GB RAM.


## 6.2. Configuration Examples

To save space, we only demonstrate examples in four periods, 12:00–14:00, 16:00–18:00, 20:00–22:00, and 22:00–24:00. The period between 20:00–22:00 is the busiest period according to the recorded data.We show the four nondominated configurations in the nondominated configuration link with the largest gross similarity (see Section6.7) in Figure 12. The airspace is split into 9, 10, and 12 sectors in 12:00–14:00, 16:00–18:00, and 20:00–22:00 along with the increasing flow density. Sector number in the last period is also 12. However the traffic flow in this period is not so dense. In fact, at the beginning of period 22:00–24:00, the traffic continues the situation in 20:00–22:00, which is very dense and 12 sectors are necessary. The density gradually decreases and becomes much sparser near the midnight. However, as the configuration in one time period is fixed and it has to cover the busiest situation in this period because of the instantaneous flow constraint, the sector number during period 22:00–24:00 is also 12. If we can identify the flow density change more accurately, the pattern of time period division could be better. This is another critical and important problem in the DAC but not the scope of this paper.Figure 12
Configuration examples by the proposed approach.
(a) (b) (c) (d)In almost all periods except the period 20:00–22:00 with huge high density, northeastern, northern, and northwestern regions of the airspace are almost served by at most two sectors. This corresponds to the fact that the economy in northeast and northwest China is not so developed and flights from these regions are limited. On the contrary, in the southeastern and southern regions of the airspace, configurations changes frequently between two adjacent periods and much more sectors are needed during busy periods. In some situations, small sectors that even contain only one vertex are produced around the main stem of China’s north-south air route, “A461.” The peak flow density along “A461” occurs because flights from Guangzhou, Changsha, Wuhan, Zhengzhou and other big cities to Beijing are all probable to fly along “A461.” Hence, we have to limit the control area of any single sector so as to limit the controllers’ workload during high flow density periods.


## 6.3. Comparisons

In Table1, we summarize the majority of the metrics of configurations using six different methods, respectively, including current configuration in practice, the configuration generated only by the NSCA-based method and four evolution-based methods, that is, RI, CI, NI, and MI. We compare these results in terms of values of the three objective functions (Obj1 toObj3) and violating times of the three constraints (Con1 toCon3). “Nos” in Table 1 stands for the number of sectors. Density gap stands for the gap of the maximum and the minimum count of trajectory points in different sectors. For last four evolution-based methods, the values from the third column to the values in the fifth column are all averaged values of nondominated individuals in the final generation; values from the sixth to the eighth column and values in the column “density gap” are the biggest related values in the nondominated set. Column “Nos” gives the smallest sector number of single configuration in the nondominated set. Column “RT” gives the running time of evolution-based methods.


### 6.3.1. Comparison with Current Configuration in Practice

In Table1, we notice the reductions in values of all three objective functions as well as the reduction of the density gap from configuration in practice to those by our approach. In addition, there are several occurrences of constraint violations in current configuration while any kind of violations are absent in our approach. Furthermore, much less sectors are needed in the proposed method comparing to the configuration in practice.


### 6.3.2. Approximation of the NSCA

The values in rows named with “Only NSCA” show that the NSCA is just an approximation of the graphN-cut. Even if we do not consider constraint violations, the configuration obtained by only applying the NSCA shows no superiority to the current configuration in practice, especially from the aspect of balancing intrasector workloads. However, “Only NSCA” often performs well in terms ofObj2 orObj3. Hence, we may expect “nice genes” from the NSCA-based method and this is why we execute NI and MI to produce the initial generation.


### 6.3.3. Comparison of Different Initial Generations

For every period, we set four cases about the initial generation of the CNSGA-II and list related results in last four rows of that period in Table1. Solutions obtained in RI are the worst considering even no feasible individuals can be produced (We have tried much more iterations but no feasible individual can be produced yet.). All individuals in the final generation in CI, NI, and MI are feasible, but there is no obvious dominance considering the averaged values of three objective functions. Although in the period between 12:00–14:00, the averaged value in CI dominates the averaged results in NI and MI, and several individuals in NI and MI are still not dominated by individuals in CI (see Figure 13). Similar situations occur in other periods. That is, if we unite the nondominated sets by CI, NI, and MI together and select once again, the final approximate Pareto Front is always made up by solutions from all these three cases.Figure 13
Nondominated sets of the period 12:00–14:00.


## 6.3.1. Comparison with Current Configuration in Practice

In Table1, we notice the reductions in values of all three objective functions as well as the reduction of the density gap from configuration in practice to those by our approach. In addition, there are several occurrences of constraint violations in current configuration while any kind of violations are absent in our approach. Furthermore, much less sectors are needed in the proposed method comparing to the configuration in practice.


## 6.3.2. Approximation of the NSCA

The values in rows named with “Only NSCA” show that the NSCA is just an approximation of the graphN-cut. Even if we do not consider constraint violations, the configuration obtained by only applying the NSCA shows no superiority to the current configuration in practice, especially from the aspect of balancing intrasector workloads. However, “Only NSCA” often performs well in terms ofObj2 orObj3. Hence, we may expect “nice genes” from the NSCA-based method and this is why we execute NI and MI to produce the initial generation.


## 6.3.3. Comparison of Different Initial Generations

For every period, we set four cases about the initial generation of the CNSGA-II and list related results in last four rows of that period in Table1. Solutions obtained in RI are the worst considering even no feasible individuals can be produced (We have tried much more iterations but no feasible individual can be produced yet.). All individuals in the final generation in CI, NI, and MI are feasible, but there is no obvious dominance considering the averaged values of three objective functions. Although in the period between 12:00–14:00, the averaged value in CI dominates the averaged results in NI and MI, and several individuals in NI and MI are still not dominated by individuals in CI (see Figure 13). Similar situations occur in other periods. That is, if we unite the nondominated sets by CI, NI, and MI together and select once again, the final approximate Pareto Front is always made up by solutions from all these three cases.Figure 13
Nondominated sets of the period 12:00–14:00.


## 6.4. Final Nondominated Set

From the discussion in Section6.3, we propose the idea of parallel computation with CI, NI, and MI when executing the evolution algorithms. Then we merge their final generation together and select nondominated individuals from this merged set so as to get the final solutions of the 2D ERSP in one time period.


## 6.5. Applicability of the Automatic Boundary Depiction Method

The most important factor that affects the applicability of the automatic boundary depiction method is the conciseness of the depicted boundaries. If the boundaries are composed of too much segments or even jagged, the method is undesirable. We randomly select 10 optimal configuration links (that correspond to 100 configurations and 8500 sectors in total) and draw the boundaries using the proposed boundary depiction method. 59284 segments are generated in total and the boundary of single sector is made up by 6.97 segments in average. Among the 59284 segments, only 34 are produced by theOABD. TheOABD is called 17 times, which means that one segment in the concave hull is almost replaced by a two-part polygonal line in theOABD. This phenomenon may attribute to the advanced graph structure proposed in this paper that involves edges representing “Vonoroi neighborhood.” With these newly added edges, the distance between two adjacent vertices on the concave hull is limited so that the line segment connecting them seldom walks through too much “obstacle circles.” In all, the limited average number of segments consisting of the boundary of single sector and the fact that no jagged boundaries are produced by theOABD show the applicability of our automatic boundary depiction method.


## 6.6. Running Time

Because the sectorization problem is not a real-time decision problem and almost all solving methods need iterative optimization, few previous studies have reported the running time of their methods. However, in our mind, the running time is also important because short running time is helpful in utilizing latest traffic information.


### 6.6.1. Time for Data Preparation

In Section3.3, we propose multiple data forms to store the traffic information. In our simulation, time spent to obtain traffic information for single time period ranges from 169 seconds to 324 seconds. The heavier the traffic, the more the time needed to calculate the metrics.


### 6.6.2. Time for Optimization

The values in the last column of Table1 show that the evolution algorithm-based optimization is around 15 minutes for any single time period.


### 6.6.3. Time for Boundary Depiction

The time spent to draw the sector boundaries for single configuration is about 1.5 minutes. The most time-consuming part is judging whether segments on boundary polygons intersect with inner-edges of sectors or whether they intrude obstacle circles. Besides, as we do not need to depict the boundaries during the optimization process and the depiction is only needed once after the evolution, the 1.5 min running time is acceptable and it is expected to have advantages over methods with manual operations.In all, in our simulation, it takes about 20 minutes from the data preparation to the sector boundary depiction for single configuration during a single time period. This implies if we want to get a configuration for next 2 hours, we can start just about 20 minutes before the starting moment of the period, which will benefit our approach with latest air traffic information so that to minimize the impact of uncertainty. In addition, as the approaches are all executed on a PC with Core-i3 CPU, we can expect much shorter running time if we use a computer with more computing power.


## 6.6.1. Time for Data Preparation

In Section3.3, we propose multiple data forms to store the traffic information. In our simulation, time spent to obtain traffic information for single time period ranges from 169 seconds to 324 seconds. The heavier the traffic, the more the time needed to calculate the metrics.


## 6.6.2. Time for Optimization

The values in the last column of Table1 show that the evolution algorithm-based optimization is around 15 minutes for any single time period.


## 6.6.3. Time for Boundary Depiction

The time spent to draw the sector boundaries for single configuration is about 1.5 minutes. The most time-consuming part is judging whether segments on boundary polygons intersect with inner-edges of sectors or whether they intrude obstacle circles. Besides, as we do not need to depict the boundaries during the optimization process and the depiction is only needed once after the evolution, the 1.5 min running time is acceptable and it is expected to have advantages over methods with manual operations.In all, in our simulation, it takes about 20 minutes from the data preparation to the sector boundary depiction for single configuration during a single time period. This implies if we want to get a configuration for next 2 hours, we can start just about 20 minutes before the starting moment of the period, which will benefit our approach with latest air traffic information so that to minimize the impact of uncertainty. In addition, as the approaches are all executed on a PC with Core-i3 CPU, we can expect much shorter running time if we use a computer with more computing power.


## 6.7. Results of the Nondominated Configuration Links

We randomly select 20 nondominated solutions (configurations) for each time period; that is, the number of links of configurations in the forward network is about 2010. Among them, we get 5284 nondominated links. One of the links with the largest gross similarity is C113→C222→C37→C497→C556→C6143→C7215→C897→C9112→C109. In fact, we have demonstrated the results of C556, C7215, C9112, and C109 in Figure 12.We must point out that these planned links can serve as candidates and they however may not be employed definitely. Air traffic management officers can choose one such link at the beginning of a day in order to start some necessary preparations for subsequent periods, such as determining the schedule of controllers and previewing the sectors. Although exact configurations applied may change because of uncertainties in the air traffic system, such changes should not be too serious if no emergency occurs.


## 7. Conclusion and Future Work

This study solves the 2D ERSP with approaches based on weighted undirected graph cuts. We comprehensively considered the multiple objectives and constraints of the 2D ERSP. We employed memetic local search-embedded CNSGA-II to evolve the initial clusters and proposed several methods to produce initial generation. We also proposed a concave hull-based method to draw exact sector boundaries. In addition to the sectorization problem, we also solved the CTP, which generated a series of nondominated links of configurations for several consecutive time periods.The simulation results showed significant improvement from the configuration applied in practice in almost all important aspects of the sectorization problem. Along with the results in depicting boundaries and in the CTP, we can claim that the proposed approaches are promising practical decision support tools for DAC considering their comprehensiveness, practicality, automatization, and efficiency.However, there are still several important issues we should address in the future:(1) Other Automatic Boundary-Depicting Algorithms. As we stated in Section 1, powerful boundary-depicting algorithms for undirected graph cut-based approaches are very limited. Searching for more efficient depicting algorithms is very important.(2) Interaction between Traffic Flow Management (TFM) and DAC. Both TFM and DAC are key problems in air traffic management, and their outputs seem to act as inputs to each other. Can we find a logic to deal with them simultaneously? Research in this direction is already underway, but no general architecture or principle has been proposed.(3) Local Adjustment of Airspace Configuration. We mentioned in Section 6.7 that we can plan a configuration link at the beginning of a day and follow the link if no emergency occurs. However, when an emergency occurs that sharply changes the traffic flow pattern, we must adjust the configuration temporarily. This local adjustment must be fulfilled very soon after the emergency; therefore, algorithms dealing with local adjustment must be very efficient.



---

*Source: 6048326-2016-06-26.xml*