Multi-Attributes Decision-Making for CDO Trajectory Planning in a Novel Terminal Airspace

: Continuous Descent Operations (CDO) has been recognized as an effective way to signiﬁcantly reduce fuel burn and noise impact. Designing efﬁcient and ﬂexible arrival routes for generating conﬂict-free and economical trajectories is a cornerstone for fully achieving CDO by high-level automation in high-density trafﬁc scenarios. In this research, inspired by the Point Merge (PM), we design the Inverted Crown-Shaped Arrival Airspace (ICSAA) and its operational procedures to support Omni-directional CDO. In order to generate optimal conﬂict-free trajectories for upcoming aircraft in an efﬁcient manner, we established a multi-objective trajectory optimization model solved by Non-dominated Sorting Genetic Algorithm with Elitist Strategy (NSGA-II). The Pareto solutions of minimal fuel consumption and trip time were achieved in single aircraft and highly complex multi-aircraft scenarios. Among all the elements of Pareto front, we obtained an unique solution with Entropy-Weights Method and TOPSIS (Technique for Order Preference by Similarity to an Ideal Solution) to strike a better trade-off among collision probability, fuel consumption, and trip time, which incorporates both air trafﬁc controller’s and pilot’s interests. The effectiveness of CDO performance improvement and computational efﬁciency in different scenarios were veriﬁed. The ICSAA would be a promising structure that promotes the application of automated and ﬂexible CDO.


Introduction
Worldwide Air Traffic Management (ATM) system is undergoing the process of upgrading and transformation, in order to cope with increasing air traffic demand and congestion, as well as diverse expectations of stakeholders for safety, efficiency, economy, and pro-environment, especially in high-density airports and terminal airspace [1]. Trajectory-Based Operation (TBO) and Performance-Based Operation (PBO) have been identified as the core concept of the future ATM system [2].
Continuous Descent Operation (CDO), one of the key elements in TBO and PBO, has been recognized as an effective procedure that may improve operation efficiency and environmental benefits in terminal airspace, which was initially designed to abate noise [3]. A series of CDO flight trials validated that CDO could significantly reduce fuel burn and noise impact during arrival phase by keeping arriving aircraft at their cruising altitude for longer and then executing continuous descent with no level-flight segments [4]. In 2010, the International Civil Aviation Organization (ICAO) published Doc 9931 Continuous Descent Operation Manual, which provides guidelines for CDO procedure design [5]. At present, CDO has become one of the building blocks for Global Air Navigation Plan (GANP), Single European Sky ATM Research (SESAR), and Next Generation Air Transportation System (NextGen).
Previous studies demonstrated that some factors, like vertical profile [6], speed profile [6], flight path angle [7], meteorology [8], sequencing and flight scheduling [9], and capacity [10] have great effects on the performance of CDO. In view of such factors, pre-tactical CDO trajectory optimization with various constraints and objectives in singleaircraft or multi-aircraft scenarios attracted great attention [11,12]. In most of the literature, the objectives of trajectory optimization include fuel consumption [13], flight time [14], emission [15] and noise impact [16], etc. A few scholars also assessed trajectories based on optimization in terms of safety, efficiency, airspace capacity, and ecological compatibility [14,17]. Samà et al. optimized CDO trajectories in a busy Terminal Maneuvering Area (TMA) in terms of flight time and fuel consumption using lexicographic method to make a decision among alternative approaches, by presetting the primary and secondary performance indicator [14]. It was essentially a greedy strategy with subjectively prioritizing the importance of objectives. In addition, real-time CDO trajectory control strategies have been investigated in the presence of abnormal operation deviating from the planned path [18], as well as adverse meteorological factors, such as storm [19] and gust [20].
With the deepening of studies, it suggested that traditional Standard Terminal Arrival Route (STAR) is too rigid to exploit the potential of the CDO benefits. Alam et al. established a new type of transition airspace to generate dynamic CDO trajectories in order to achieve lateral optimization [21]. However, it is difficult to execute path recovery once the approach aborted. Furthermore, in Reference [21], the aircraft has to adjust headings continuously in low-altitude area, which would increase the operational complexity for both controllers and pilots. From the practical perspective, it is suggested that, without advanced Air Traffic Control (ATC) automations or flexible flight procedures, 4D CDO trajectory operation that mostly relies on Air Traffic Controller's mental cognition could only be applied in some less busy airports or off-peak hours [12].
The Point Merge System (PMS) is a promising procedure to improve the performance of CDO, which was proposed by the EUROCONTROL Experimental Center (EEC) [22]. It comprises sequencing legs and a merge point, as shown in Figure 1. The PMS enables flights maintain altitude and fly along the sequencing legs, and then descend and fly to the merge point once a "direct-to" instruction is issued. The landing sequence and the longitudinal inter-aircraft separations are controlled during operation. Research has validated that the Point Merge procedure possesses higher flexibility and lateral predictability, which brings benefits to operator's workload, arrival efficiency [22], and environmental impact [23]. Nevertheless, the delay absorption ability of the sequencing legs in the PMS is constrained by their limited length. Once the sequencing legs became saturated or in some emergent situations, inbound flow should be vectored off from the PMS. Such abnormal maneuver would significantly increase the traffic complexity, as well as degrade operation safety and efficiency [24]. Therefore, designing a new CDO airspace system and automated 4D conflict-free trajectory planning, would be a practical breakthrough to improve the flexibility, controllability, and predictability of CDO operation. In this paper, inspired by PMS, we propose a novel terminal airspace structure named Inverted Crown-Shaped Arrival Airspace (ICSAA) that supports Omni-directional arrival.
Based on such an airspace system, the conflict-free, economical, and efficient CDO trajectories are generated and assigned to each aircraft automatically using Multi-Objective Decision-Making(MODM) methods. Simulation experiments validate overall performance of optimized trajectories within the ICSAA in terms of collision probability, fuel consumption, and trip time. The main work of this paper are illustrated as follows: • The design of ICSAA with intuitional operation modes: In order to increase the elasticity of traditional STAR, and the capacity of sequencing legs in the PMS, a flexible terminal airspace called ICSAA was designed to merge traffic streams from all the directions and accommodate more aircraft using circular design. Based on ICSAA, two concise and intuitive operation modes: Mode H and Mode V, were designed to ensure the predictability, safety, and transparency of automated CDO operations. • Pareto optimal front generation for Conflict-free CDO trajectories in the ICSAA: Inside the ICSAA, a multi-objective optimization model that aims at generating conflict-free CDO trajectories with minimal fuel consumption and trip time is proposed by simulating flight dynamics. By comparing with some state-of-the-art multiobjective algorithms, the Non-dominated Sorting Genetic Algorithm with Elitist Strategy (NSGA-II) shows its best performance in convergence near the true Pareto optimal front efficiently, thus is selected to solve the proposed optimization problem. • Multi-attribute decision-making for optimal CDO trajectories selection: In order to select the optimal trajectories among the Pareto optimal front, from the perspective of both pilots' and controllers' preferences, we additionally take collision risk as an attribute besides fuel consumption and trip time. The entropy-based TOPSIS method is adopted to select an unique solution by multi-attribute decision-making strategies.
To validate the proposed model and algorithm, the overall performance in singleaircraft, low-density and high-density scenarios are investigated. The results verified that proposed automated CDO trajectory planning in the ICSAA are of high efficiency and strike a better trade-off in terms of economy, efficiency, and safety.
The remainder of this paper is organized as follows: We first illustrate the structure and operational procedures of the ICSAA in Section 2; then, the descriptions are made to explain the detailed trajectory planning problems to be solved in the novel airspace in Section 3. In order to derive optimization objectives, the aircraft performance model and fuel consumption model are established in Sections 4 and 5, respectively, followed by a multi-objective optimization model in Section 6 and multi-attribute decision-making strategies in Section 7. Section 8 shows the numerical experiments of verifying the effectiveness of proposed methods in different scenarios. Conclusion and the future work are discussed in Section 9.

Basic Structure of the Inverted Crown-Shaped Arrival Airspace
According to Doc 9931 Continuous Descent Operation Manual, CDO is an aircraft operating technique aided by appropriate airspace and procedure design and reliable ATC clearances enabling the execution of an optimized trajectories with low engine thrust settings and, where possible, a low drag configuration, thereby reducing fuel burn and emissions during descent. Based on that ideology and enlightened by PMS, we proposed a novel flexible terminal airspace structure and operation procedure, named the ICSAA.
In the ICSAA, we still instruct flights with "direct to" instructions instead of radar vectoring, which could reduce the pilot-ATC communication. Given that the length and number of sequencing legs impose constraints on the performance of the PMS, we extend the arc-shaped sequencing leg to a circle. As shown in Figure 2, we define the artificial terminal area as a series of concentric sequencing rings with the merge point in the center. Each ring has preset waypoints spaced 10 • apart from each other along the circumference to ensure safety separation between arriving aircraft. The artificial ICSAA has five levels of rings at different altitude. The altitude of the outermost ring is 6000 m, on which waypoints represent the entry points, while the innermost ring is at 4800 m, on which each waypoint represents the beginning of the CDO trajectory. The vertical separation between adjacent rings is 300 m. It is worth noting that, for the fifth-level ring, there are three layers at the same altitude. Departing from different waypoints on different layers, aircraft could fly to the merge point with different flight path angles ranged from 3 • to 4 • . In addition, to assure a smooth and efficient interception to the glide slope, the height of the merge point and the central angle of the fan-shaped envelope is set as 900 m and 120 • , respectively. For the sake of simplicity, we numbered the waypoints along each ring, as shown in Figure 2b. Please note that, at the fifth level, the waypoints on the outermost layer, the second layer, and the third layer are numbered from 1 to 13, 14-26, and 27-39, respectively.

The Operational Procedure
Based on flight dynamics of aircraft in the air, two switchable operational modes are defined in the ICSAA: Mode H and Mode V, a multiphase mixed-integer optimal control approach. Definition 1 (Mode H). The aircraft is considered as flying into a horizontal plane in a clockwise direction (or anti-clockwise direction). Thus, flight path angle γ, the derivative of flight path anglė γ and the derivative of heightḣ are set to 0. In this mode, the equation also holds: Lcosφ = mg, where L is lift; φ is bank angle; m is the mass of aircraft, and g is the gravitational acceleration.
Definition 2 (Mode V). The aircraft is considered as flying into the vertical plane, performing a leveled-wing descent. Therefore, flight path angle γ, the derivative of flight path angleγ and the derivative of heightḣ are not equal to 0, while the bank angle φ is set to 0.
Furthermore, to assure trajectory predictability and operation order, aircraft should descend level by level and cannot go through the interior of ICSAA. In other words, the aircraft must always fly on the wall of the ICSAA. Therefore, we define Separation Z. If and only if the separation Z is less than a certain value, the aircraft could descend with Mode V; otherwise, the aircraft had to fly along the circumference of rings in a clockwise direction (or anti-clockwise direction) with Mode H.

Definition 3 (Separation Z).
Suppose that the present position of aircraft is at the waypoint i(x i , y i , h i ), for any waypoint j(x j , y j , h j ) on the adjacent inner ring, which the aircraft is going to fly to, the Separation Z is presented in Figure 2b and defined by Z(i, j) = (x i − x j ) 2 + (y i − y j ) 2 .
In the multi-aircraft scenario, the situation would be more complex. Not only should Separation Z be less than the certain value but also no potential conflicts exist during descending, when the aircraft could start to fly to the next inner ring. Otherwise, it would fly with Mode H until meeting descent conditions. Obviously, when flying along the circumference of rings, aircraft should keep the safety separation with preceding aircraft by adjusting speed.
Compared with traditional PMS, thanks to the circular design, the ICSAA could integrate inbound flow from all directions, and the sequencing legs (i.e., the rings of ICSAA) would accommodate more aircraft. The proposed intuitive flight modes also possess positive potentials in reducing operational complexity while simplifying automated conflict-free trajectory planning. To some extent, this design also renders more flexible segregation of the inbound and outbound flow. Furthermore, the fifth-level rings provide more continuous descending options for aircraft with different flight path angles, which would significantly affect the fuel consumption and trip time. Therefore, the proposed ICSAA has higher structural and procedural flexibility that would provide larger solution space for CDO trajectory optimization.

Problem Definition
In order to search optimal CDO trajectories, the airspace is modeled as a network G = (V, E), and each vertex (i, j) represents waypoint of ICSAA. The adjacency matrix is constructed, of which the element, Θ = e[(i n , j n ), (i m , j m )], represents connectivity of the directed network, in accordance with the operation rules of ICSAA mentioned above. If the edge, (i n , j n ) → (i m , j m ), is connected, e[(i n , j n ), (i m , j m )] = 1; otherwise, e[(i n , j n ), (i m , j m )] = + ∝. Therefore, the path P ath of the aircraft α is expressed as a set of connected vertices in network G. The problem we studied can be stated as follows: based on the expected entry points and expected arrival time of aircraft, we would generate optimal CDO trajectories automatically in the ICSAA, aiming at achieving best performance in economy, efficiency, and safety on the promises of acceptable computation time, to support the pre-tactic operation of air traffic management.

Aircraft Performance Model
Aircraft performance model is used to control the flight dynamics in accordance with the flight modes in the ICSAA. Glover and Lygeros [25] came up with the Point Mass Model (PMM), derived from the Newtonian dynamics. EUROCONTROL [26] proposed the Total-Energy Model (TEM) in Base of Aircraft Data (BADA), derived from the workenergy theorem for continuous descent. In this research, in order to simplify the trajectory optimization problem, a Point-Mass performance model, adapted from the work mentioned above, is adopted to simulate the CDO flight while measuring the fuel consumption and trip time. The model can be stated as follows: where the dotted terms are derivatives with respect to time; x, y, and h represent the position of aircraft in the three-dimension coordinate system; v is true airspeed; φ is bank angle; τ is attack angle; γ is flight path angle; ψ is course angle; W = (ω 1 , ω 2 , ω 3 ) ∈ R 3 is wind speed; T is engine thrust; D and L is the aerodynamic drag and lift, which are modeled, respectively, as: where ρ is air density; S is the wing reference area; and C D and C L are the drag coefficient and the lift coefficient, respectively, modeled as: where C D0 and C D2 are BADA coefficients. The drag coefficient C D is expressed as a function of the lift coefficient C L . When simulating the aircraft descent, we use the standard speed profile recommended by BADA and perform curve fitting for those data, as shown in Figure 3. The goodness of curve fitting R 2 = 0.9820. The result is stated as follows: In the above expressions, velocity is the function of height. In order to obtain the derivative of velocity versus time, we made the following deformations.
It is worth noting that, in order to improve the accuracy of the flight descent path modeling, some atmospheric properties (temperature, pressure, density) are expressed as a function of altitude, based on the BADA reversion 3.11 issued in 2013, rather than fixed values [27].
where R is real gas constant for air, the value of which is 287.05287 m 2 /(K · s 2 ); T em,0 is the standard atmospheric temperature at mean sea level, of which the value is 288.15 K; T em is temperature differential at mean sea level, which is the difference in atmospheric temperature at mean sea level between a given non-standard atmosphere and International Standard Atmosphere (ISA), and the value is set to 0; β is ISA temperature gradient with altitude below the tropopause, the value of which is −0.0065 K/m; p 0 is standard atmospheric pressure at mean sea level; and g is gravitational acceleration, the value of which is 9.80665 m/s 2 .
Based on the aircraft performance model (those equations are denoted as M), once the starting point (i s , j s ), ending point (i e , j e ), flight path P ath , and initial speed v 0 are determined, the real-time position of aircraft (x, y, h) at time t could be obtained following speed-altitude profile recommended by BADA, which can be expressed by the Equation (7). Such an equation enables the conflict detection of CDO trajectories.

Fuel Consumption Model
Thrust Specific Fuel Consumption (TSFC) model in BADA is applied for fuel consumption computation. For the jet engines, the thrust specific fuel consumption is specified as a function of the true airspeed: where C f 1 , C f 2 are the fuel flow coefficient, and C f cr is the fuel flow correction factor, the value of which are given by BADA; f c is the fuel flow and the fuel consumption, and F c is the integral of fuel flow with respect to time. In the network of ICSAA, the integral can be also approximately considered as a sum of discretized fuel consumption on each link along the whole flight path.

Decision Variables
In order to identify whether the aircraft α fly through the j th waypint of the i th ring in the ICSAA, the 0-1 variable, χ(α) (i,j) , is introduced, as shown in Formula (9): Furthermore, in order to indicate the connectivity among nodes determined by Formula (9), a variable is introduced to denote whether the aircraft chooses the edge between two nodes as a part of the CDO path. Assume two nodes are (i 1 , j 1 ) and (i 2 , j 2 ), respectively,

Objectives
Total fuel consumption and trip time that aircraft take from entry point to Final Approach Fix (FAF) are set as two objectives for the CDO trajectory optimization. Please note that the noise abatement here is not chosen as an objective in an artificial scenario be-cause the noise impact evaluation highly relies on the geographical features of the vicinity of airport [28].
where α names different aircraft; f (α) and t(α) are the fuel consumption and trip time, which aircraft α spent during flight from entry point to the FAF; and N(i, j) is introduced to represent the next waypoint connects with waypoint (i, j) in the planed CDO trajectory.

• Uniqueness Constraints
The trajectory of aircraft α should be unique and definite, that is, for any vertex (i p , j p ) along a trajectory, one and only one following waypoint directly connects to it.
• Descent Constraints In order to assure that aircraft could descend level by level, thus, the sum of χ(α) (i,j) in the each level should be equal to or greater than 1. •

Connection Constraints
The trajectory of each aircraft should be in compliance with operation mode of IC-SAA mentioned above. The edge [(i, j), N(i, j)] along the planned trajectory shall satisfy that e[(i, j), N(i, j)] = 1.
• Acceleration Rate The accelerate rate is a value that relates to the performance of aircraft. In general, in view of the comfort for passengers, the accelerate rate could not be too large.
The BADA manual recommends a maximum longitudinal acceleration: • Final Approach Speed When aircraft start to land, the velocity should be decreased to intercept the glide slope. But there is a minimum speed limit. The minimum landing speed v min is now determined as follows: where 1.3 is a factor recommended by the BADA manual for all aircraft operations, v stall is the stall speed at the reference mass m re f , and m is the simulated aircraft mass which is updated with fuel consumption.

• Safety Separations
In order to assure safety, the separation between aircraft shall be larger than the minimum. The position of aircraft α 1 and aircraft α 2 at time t could be obtained with the Equation (7), then we could get the horizontal separation at time t, S hor (t, α 1 , According to regulation, two aircraft are considered to be in a conflict if their horizontal separation is less than 10 km and vertical separation is less than 300 m. Therefore, during approach, C(t, α 1 , α 2 ) should always be 1.

Solution Algorithm for Generating Pareto Optimal Front
In order to solve the multi-objective optimization problem, we introduce the Nondominated Sorting Genetic Algorithm with Elitist Strategy (NSGA-II) algorithm, which has the best performance for this particular problem, compared with the state-of-the-art multi-objective optimization algorithms, as illustrated in Section 8.2. Adaptions to the critical genetic operators are made to generate the Pareto solutions [29,30]. The framework of the NSGA-II algorithm is illustrated in Figure 4. (1,1) , µ(α) (2,1) (1,1) , µ(α) (2,2) (1,1) , µ(α) (2,3) (1,1) , µ(α) (2,4) (1,1) , ...µ(α) Merge Point (5,39) }, while, in the multi-aircraft scenario, the chromosome is the splicing of chromosome for the each aircraft, that is, C h = {X(α 1 ), X(α 2 ), X(α 3 ), . . . , X(α n ), F, T ime , r, D cro }, where F and T ime are the fitness of the individual, defined as the reciprocal of objective functions, and r and D cro are the rank and crowding distance, respectively, which can be calculated based on Non-Dominated Sort introduced in the following part. Based on above coding schema, we can easily decode the chromosome into a trajectory of aircraft α. • Initial Population Initial population is the starting state of iteration for heuristic algorithm, which has a great effect on optimized results and computational efficiency. In the single aircraft scenario, the trajectory is optimized with a randomly generated initial population, while, in the multi-aircraft scenario, the initial population is generated by randomly combining the non-determined solutions derived in the single aircraft scenario. • Non-Dominated Sort and Crowding Distance Each chromosome of the NSGA-II corresponds to a feasible solution, by which we could get the trajectory of the aircraft inside the ICSAA through decoding. The population, composed of multiple chromosomes, is the set of feasible solutions. Based on the objective function, we could get the value of fuel consumption and trip time for each solution, as well as the fitness of the chromosomes. In order to achieve faster convergence, each population can be layered and sorted to create multiple Pareto fronts with different rankings based on fitness. Solutions with lower rankings (higher fitness) are preferred. If any two solutions belong to the same front, then the solution located in a less crowded region is regarded as a better one, by introducing crowding distance mechanism to improve diversity of the population. Obviously, solutions with larger crowding distance imply that there would be probably more unexplored solution space around them. Determined by aforementioned non-domination sort and crowding distance, the identified prospective individuals will be moved to the front of population and used to generate the next generation. • Genetic Operations In this research, the elite individual is selected with binary tournament selection and elitism strategy to create the offspring with crossover and mutation. To be specific, the crossover operator is implemented pairwise. The basic unit of crossover is the gene fragment, which could represent an entire trajectory of the aircraft. Suppose we pick out two individuals C h,1 and C h,2 in the current population randomly; then, two new individuals will be generated in the following way: where rp is a random integer between 1 and n. The above crossover operation is triggered by a crossover probability. Mutation is operated in a simple and intuitive way. For each of individuals, once the random number is greater than the mutation probability that we set beforehand, the individual starts to mutate at a selected gene randomly. After the mutation process, if the individual does not satisfy the constraints mentioned above, the fitness of the individual will be set to 0 and the individual would be rejected automatically during solution. After the process of crossover and mutation, combine the newly generated population and parent population, and calculate fitness value of each chromosome. Use elite strategy to keep the better chromosome to generate a new population.

• Termination Conditions
By setting the maximum number of iteration and test the convergence results of each generation population, we could terminate the heuristic algorithm. In this article, if the number of elements in the Pareto front account for 80 percent or greater, of the population size, the heuristic algorithm terminated. Otherwise, the algorithm could terminate when reaching the maximum number of iterations.

Optimal Trajectory Election Using Multi-Attribute Decision Making (OTE-MADM)
Based on the Pareto optimal fronts obtained by the NSGA-II algorithm, determining a unique solution from the Pareto optimal fronts is a crux in the application of the CDO. Thus, we assess each solution of Pareto front set with entropy-based TOPSIS method [31,32].

Attributes
The economy, efficiency, and safety metrics shall be considered in trajectory planing, which represent common interests among different stakeholders (e.g., pilots, ATC, airlines). Here, although the safety separation between aircraft is guaranteed during trajectory planning, with the consideration of uncertainty factors (e.g., positioning error) in realworld operation, we introduce the collision probability as an additional index to evaluate the potential risks and operational robustness of planned CDO trajectories. It is worth noting that there is no collision in the single-aircraft scenario (refer to Section 8.2), so the collision probability will not be included.
Collision probability is determined by the co-variance matrices of the position errors, shape and track of aircraft [33]. The position vector of aircraft could be obtained by the aircraft performance model, while the position error vector is usually obtained by the prediction error. Assume there are two aircraft α 1 and α 2 in the airspace, the real-time positions of them are denoted as (x α1 , y α1 , h α1 ) and (x α2 , y α2 , h α2 ), respectively. If the probability density function of collision is f ca , the calculation of collision probability can be transformed into an integral of the probability density function of the equivalent region, V area . The formulation is shown as follows: For each pair of aircraft, one is regarded as a reference aircraft, the other is an intruder one. When calculating collision probability, we assume the following conditions for simplicity: • The position error of aircraft satisfies the Gauss distribution. • The co-variance matrices of aircraft positions are not related. For each pair, the position error of the reference aircraft is e r , while the position error of the intruder is e v , and the relative position error of them is e f = e r + e v . Based on position error, the co-variance matrix of the relative random error of two aircraft is expressed as a diagonal matrix: Based on such assumptions, the intruder is regarded as a point and takes the center of the reference aircraft as the origin of the coordinate system, so the collision area of them regarded as a cylinder with height of H and radius of R, where H is the sum of their height, and R is the sum of radius of their envelop, as shown in Figure 5. The relative position of the intruder is (µ x , µ y ) in the horizontal plane, so the probability density function, f ca,XY (x, y), and collision probability, P ca,XY , are: In order to calculate the collision probability in two-dimensional space, a fast computation method based on compressed infinite series is applied to solve the problem of the probability integration. After derivation, the upper form is rewritten as follows [33]: where σ x and σ y are the variance of the probability density function.
Similarly, the relative position of the intruder aircraft is µ z in the vertical plane, so the probability density function f ca,Z (z)and collision probability P ca,Z are: According to the functional nature of normal distribution, we find a rational function, in which the function characteristic is consistent with that of distribution function, and its derivative characteristic is consistent with probability density function of normal distribution. The rational function is: where C ca,Z = 4 √ 2πσ z , and σ z is the variance of the probability density function.

Entropy-Weights Method
In MADM, whether weights of attributes (i.e., fuel consumption, trip time, and collision probability) are reasonable or not is critical to the decision-making accuracy. Entropy Weight Method (EWM) is an important information weight model that has been extensively studied and practiced. Compared with various subjective weighting models, the most significant advantage of the EWM is the avoidance of the interference of human factors on the weight of indicators, thus enhancing the objectivity of the comprehensive evaluation results. Therefore, the EWM has been widely used in multi-criteria decisionmaking [34,35]. Recently, the benefits of the EWM for multi-objective optimization has been also proven [36], although the results are prone to distortion when too many zero values are in the measured indicators [37]. Considering the non-zero features of fuel consumption, trip time, and collision probability, EWM is chosen to calculate the weights of attributes in Pareto's solutions.
Suppose there are m solutions in the Pareto optimal front; each solution is evaluated with respect to n attributes, i.e., fuel consumption, trip time, and collision probability, in which values constitute a decision matrix Z = (z ij ) m×n . Since such attributes in this paper belong to cost value, we normalize each attribute value z ij in decision matrix into a corresponding element r ij with the formulas: r ij = max i (z ij )−z ij max i (z ij )−min i (z ij ) , to facilitate interattribute comparisons. Note that, among classic normalization methods, the minmax standardization method was validated to be a better option [38][39][40]. The normalized decision matrix is denoted as: Therefore, entropy and weight of each attribute are formulated as: where P ij = r ij ∑ m i=1 r ij , if r ij = 0, P ij · ln P ij = 0.

TOPSIS Method
According to the multi-criteria method selection strategy proposed by Wątróbski [31], TOPSIS proposed by Hwang and Yoon, in 1981 [41], was proven to be the most suitable one for this particular decision situation, where the weight of individual criteria has to be considered quantitatively and complete ranking is required. The principle of this method is to rank the alternatives by calculating the distance of each alternative from the positive ideal solution and the negative ideal solution for problems in decision-making, thus determining the optimum alternative. By applying the weights calculated by EWM, entropy-based TOPSIS is adopted to select the optimal CDO trajectories.
As mentioned above, in this paper, we suppose a decision matrix Z = (z ij ) m×n , where there are m solutions with n attributes, i.e., fuel consumption, trip time, and collision probability. The TOPSIS method consists of the following steps: • Normalize the decision matrix. Given that all the attributes of each solution, i.e., fuel consumption, trip time, and collision probability, are cost attributes in the MADM problems, we introduce the following formulas [41] to normalize each attribute value z ij in decision matrix Z = (z ij ) m×n into a corresponding element u ij : • Calculate the weighted normalized decision matrix. Suppose that weights of attributes, W = (w 1 , w 2 , w 3 , · · · , w n ) T , is obtained with Equation (29), where w j ≥ 0, ∑ n j=1 w j = 1, we can construct the weighted normalized decision matrix as • Determine the positive and negative ideal solutions. The Positive Ideal Solution (PIS) A + and the Negative Ideal Solution (NIS) A − are determined, respectively, as follows: where k + j = max 1≤i≤m {k ij }, and y − j = min 1≤i≤m {k ij }. • Measure the distance from each solution to PIS and NIS. The separation, S + i and S + i , are given as: • Calculate the closeness coefficient to the ideal solutions. The closeness coefficient of the i th solution in the Pareto optimal front with respect to the ideal solutions, C i , is defined as: • Rank the preference order. The solutions of the Pareto optimal front can be ranked according to the descending order of C i ; the larger C i means the better solution.

Numerical Tests
In this section, numerical results with different settings of (user-defined) algorithm parameters and various scenarios were presented and discussed to validate the performance of proposed trajectory optimization method in the new airspace. The overall process is run on a 2.00 GHz, CORE i7 CPU, under Windows 10 operating system laptop based on MATLAB code.

Parameter Setting
For simplification, we assume the geopotential pressure altitude is equal to geodetic altitude; true airspeed is equal to ground speed, which implies the international standard atmosphere and no wind. Airbus A320, one of main aircraft types in the aviation industry, is selected as specific aircraft for the CDO trajectory optimization in the experiments. Performance parameters of A320 in BADA version 3.9 are presented in Table 1.

Algorithm Comparison
In order to verify the superiority of NSGA-II in solving proposed trajectory planning problem, two kinds of multi-objective optimization algorithms with good performance and wide application at present [30] have been selected: Speed-constrained Multi-objective Particle Swarm Optimization (SMPSO) and Improving the Strength Pareto Evolutionary Algorithm (SPEA2). Based on the two-dimensional objective function space where the fuel consumption and flight time are taken as the coordinate axes, five evaluation indexes are used to comprehensively demonstrate the performance of algorithms: general distance (GD), inverted generational distance (IGD), hypervolume (HV), Spacing, and maximum Pareto front error (MPFE) [30]. In order to compare the performance characteristics of NSGA-II, SMPSO, and SPEA2 algorithms in different problem scales, three traffic scenarios with different densities, i.e., 1/6/12 aircraft arrive at the entry points simultaneously, are designed. In order to enhance the comparability, the three algorithms have the same convergence conditions. Moreover, the population number of NSGA-II and SPEA2 is 2000, the maximum number of iterations is 80, the crossover probability is 0.3, and the mutation probability is 0.7; in SMPSO, the number of particles and archive size is 2000, and the maximum number of iterations is 80, as well. The parameters corresponding to the optimal solution performance are set as follows: initial velocity influence factor is 0.2, individual optimal influence factor is 0.2, and global optimal influence factor is 0.5. Under the above benchmark parameters, the performance of the three algorithms is shown in Table 2.   2274 1.4045 1.4019 1.5743 1.639 1.6319 1.6432 1.9043 1.8377 The results show that, although the performance of the three algorithms decreases with the increase of the problem size, NSGA-II is always better than the other two algorithms in various performance indicators. In addition, taking the single-aircraft scenario as an example, the convergence rate of NSGA-II is obviously better than the other two algorithms, as shown in Figure 6. Therefore, NSGA-II is capable of converging to the Pareto optimal solution set more efficiently and accurately. Based on NSGA-II algorithm, we will discuss the optimization of continuous descent trajectory in the novel flexible airspace from the perspectives of single aircraft, multi-aircraft simultaneous, and continuous arrival, respectively.

Single Aircraft Scenario
In the single aircraft scenario, the entry point is set to No.5. Some feasible solutions are shown in Figure 7a, where the Pareto-Optimal solutions are in the red diamond. The trajectories of such solutions are shown in Figure 7b. The convergence of Pareto front is shown in Figure 6.
For each Entry Point (EP), we search for the Optimal Trajectory with the least Fuel Consumption (OTFC) and the least Trip Time (OTTT) and get the Ideal Trajectory with Entropy-based TOPSIS method (ITET), respectively, as shown in Table 3, where the "Trajectory" is named by the way-point number of each ring from the second-level to the fifth-level in the ICSAA. Please note that, for simplification, the number of level is omitted.
It is shown that fuel consumption gets higher with the trip time increasing. In view of that, we begin to gain insights into the reasons by investigating the dynamics parameters of aircraft during descent. As shown in Figure 7b, two different trajectories, Trajectory 1 (5-5-5-5) and Trajectory 2 (5-5-5-31), are generated. It is clear that the difference comes from the selection of the way-points at the fifth level, resulting in different flight path angles as mentioned in Section 2. Figure 8a shows the altitude evolves with time. It can be seen that the altitude of Trajectory 2 is always higher than that of Trajectory 1, and so is the velocity (Figure 8b), since we use the standard speed profile recommended by BADA, which shortens the trip time of Trajectory 2. In addition, for Trajectory 2, the time that the aircraft starts executing the CDO approach is around 250 s, as shown in Figure 8a, while the time for Trajectory 1 is about 150 s. Apparently, staying longer in the outside rings need more power according to the operation mode. When aircraft are executing continuous descent, the engine is idle or near-idle and the fuel burn is zero or near zero, as presented in Figure 8c,d. Therefore, more fuel burn of Trajectory 2 is observed.

Multi-Aircraft Scenario
Safe and efficient operation in the high-density terminal area is always one of the great challenges in application of the CDO. In order to further investigate conflict-free trajectory planning based on ICSAA in the high-density situation, two traffic scenarios are designed: simultaneous arrival and continuous arrival, to explore the instantaneous and durative capacity performance of ICSAA.

Simultaneous Arrival
Four traffic scenarios with increasing density of 4, 6, 9, and 12 are designed to explore the transient operation performance in the ICSAA. For each scenario, the entry points are evenly spaced along the outermost ring. By estimating the time and speed of aircraft flying over the entry point in advance, optimal trajectories can be generated within 2 min. For each scenario, the Pareto optimal fronts are obtained with NSGA-II, and the collision probability of each solution is computed, as shown in Figure 9, where the solutiones of ITET, marked with red diamond, are obtained by OTE-MADM. In general, as the number of aircraft increases, average fuel consumption, average trip time and collision probability increase gradually. The results also suggest that, with the deepening of airspace congestion, the trade-off relationship is getting convergent. Taking the scenario with 4 aircraft as an example, the trajectories corresponding to OTFC, OTTT, and ITET are shown in Figure 10. The main difference is still the choice of way-points at the fifth level of ICSAA, which also means the different flight path angles of trajectories. In the multi-aircraft scenario, choosing different flight path angle is helpful to resolve conflicts between aircraft, while having significant impact on the fuel consumption and trip time, as presented in Table 4.
In order to further evaluate the optimization effect, we set ITET in the single aircraft scenario based on the Table 3 as reference. It means reference total fuel consumption is the sum of fuel consumption of each ITET and so is the reference total trip time. Then, we compare the reference value with the corresponding value of OTFC, OTTT, and ITET, as shown in Table 4, and the increase rate is shown in Figure 11. The results demonstrate that, as the number of aircraft increases, the average fuel consumption, average trip time and collision risk show an increasing trend. Under different congestion, when the ITET is picked, the increase rate of average fuel consumption is less than 20%, while the increase rate of average trip time is less than 10%. It is worth to note that there are some negative increases of fuel consumption of OTTT in some scenarios. The reason is that the ITET in the single aircraft scenario is selected as the reference, and, in some cases, the ITET in single aircraft scenario actually is the trajectory with the minimum trip time, so the fuel consumption is relatively large. But, in the multi-aircraft scenario, it is generally not possible to follow the trajectory with the maximal fuel consumption of the single aircraft scenario, since they have to take some maneuver to avoid conflicts. And such a maneuver does not take too much cost (i.e., fuel consumption and trip time) because the density of them is relatively low. Therefore, there is a path solution space with lower fuel consumption, and their fuel consumption could be less than the reference value.  In summary, when conducting CDO based on ICSAA in the simultaneous arrival, the optimization model we proposed could generate and get a set of optimized conflict-free four-dimensional trajectories for each aircraft. Based on a set of optimization results, we get a unique ideal solution by OTE-MADM, which makes accurate trade-offs among fuel consumption, trip time, and collision risk. In addition, compared the reference value, the increase of fuel consumption and trip time is optimized to an acceptable level.

Continuous Approach
In order to explore the continuous performance of CDO based on ICSAA, this paper has designed 4 traffic scenarios that four entry points are distributed uniformly in the outermost ring, and the number of aircraft increases from 4 to 16 by the step of 4, and the arrival time interval at the same entry point follows the Poisson distribution with the expectation of 2 min. In addition, the minimum separation at the entry points should be also guaranteed for safety.
We firstly generate planned trajectory for each aircraft based on ITET trajectory in Table 3. In order to further investigate the optimization ability of proposed models and algorithms, a Monte Carlo simulation was adopted to pick out the worst traffic pattern that contains the highest collision risk for each scenario. Flight trajectories in the ICSAA are then efficiently and optimally selected by proposed methods.
The Pareto optimal fronts and collision probability of each solution were obtained within 3 min, as presented in Figure 12, where the solution of ITET is marked with the red diamond. In general, with the increase of traffic density, the trade-off decision space between fuel consumption and trip time compresses. In the low-density scenario, the potential conflict between aircraft is less, so they have more choices, while, in the high-density scenario, the coupling degree among aircraft is high, and the airspace situation is complex, the solution space of conflict-free trajectories is reduced. For such optimization results, ITET is obtained with OTE-MADM, as presented in Table 5. The results demonstrate that the average fuel consumption and the average trip time present an increasing trend with the number of aircraft rises. Similar to the previous subsection, compared with the reference value (i.e., reference total fuel consumption, reference total trip time, etc.), as shown in the Figure 13, all the increase rate of the fuel consumption and the trip time is less than 30% and 15%, respectively, with the number of aircraft increasing.   In conclusion, in the continuous approach scenario, the proposed algorithm is still efficient and reliable. Although the additional ratio of average fuel consumption and trip time after optimization continue to increase with the increase of traffic density, it is still at an acceptable level. OTE-MADM could also make a good trade-off among fuel consumption, trip time, and collision probability. It can effectively support the highperformance operation requirements of high-density terminal area and help to promote the level of coordination between ATC and airlines under TBO mode.

Conclusions and Future Work
Improving flexibility of approach airspace is one of cornerstone for implementing CDO especially in high-density traffic scenarios. Inspired by traditional PMS, this paper broke through the traditional rigid terminal airspace structure by designing ICSAA and its operational procedures to deliver Omni-directional CDO. A multi-objective trajectory optimization model aiming at minimizing fuel consumption and trip time was solved by NSGA-II. In order to provide decision-making support to air traffic controller, the best trade-off was determined using entropy-based TOPSIS from the perspective of safety, efficiency and economy.
The proposed models and algorithms were validated in the single aircraft scenario, the multi-aircraft scenarios with simultaneous and continuous arrivals to gain insights into the interplay of traffic demand and given objectives. The results verified the effectiveness and efficiency of optimal CDO operation in ICSAA dealing with complex traffic patterns, which demonstrates the potential of CDO operation in high-density terminal area and enables the application of TBO. This paper is a preliminary exploration of CDO trajectory planning in the novel airspace. Further research shall be carried out to solve the limitations and to facilitate the transition from theoretical studies to real-world practice.

•
The ICSAA was designed based on an artificial airspace irrespective to terrain, obstacles, noise abatement, airspace restriction and runway configurations, etc. The very essential work is to adapt the ICSAA to a real airport and thoroughly validate the effectiveness by both fast-time simulations and Human-In-The-Loop (HITL) experiments. • Aircraft positioning error was considered in trajectory selection in this paper. However, more realistic uncertain factors, like flight control error, meteorological conditions, stochastic runway scheduling, etc., shall also be included in robust CDO trajectory planning and implementation. • Human-machine system is another emerging topic in future autonomous ATM. The decision-making by machine shall be transparent and understandable, which allows air traffic controllers and pilots to take over automated equipment at any time safely. It means the cognition of human and machine needs to be strongly synchronized [42]. For automated CDO trajectory planning and negotiation in flexible airspace, designing transparent Decision Support System (DSS) based on explainable Artificial Intelligence(AI) would be an indispensable approach to enhance human cognition, operational efficiency, and safety, especially for a high-density traffic scenario. Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
The following abbreviations are used in this manuscript: