# Complexity Optimal Air Traffic Assignment in Multi-layer Transport Network for Urban Air Mobility Operations

Daniel Delahaye, Sameer Alam, Jean-Loup Farges, Zhengyi Wang

► To cite this version:

Daniel Delahaye, Sameer Alam, Jean-Loup Farges, Zhengyi Wang. Complexity Optimal Air Traffic Assignment in Multi-layer Transport Network for Urban Air Mobility Operations. Transportation research. Part C, Emerging technologies, 2022, 142, pp.103776. 10.1016/j.trc.2022.103776. hal-03721341

HAL Id: hal-03721341

https://enac.hal.science/hal-03721341v1

Submitted on 1 Aug 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Complexity optimal air traffic assignment in multi-layer transport network for Urban Air Mobility operations

Zhengyi Wang<sup>a,*</sup>, Daniel Delahaye<sup>a</sup>, Jean-Loup Farges<sup>b</sup>, Sameer Alam<sup>c</sup>

## ARTICLE INFO

**Keywords:**

- Urban Air Mobility
- Unmanned traffic management
- Air Traffic Management
- Air traffic assignment
- Air traffic complexity
- Linear Dynamical System

## ABSTRACT

Large numbers of Urban Air Mobility (UAM) vehicles are expected to operate in urban airspace in the near future, exceeding the capacities of current airspace and Air Traffic Management (ATM) systems. This paper presents an air traffic assignment framework for 3D air transport networks in urban airspace to enable UAM operations at future demand levels. The individual vehicle dynamics are aggregated to describe the overall flow feature in this macroscopic model. Firstly, UAM operations are modeled as flows and structured in a three-dimensional two-way air transport network. Then, a complexity optimal air traffic assignment in urban airspace is formulated as an optimization problem. Based upon the Linear Dynamical System (LDS), a novel complexity metric is defined as objective function, which takes into account dynamic flow structure, congestion, and operational efficiency. A two-phase approach combining Simulated Annealing (SA) and Dafermos' Algorithm (DA) is introduced to efficiently solve this problem. To validate the proposed model, a case study of an air transport network in Singapore's urban airspace with two different demands is conducted. Comparative studies are carried out between the proposed algorithm and other widely used traffic assignment algorithms. The results show that the proposed approach is capable of assigning flows in an efficient and effective manner, reducing the complexity of the air transport network significantly. The results also show that optimizing the flow pattern reduces total complexity by 90.44% ± 0.53% and 92.12% ± 0.35% with 95% confidence interval, respectively in two scenarios. The framework may be useful for Air Navigation Service Providers (ANSP) in strategic planning for UAM operations and urban airspace design.

## 1. Introduction

The urbanization process gradually shifts the population from rural to urban areas. The United Nations (UN) estimates that by 2050, 68.4% of the world population is expected to be urban, and the urban population will increase by 2.5 billion (UN, 2018). This growth will put strain on transportation infrastructure, resulting in undesirable outcomes such as traffic congestion, complex traffic situations, delays, etc. (Bauranov and Rakas, 2021). As a result, innovative mobility solutions are required to ensure safe and efficient urban transportation. Urban Air Mobility (UAM) is one of the most effective concepts for achieving high levels of automation in metropolitan areas. UAM provides on-demand or scheduled operations using manned and unmanned aircraft systems

E-mail addresses: zhengyi.wang@alumni.enac.fr (Z. Wang), delahaye@recherche.enac.fr (D. Delahaye), jean-loup.farges@onera.fr (J.-L. Farges), sameeralam@ntu.edu.sg (S. Alam).

<footnote><sup>a</sup> OPTIM Group, ENAC Lab, Ecole National de l'Aviation Civile, Toulouse, France</footnote>

<footnote><sup>b</sup> ONERA, Toulouse, France</footnote>

<footnote><sup>c</sup> ATMRI, Nanyang Technological University, Singapore</footnote>

<footnote>* Corresponding author.</footnote>
for intra- or inter-urban transport (Thipphavong et al., 2018). In particular, Unmanned Aerial Vehicles (UAVs) play an important role in UAM operations in low-level urban airspace. Services such as passenger transportation, parcel delivery, and traffic monitoring are expected to be provided by UAVs.

To make the best use of urban airspace, some preliminary projects of urban airspace design have been conducted (Singireddy and Daim, 2018; FAA, 2020; Schneider et al., 2015; Jang et al., 2017; Lascara et al., 2019; Geister, 2017; Sunil et al., 2015). These projects envisioned UAM in the low-altitude urban area in class G airspace, which is uncontrolled airspace below 1200 ft and not equipped with any ATM services (Bauranov and Rakas, 2021). It is further classified into different categories and layers. Aircraft operate within a fixed volume and must adhere to specific procedures, rules, and performance requirements. These volume segments create a multi-level air transport network with high throughput, efficiency, and safety. To handle excessive future demands, additional capacity can be easily enhanced by expanding the internal structure of the volume segment. The separation of aerial vehicles within a volume segment can be assured in the tactical planning phase by Conflict Detection and Resolution (CD&R).

However, a major challenge is integrating UAM with the Air Traffic Management (ATM) system for ever-increasing and high-density UAM operations in class G airspace. To this end, a variety of research related to air traffic assignments are underway in order to decongest the air transport network under current airspace capacity (Deschinkel et al., 2002; Delahaye et al., 2005; Nosedal et al., 2014; Ho-Huu et al., 2019). Air traffic assignment remedies the deficiencies in the existing transportation system by assigning estimated future demands on the existing systems. The flow structure is expected to be strategically optimized to ensure low-complexity operation in the specified urban airspace. In addition, the remaining conflicts will be substantially reduced and much easier to solve in the tactical planning phase (FAA, 2020).

In order to assess the inter-dependency of congestion and the degree of difficulty in controlling UAM operations, it is essential to define an intrinsic air traffic complexity measure that depends on the geometry of airspace and is independent of the traffic control system (Treimuth et al., 2015). In the literature, several metrics have been developed and applied to measure air traffic complexity, including convergence metric (Delahaye and Puechmorel, 2000), geometric metric (Delahaye and Puechmorel, 2000), proximity metric (Delahaye and Puechmorel, 2000), clusters metric (Histon et al., 2002), Grassmannian metric (Delahaye and Puechmorel, 2013), König metric (Essén, 1993; Delahaye et al., 2002; Juntama et al., 2020). Although most of these metrics are proven to be effective on basic examples, they are not well suited for large-scale or complex applications. Additionally, spatial-temporal information is not taken into account. To overcome these limitations, an air traffic complexity metric based on the Linear Dynamical System (LDS) has been introduced (Delahaye and Puechmorel, 2000). LDS has been successfully applied in many studies (Delahaye and Puechmorel, 2000; Delahaye et al., 2002, 2004; Delahaye and Puechmorel, 2010; Treimuth et al., 2015; García et al., 2020; Wang et al., 2021). With a concise mathematical formulation, it is effective as an intrinsic air traffic complexity measurement. The property of LDS allows us to measure the complexity of air traffic by considering the disorder, interaction, and evolution of aircraft.

To alleviate the congestion and organize the air traffic flow structure in urban airspace, a preliminary air traffic assignment model for a planar urban air transport network has been conducted by authors (Wang et al., 2021). Based on previous researches and the aforementioned discussion, we develop a complexity optimal air traffic assignment model to organize intensive UAM traffic in a multi-layer two-way air transport network. The model is macroscopic that aggregates individual vehicle dynamics to describe the overall stream features in the network, such as dynamic flow structure, congestion, and efficiency. It is one of the first attempts to conduct air traffic assignments for UAM operations from a macroscopic perspective. The contributions of this paper are as follows:

* A novel macroscopic problem formulation for air traffic assignment for UAM operations.

* An efficient representation of air transport network for UAM.

* An intrinsic air traffic complexity metric based on the LDS.

* A two-phase optimization procedure on the basis of SA and DA.

* Application and comparative study on an air transport network in Singapore's urban airspace.

## 2. Air transport network design for UAM operations

In this section, we introduce the design and representation of the air transport network in UAM operations. For ease of illustration, Table 1 lists the main notations of sets and variables used throughout this paper.

### 2.1. Air transport network modeling

An air transport network $\mathcal{G}$ consists of nodes $\mathcal{V}$ and links $\mathcal{E}$. The nodes represent the waypoints, droneports, and delivery stations. The position of a node $\mathbf{v}$ is given as $(x_v, y_v, z_v)$. A link $e \in \mathcal{E}$ represents the airway that joins a pair of nodes. All links are considered to be bidirectional. Each link $e$ is characterized by its length $L_e$ and the nominal speed $s_e$ of the flow traveling on it. For simplicity, we assume that all UAM vehicles travel at free-flow speed, and all upward and downward links are vertical.

Volume segments are designed to make up the links in the air transport network (FAA, 2020; Geister, 2017). It supports point-to-point UAM operations and efficient networks with routing. Fig. 1(a) illustrates an example of volume segments for two horizontal links in opposite directions. The vertical view is shown in Fig. 1(b). The bi-directional links are made up of several volume segments with length $L_e$ and width $L_w$, located on opposite sides. The number of lanes on each side is $N_s$. Each volume segment is divided into 3D grid segments with the length $L_g$. All grid segments have a length of $L_l$, except that the grid length in the middle is in $[L_l, 2L_l[$. This is because the length varies for each link and the grid cannot be divided equally. The height of the grid segment, length $L_l$, and width $L_w$ should be determined according to the size of the protected zone of aerial vehicles.
<table><caption>**Table 1**
Notations and explanations.</caption><thead><tr><th>Notations</th><th>Explanations</th></tr></thead><tbody><tr><td><strong>Sets</strong></td><td></td></tr><tr><td>G</td><td>Air transport network</td></tr><tr><td>ε</td><td>Set of links</td></tr><tr><td>γ</td><td>Set of nodes</td></tr><tr><td>W</td><td>Set of OD pairs</td></tr><tr><td>P<sub>w</sub></td><td>Set of feasible paths between OD pair w</td></tr><tr><td>A<sub>v</sub></td><td>Cylindrical airspace around node v</td></tr><tr><td>ρ<sup>i</sup><sub>v</sub>, ρ<sup>d</sup><sub>v</sub></td><td>Convergent, divergent density points in A<sub>v</sub> of node v</td></tr><tr><td>ρ<sub>e</sub></td><td>Set of density points of link e</td></tr><tr><td>γ</td><td>Set of path flow proportion</td></tr><tr><td>f</td><td>Set of link flow</td></tr><tr><td><strong>Variables</strong></td><td></td></tr><tr><td>F<sup>i</sup><sub>w</sub></td><td>Flow on i-th path of OD pair w</td></tr><tr><td>F<sup>d</sup><sub>v</sub></td><td>Flow on path p</td></tr><tr><td>f<sub>e</sub></td><td>Flow on link e</td></tr><tr><td>d<sub>w</sub></td><td>Demand between OD pair w</td></tr><tr><td>w<sub>g</sub></td><td>Density on density point g</td></tr><tr><td>γ'<sub>w</sub></td><td>The proportion of flow assigned to i-th path for OD pair w</td></tr><tr><td>R<sub>w</sub></td><td>The number of alternative paths |P<sub>w</sub>| between OD pair w</td></tr><tr><td>δ<sub>e,w,i</sub></td><td>= 1 if link e is contained in the i-th path of OD pair w, otherwise 0</td></tr><tr><td>δ<sub>e,p</sub></td><td>= 1 if link e is contained in path p, otherwise 0</td></tr></tbody></table>

<figure><img src="image_2.png" /><figcaption>**Fig. 1.** An example of volume segments of horizontal links with opposite directions.</figcaption></figure>

Density points of the set $\rho_e$ are introduced as virtual representative points for each volume segment. As a representation of the amount of flow in each grid segment, the density point is located in the center of each grid, which is shown in Fig. 1(b). Each density point $g \in \rho_e$ is characterized by position $(x_g, y_g, z_g)$, heading $\theta_g$, speed $s_g$, and density $w_g$. The position is determined based on the position of endpoints of the link. The heading and speed information is derived from the flow traveling on the corresponding volume segments. The average amount of air traffic flow distributed on each grid segment $g$ is referred to as the density $w_g$, which is defined as follows:

$$w_g = \frac{f_e L_g}{N_s s_e} \qquad (1)$$

where $f_e$ is the link flow.

## 2.2. Feasible paths generation

A path consists of a finite sequence of links that connects an OD pair. For each OD pair $w \in W$, the feasible paths $P_w$ in a bidirectional network are generated without cycles based on depth-first search (Migliore et al., 1990). Furthermore, the range and endurance of any given flight are regarded as safety concerns. For any path $p$, its length $L_p$ must be within the range of UAV $R_{max}$:

$$L_p = \sum_{e \in p} L_e < R_{max}, \quad p \in P_w, w \in W \qquad (2)$$
The flight time should not exceed the endurance $T_{max}$ of the UAV:

$$ \sum_{e \in \mathcal{P}} \frac{L_e}{s_e} < T_{max}, \quad p \in \mathcal{P}_w, w \in \mathcal{W} \qquad (3) $$

Maintenance and battery recharging (or refueling) are assumed to be available not only at droneports but also at delivery stations. Note that, if not available, the values of $R_{max}$ and $T_{max}$ should also consider a return trip.

### 2.3. Demand modeling

Delivery services are launched from droneport to delivery stations, and the drone will also return to droneport for recharging and maintenance purposes. Based on this assumption, both a droneport and delivery station can be considered as origin or destination. Each pair of Origin-Destination (OD) $w \in \mathcal{W}$ is given with an associated demand $d_w$. This traffic demand will be distributed among paths in $\mathcal{P}_w$. Note that, the demands between an OD pair $w = (v_i, v_j)$ and the inverted OD pair $w' = (v_j, v_i)$ are frequently imbalanced.

## 3. Air traffic assignment problem formulation

Given a complete description of the proposed or existing air transport network, air traffic assignment aims to allocate the demand of a given set of OD pairs that minimizes the complexity of the air transportation system while keeping flight efficiency and operational efficiency within acceptable limits. Two main approaches are possible for modeling flow assignment problems. In the first approach, the variables are the number of flows of UAM vehicles traveling on the links for each OD pair. The constraints at each node enforce the traffic flow conservation. In the second approach, the variables are defined as the flow distributed on a set of paths that connect each OD pair. Here, the second approach is chosen. Thus, the proposed method relies on an intelligent choice of the set of paths $\mathcal{P}_w$ for each OD pair $w$.

### 3.1. Decision variables

The decision variables are the proportion of flow $\gamma = (\gamma_w)_{w \in \mathcal{W}}$ assigned to paths for each OD pair $w \in \mathcal{W}$, where $\gamma_w = (\gamma_w^1, \dots, \gamma_w^{R_w})$.

### 3.2. Constraints

The decision variables are non-negative and satisfy the conservation law:

$$ \gamma_w^i \ge 0, \quad i \in \{1, \dots, R_w\}, \quad w \in \mathcal{W} \qquad (4) $$

$$ \sum_{i=1}^{R_w} \gamma_w^i = 1, \quad w \in \mathcal{W} \qquad (5) $$

The flow of ith path $F_w^i$ is derived from the demand and the proportion of path flow:

$$ F_w^i = d_w \gamma_w^i, \quad i = 1, \dots, R_w, w \in \mathcal{W}, \qquad (6) $$

In addition, the air traffic flow should be meaningful for assignment. If the flow on a path is positive, it should be greater than a path flow threshold $F_e$:

$$ F_w^i \in \{0\} \cup [F_e, d_w], \quad i = 1, \dots, R_w, w \in \mathcal{W} \qquad (7) $$

Then, link flow $f_e$ can be calculated from path flow:

$$ f_e = \sum_{w \in \mathcal{W}} \sum_{i=1}^{R_w} \delta_{e,w,i} F_w^i, \quad e \in \mathcal{E}, i = 1, \dots, R_w, w \in \mathcal{W}, \qquad (8) $$

To model the operational efficiency, it is assumed that air traffic flows on all links have the same characteristics. In a multi-level network, there are three flight phases: climb, descend, and cruise. Therefore, an energy consumption coefficient $\alpha_e$ is associated with each link according to the flight phase. Then, $\alpha_e$ is introduced to measure the operational efficiency on link $e$ in terms of distance and flight phases:

$$ \alpha_e = \begin{cases} L_e \varphi_h, & e \in \mathcal{E}_h \\ L_e \varphi_u, & e \in \mathcal{E}_u \\ L_e \varphi_d, & e \in \mathcal{E}_d \end{cases} \qquad (9) $$

where $\varphi_h$, $\varphi_u$, and $\varphi_d$ are efficiency coefficients per distance unit depending on whether the link is horizontal, upward, or downward.

Based on $\alpha_e$, the energy consumption $E$ for flights through air transport network is defined as follows. It measures link congestion by involving link flow $f_e$.

$$ E = \sum_{e \in \mathcal{E}} \alpha_e f_e \qquad (10) $$
An efficiency threshold $\beta \hat{E}$ is set to ensure the network efficiency remains within an acceptable range:

$$E \le \beta \hat{E} \qquad (11)$$

where $\beta$ is a factor to describe the proportional increase and $\beta > 1$. $\hat{E}$ is the minimum network efficiency, that is, the minimum value of Eq. (10) when the demand $d_w$ is completely assigned to the path with the lowest efficiency for all OD pairs $w \in W$:

$$\hat{E} = \sum_{w \in W} d_w \min_{\rho \in \mathcal{P}_w} \left( \sum_{e \in \rho} \alpha_e \right) \qquad (12)$$

### 3.3. Objective function based on linear dynamical systems

The objective function aims to minimize the complexity of flow patterns in a UAM air transport network. Because volume segments organize the traffic in parallel, the complexity arises mainly at nodes. Therefore, the objective function is defined as the sum of node complexities:

$$\mathcal{C}(f) = \sum_{v \in V} C_v(f) \qquad (13)$$

where $f$ is the complete flow pattern $\{f_e|e \in E\}$. $C_v$ depends on $f$ if the network configuration is fixed.

The procedure for computing the node complexity $c_v$ is based on LDS, which can be divided into three steps:

#### (1) Measure the node complexity by traffics within the surrounding airspace

To capture the traffic situation around each node $v \in V$, we define a cylindrical airspace $\mathcal{A}_v$ with centroid $v$, radius $r_c$, and height $h_c$. The density points within the cylinder are used to measure the complexity. Each cylinder $\mathcal{A}_v$ includes a set of density points. It is worth noting that the traffic flow in each volume segment can be organized into segregated forward or inverted directions in terms of nodes. When all traffic within a cylinder is considered together, there is a risk of canceling each other out. Thus, depending on the mode of operation, they are divided into convergent density points $\rho_v^c$ which are inwards towards the node $v$, and divergent density points $\rho_v^d$ which are outward from the node $v$, where $\rho_v^c \cap \rho_e \neq \emptyset$ and $\rho_v^d \cap \rho_e \neq \emptyset$ if node $v$ is connected by link $e$. It should also be noted that $\rho_v^c \notin \rho_e$ and $\rho_v^d \notin \rho_e$ since a cylinder $\mathcal{A}_v$ may also include density points $\rho_e'$ located on other links $e' \neq e$. Fig. 2 illustrates an example of the cylindrical airspace around a node connected by six links. Except for density points outside the cylinder (in blue), convergent (in green), and divergent (in red) density points inside the cylinder are used separately to measure the air traffic complexity.

#### (2) Estimate the parameters of LDS and calculate the eigenvalues of the associated coefficient matrix

LDS is proposed to capture the intrinsic air traffic complexity by describing the evolution of a given traffic situation. LDS is defined by the following equation:

$$\dot{\mathbf{x}} = \mathbf{A}\mathbf{x} + \mathbf{b} \qquad (14)$$

where $\mathbf{x}$ denotes the state vector, $\mathbf{x}$ represents the speed vector, $\mathbf{A}$ stands for a coefficient matrix, $\mathbf{b}$ is a coefficient vector. In the conventional representation of LDS, the state vector is generated by directly concatenating the positions of all aircraft in a given airspace (Delahaye and Puechmorel, 2000). To improve the computation efficiency, we propose a Weighted Minimum Mean Square Error (WMMSE) estimator to derive the coefficient matrix of a 3-dimensional LDS. WMMSE solves the following optimization problem:

$$\min_{\mathbf{A},\mathbf{b}} \sum_{g \in \mathcal{G}} \omega_g \left\| \dot{\mathbf{x}}_g - (\mathbf{A}\mathbf{x}_g + \mathbf{b}) \right\|^2 \qquad (15)$$

where $\rho$ are divergent density points $\rho_v^d$ or convergent density points $\rho_v^c$. $w_g$ is the density at point $g$ and $\mathbf{x}_g = (x_g, y_g, z_g)^\mathrm{T}$ is the three-dimensional position vector of point $g$. The speed vector of point $g$ depends on the phase of associated link:

$$\dot{\mathbf{x}}_g = \begin{cases} s_g(\cos(\theta_g), \sin(\theta_g), 0)^\mathrm{T}, & g \text{ attached to } e \in \mathcal{E}_h \\ (0, 0, s_g)^\mathrm{T}, & g \text{ attached to } e \in \mathcal{E}_u \\ (0, 0, -s_g)^\mathrm{T}, & g \text{ attached to } e \in \mathcal{E}_d \end{cases} \qquad (16)$$

where $\theta_g$ is the heading at point $g$ and only exists in cruise phase, $s_g$ is the speed of flow on the associated link. The details for calculating the coefficient matrix $\mathbf{A}$ that minimizes the WMMSE is given in Appendix A.

**Proposition 2** (see Appendix B) enables us to partially characterize the qualitative behavior of the system by eigenvalues of coefficient matrix $\mathbf{A}$. The proof is given as well in Appendix B. Fig. 3 shows the evolution of LDS in terms of eigenvalues in the complex plane. The organized traffic situation corresponds to the set of eigenvalues within the blue strip. For more details, four typical traffic situations and their eigenvalue loci are summarized in Appendix D.

Based on the aforementioned discussion and Proposition 3, the dominant eigenvalue is then introduced as a representative of all eigenvalues to describe the stability of the whole LDS:

**Definition 1.** Given $\lambda_1, \dots, \lambda_N$ as the eigenvalues of a matrix, if $\hat{\lambda} \in \{\lambda_1, \dots, \lambda_N\}$, and $\mathrm{Re}(\hat{\lambda}) = \max(\mathrm{Re}(\lambda_1), \dots, \mathrm{Re}(\lambda_N))$, then $\hat{\lambda}$ is the dominant eigenvalue.
<figure><img src="image_3.png" /><figcaption>Fig. 2. A example of cylindrical airspace, density points and eigenvalue calculation of the associated LDS. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

The proof of Proposition 3 is given as well in Appendix C.

To summarize, Fig. 2 presents the framework for calculating the dominant eigenvalues $\lambda_v^c$ and $\lambda_v^d$ of the coefficient matrix in terms of different flow modes.

### (3) Formulate the cost function for each node

To model the link congestion and operational efficiency in the cost function of node $v \in V$, the efficiency factor $\alpha_v$ and flow $f_v$ of associated link $e \in \mathcal{E}_v^m$ is also taken into account, where $m \in \{c, d\}$ is the operational mode, $\mathcal{E}_v^c$ denotes the set of links with head $v$, and $\mathcal{E}_v^d$ denotes the set of links with tail $v$. The node complexity is defined as follows:

$$C_v(f) = \sum_{m \in \{c,d\}} |\text{Re}(\lambda_v^m)| \sum_{e \in \mathcal{E}_v^m} \alpha_e f_e \quad (17)$$

Finally, the optimization problem can be formulated as follows:

$$\min_{f} \mathcal{C}(f) = \sum_{v \in V} C_v(f)$$

s.t. $E \le \beta \hat{E}$

$$\sum_{i=1}^{R_w} \gamma_w^i = 1, \quad w \in W \\ F_w^i = d_w \gamma_w^i, \quad w \in W, i = 1, \dots, R_w \quad (P1)$$

$$\sum_{w \in W} \sum_{i=1}^{R_w} \delta_{e,w,i} F_w^i = f_e, \quad e \in \mathcal{E}, w \in W, i = 1, \dots, R_w \\ F_w^i \in \{0\} \cup [F_t, d_w], \quad w \in W, i = 1, \dots, R_w \\ \gamma_w^i \ge 0, \quad i = 1, \dots, R_w, \quad w \in W$$

Obviously, the problem (P1) involves large numbers of continuous variables, alternative paths, and air traffic demands. Moreover, the objective function is non-separable on each link, because the node complexity depends on the flow of connected links. It can be proved that problem (P1) is NP-hard (Jahn et al., 2000). We then seek efficient heuristic algorithms to produce good-quality solutions in reasonable computation time for this complex, high-dimensional, and NP-hard problem.

## 4. A two-phase algorithm for air traffic assignment problem

### 4.1. Algorithm overview

We propose a two-phase optimization algorithm for solving the problem (P1). Simulated Annealing (SA) is used in the first phase to explore a good-quality solution as initial values for the next phase. Using Dafermos' algorithm (DA), the second phase aims to further optimize the problem by exploiting the search space based on the initial solution generated by SA.
<figure><img src="image_4.png" /><figcaption>Fig. 3. Evolution in terms of eigenvalues of coefficient matrix A. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

The resulting optimization architecture is shown in Fig. 4. It is made up of three main modules: initialization, simulation, and optimization module. The simulation module is firstly built by constructing the air transport network in UAM architecture and graph representation. The initialization model consists of an efficient SA algorithm. It aims to provide a good quality initial solution for the optimization module. In the optimization module, DA directly controls the decision variables $\gamma$. Based on $\gamma$, the simulation module generates the associated flow pattern in the air transport network for computing the objective function $C$, which leads the optimization module to search for solutions that are operational in UAM operations with the minimum complexity. In addition, performance indicators are introduced to control the selection of path to transfer flow for each path in the optimization module to adjust the decision variables $\gamma$.

## 4.2. Simulated annealing

SA is a powerful metaheuristic algorithm. The advantage of SA has been thoroughly discussed in Delahaye et al. (2019). Due to the high-dimensional continuous search space in the problem (P1), it is not appropriate to use SA directly to solve this problem. Thus, the SA algorithm is used to initialize the DA.

Fig. 5 illustrate the mechanism for each transition of SA. Based on the current state $\gamma_i$, the neighbor state $\gamma_j$ is generated. A proper neighborhood selection strategy can significantly speed up the convergence of the SA algorithm (Ming, 2008). In this study, we propose a probabilistic performance-based neighborhood generation strategy. The goal is to prevent DA from falling into local optimum by initializing with the flow result of SA. The performance indicator aims to improve the convergence performance of SA, and the probabilistic selection of paths and flow is designed to improve the assignment of flow patterns. The pseudocode of neighborhood generation is presented in Algorithm 1. The path complexity is selected as the performance indicator. It is defined as the sum of complexity of the nodes on a path:

$$C_p(f) = \sum_{v \in p} C_v(f), \quad p \in P_w, \quad w \in W \qquad (18)$$

We aim to generate an initial flow pattern such that paths in each OD pair have similar complexities. To this end, for each pair of OD $w$, a random flow $\delta$ is transferred from a path of higher complexity to a path of lower complexity. To limit the variation between the current state and generated neighbor, $\delta$ is set to be a small value. The selection probability of the outflow path in each OD pair $w \in W$ is given by:

$$\mathcal{P}_{w,i}^{\text{out}} = \begin{cases} \frac{C_p(f) - \min_{p \in P_w, F_p > 0} C_p(f)}{\sum_{p \in P_w, F_p > 0} (C_p(f) - \min_{p \in P_w, F_p > 0} C_p(f))}, & F_p > 0 \\ 0, & \text{Otherwise} \end{cases} \qquad (19)$$
<figure><img src="image_5.png" /><figcaption>Fig. 4. The optimization framework of two-phase proposed algorithm on a simulation environment.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>Fig. 5. Probabilistic performance-based neighborhood generation and evaluation of objective function for each transition in SA.</figcaption></figure>

and the selection probability of the inflow path is given as:

$$ \mathcal{F}_{w,i}^{\text{in}} = \frac{1/C_p(f) - \min_{p \in P_w} (1/C_p(f))}{\sum_{p \in P_w} (1/C_p(f) - \min_{p \in P_w} (1/C_p(f)))} \quad (20) $$

Metropolis principle (Metropolis et al., 1953) is used to determine whether a neighbor state can be accepted. If yes, the state is updated by its neighbor in terms of current OD pair. If the neighbor is rejected, only the variables of current OD pair need to be recovered. When the stopping criterion of SA is met, the final state will be used as initial flow pattern for DA in the optimization module.

### 4.3. Dafermos’ algorithm

As an efficient sequential decomposition algorithm, DA has been used to efficiently solve traffic assignment problems (Dafermos and Sparrow, 1969; Dafermos, 1971, 1972). Though DA has strong theoretic support and fast convergence speed, it is computationally expensive to determine the amount of flow and the paths for flow transferring due to the marginal cost computation. In
**Algorithm 1** Probabilistic performance-based neighborhood generation strategy

**Input:**

$\gamma$: Current state

**Output:**

$\hat{\gamma}$: Generated neighborhood of $\gamma$

1: **procedure** GENERATENEIGHBOR($\gamma$)

2: for $w \in \mathcal{W}$ do

3: do

4: Calculate $\mathcal{P}_{w,i}^{\text{out}}$, $\mathcal{P}_{w,i}^{\text{in}}$, $i = 1, \dots, R_w$

5: $r_1, r_2 \leftarrow \text{UNIFORM}(0, 1)$

6: $i_1, i_2 \leftarrow 0$

7: do

8: $r_1 \leftarrow r_1 - \mathcal{P}_{w,i_1}^{\text{out}}$

9: $i_1 \leftarrow i_1 + 1$

10: while ($r_1 > 0$)

11: do

12: $r_2 \leftarrow r_2 - \mathcal{P}_{w,i_2}^{\text{in}}$

13: $i_2 \leftarrow i_2 + 1$

14: while ($r_2 > 0$)

15: $\delta \leftarrow \min(\text{UNIFORM}(0, 5F_i/d_w), \gamma_w^{i_1})$

16: $\hat{\gamma}_w^i \leftarrow \gamma_w^i$, $i \in \{1, \dots, R_w\} \setminus \{i_1, i_2\}$

17: $\hat{\gamma}_w^{i_1} \leftarrow \gamma_w^{i_1} - \delta$

18: $\hat{\gamma}_w^{i_2} \leftarrow \gamma_w^{i_2} + \delta$

19: while (constraint in problem (P1) is not satisfied)

20: end for

21: return $\hat{\gamma}$

22: end procedure

addition, the convergence speed of DA highly depends on the choice of the initial value, especially when there are a lot of route choices or the cost function is not strictly convex. That is why DA will perform better if it is initialized by SA. DA requires successive application of the equilibrium operator to decompose the optimization problem into a sequence of simpler ones. The resulting flows satisfy the system equilibrium property.

DA enumerates all OD pairs. For each OD pair $w \in \mathcal{W}$, the most expensive route $p_{w,max}$ and the least expensive route $p_{w,min}$ in terms of the performance indicators are identified. Marginal cost $C'_p(f)$ is introduced as the performance indicator that measures the variation of the cost function $\mathcal{C}(f)$ in response to a small change of the flow $F_p$ on a path $p \in \mathcal{P}_w$. It is defined as follows:

$$C'_{p}(f) = \frac{\partial \mathcal{C}}{\partial F_{p}} \qquad (21)$$

To calculate the marginal cost, we firstly introduce the link cost function $C_e$. Given a link $e = (v, v')$, the associated link cost function $C_e$ is given as follows:

$$C_e(f) = (|\text{Re}(\lambda_v^d)| + |\text{Re}(\lambda_{i'}^e)|)\alpha_e f_e \qquad (22)$$

Obviously, the sum of link cost functions is equal to the total complexity of the network:

$$\mathcal{C}(f) = \sum_{e \in \mathcal{E}} C_e(f) \qquad (23)$$

Therefore, the marginal cost can be expressed as:

$$C'_{p}(f) = \sum_{a,b \in \mathcal{E}} \frac{\partial C_{a}(f)}{\partial f_{b}} \delta_{b,p} \qquad (24)$$

where $\delta_{b,p} = 1$ if link $b$ is contained in path $p$, otherwise $\delta_{b,p} = 0$.

In our problem, the dominant eigenvalue of LDS involved in the cost function is also dependent on the link flows. Thus, the derivative of the cost function does not have an analytical form. We use numerical differentiation to estimate the derivative on each link flow:

$$\frac{\partial C_a(f)}{\partial f_b} = \lim_{\epsilon \to 0} \frac{C_a(f_b + \epsilon) - C_a(f_b)}{\epsilon} \qquad (25)$$

A large value of marginal cost indicates that the flow on this path has a strong influence on the total complexity of the network and vice versa. The following proposition states that the equivalent marginal cost with positive path flow can lead to a system equilibrium solution to the optimization problem (P1) (Dafermos, 1971):
**Proposition 1.** F is a system equilibrium flow pattern if and only if for any OD pair w ∈ W with feasible routes $p_w^1, ..., p_w^{R_w}$, the marginal cost of these paths can be sorted to satisfy:

$$C'_{p_w}(f) = \dots = C'_{p_w^m}(f) \le C'_{p_w^{m+1}}(f) \le \dots \le C'_{p_w^{R_w}}(f) \qquad (26)$$

With

$$F_w^i > 0, \quad i = 1, \dots, m \qquad (27)$$

$$F_w^i = 0, \quad i = m + 1, \dots, R_w \qquad (28)$$

To balance the marginal cost, a certain quantity of flow should be transferred from the most expensive route $p_{w,min}$ to the least expensive route $p_{w,max}$ such that the total complexity is minimized, where

$$\begin{cases} p_{w,min} = \arg\min_{p \in P_w} C'_p(f) \\ p_{w,max} = \arg\max_{p \in P_w, F_p > 0} C'_p(f) \end{cases} \qquad (29)$$

The flow transfer operation is defined as an equilibrium operator $\mathcal{E}_w$:

$$\mathcal{E}_w : \mathbb{Z} \to \mathbb{Z} \qquad (30)$$

$$F_w \mapsto \mathcal{E}_w F_w \qquad (31)$$

Starting from an initial flow pattern $F_w^{(0)}$, the cost function converges through iterating $\mathcal{E}_w$:

$$F_w^{(n)} = \mathcal{E}_w F_w^{(n-1)}, \quad n \in \mathbb{N}^* \qquad (32)$$

New path flow pattern $\hat{F}_w$ is updated using the following rule:

$$\left\{ \begin{array}{l} \hat{F}_p = F_p, \quad \forall p \in (P_w - \{p_{w,min}, p_{w,max}\}) \\ \hat{F}_{p_{w,min}} = F_{p_{w,min}} + \sigma_w \\ \hat{F}_{p_{w,max}} = F_{p_{w,max}} - \sigma_w \end{array} \right. \qquad (33)$$

where $\sigma_w$ is the solution to the following optimization problem:

$$\begin{array}{ll} \min_{\sigma_w} & \mathcal{E}(\hat{f}) \\ \text{s.t.} & \hat{f}_e = f_e + (\delta_{e,p_{w,min}} - \delta_{e,p_{w,max}})\sigma_w, \quad e \in \mathcal{E} \\ & \sigma_w \in [-F_{p_{w,min}}, F_{p_{w,max}}] \end{array} \qquad (34)$$

For the efficiency of the algorithm, the constraint interval of $\sigma_w$ is uniformly discretized into a set of size $N_d$. The element $\sigma_w$ that minimizes the objective function is selected as the flow to transfer in each OD pair. The algorithm is stopped when there is no improvement on any of the OD pairs.

## 5. Experimental results

In this paper, the case study focuses on the urban airspace of Singapore. As a country with 100% urban population (Chizho et al., 2021), the development and implementation of UAM are prosperous in Singapore. Based on the study of Yixi et al. (2020) that determined the optimal number and location of droneport in Singapore under the estimated future demand of UAVs, a multi-layer air transport network is constructed as an example network covering most Singapore areas for delivery purposes. The construction of the network follows the following principles:

1. The network layer is sparser when the altitude increases. Since fewer obstacles exist at higher altitudes, nodes have more possibility to be connected directly (Pang et al., 2020).

2. A trade-off of numbers of layers. The fewer horizontal layers, the more constrained and complex the network. However, the computation can be more efficient. Conversely, more horizontal layers can provide more routing options and increased airspace capacity. The complexity level of the optimal flow pattern can be low, but the computation time will be increased.

The map view of the resulting network is plotted in Fig. 6. The graph representation of the transport network is shown in Fig. 7(a), which consists of 63 nodes, 110 links, 24 OD pairs, and 3 horizontal layers. The nodes are made up of 3 droneports, 5 delivery stations, and 55 waypoints. Adjacent layers have a separation of 15 m and are connected by vertical routes. The air traffic demand during a period is assigned to each OD pair. A given amount of UAVs carrying parcels fly from droneports to delivery stations, while the charged and maintained UAVs return to droneports from delivery stations. The two-way demands can be unbalanced, even for the same OD pair.

Fig. 7(b) depicts the complexity area and density points located along routes. The density points are located on each link and the width of the volume segment is expanded for illustration. The density points contained in each cylinder are the observations to measure the node complexity. The size of the volume segment is determined in the presence of uncertainty (positioning accuracy,
<figure><img src="image_7.png" /><figcaption>Fig. 6. A map-based visualization of air transport network in Singapore's urban airspace.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7. Graph representation of air transport network used in this study.</figcaption></figure>

safety margin, emergency response, etc.). According to the regulation of unmanned aircraft (CAAS, 2019), the parameters of the simulation module are determined, as shown in Table 2.

To evaluate the effectiveness of the proposed algorithm, two cases with different traffic demands are carried out for the framework. In the first case, the two-way demand of each OD pair is equal, and the demand gradually increases with the index of the OD pair. The second case simulates imbalanced stochastic demands. Figs. 8(a) and 8(b) plot the demand imposed on all OD pairs in these two cases.

To initialize the simulation module, the traffic demand is first randomly assigned to path $p \in P_w$ in the following way, where the constraints in problem (P1) should be satisfied and the random number generator is kept as the same for these two cases:

$$F_p = \frac{d_w \text{rand}_p}{\sum_{p \in P_w} \text{rand}_p} \qquad (35)$$
<table><caption>**Table 2**
Parameters setting.</caption><thead><tr><th>Module</th><th>Parameter</th><th>Value</th><th>Description</th></tr></thead><tbody><tr><td rowspan="10">Simulation module</td><td>L<sub>h</sub></td><td>0.3 (km)</td><td>Length of grid on horizontal volume segment</td></tr><tr><td>L<sub>v</sub></td><td>2 × 10<sup>-3</sup> (km)</td><td>Length of grid on vertical volume segment</td></tr><tr><td>L<sub>w</sub></td><td>0.03 (km)</td><td>Width of grid</td></tr><tr><td>s<sub>h</sub></td><td>120 (km/h)</td><td>Speed for horizontal link e<sub>h</sub></td></tr><tr><td>R<sub>max</sub></td><td>20 (km)</td><td>Maximum range of UAV</td></tr><tr><td>T<sub>max</sub></td><td>2000 (s)</td><td>Maximum endurance of UAV</td></tr><tr><td>s<sub>v</sub></td><td>5 (km/h)</td><td>Speed for vertical link e<sub>v</sub></td></tr><tr><td>N<sub>s</sub></td><td>3</td><td>Number of lateral points on each side</td></tr><tr><td>β</td><td>3</td><td>Efficiency threshold parameter</td></tr><tr><td>R<sub>ε</sub></td><td>5 (km)</td><td>Radius of complexity area for link e</td></tr><tr><td>F<sub>l</sub></td><td>10</td><td>Meaningful flow threshold e</td></tr><tr><td>φ<sub>h</sub></td><td>2 × 10<sup>-4</sup></td><td>Efficiency coefficient for horizontal links</td></tr><tr><td>φ<sub>s</sub></td><td>10<sup>-3</sup></td><td>Efficiency coefficient for upward links</td></tr><tr><td>φ<sub>d</sub></td><td>2 × 10<sup>-5</sup></td><td>Efficiency coefficient for downward links</td></tr><tr><td rowspan="4">Initialization module</td><td>T<sub>0</sub></td><td>1000</td><td>Initial temperature</td></tr><tr><td>N<sub>tr</sub></td><td>100</td><td>Number of iterations at each temperature</td></tr><tr><td>α</td><td>0.9</td><td>Cooling rate</td></tr><tr><td>stopping criterion</td><td>(y<sub>p</sub> - y<sub>c</sub>)/y<sub>p</sub> ≤ 1%</td><td>The best objective value of current temperature is improved less than 1%</td></tr><tr><td>Optimization module</td><td>ε</td><td>0.01</td><td>Small variation for numerical derivation</td></tr></tbody></table>

<figure><img src="image_9.png" /><figcaption>**Fig. 8.** The distribution of demand for OD pairs in two cases.</figcaption></figure>

## 5.1. Performance comparison between optimization algorithms

In this subsection, we report the experimental results regarding our proposed approach. In addition, a comparative study on the same simulation module is conducted between the proposed method and other widely used algorithms for traffic assignment, including all-or-nothing (AON) assignment, SA, DA, and Method of Successive Averages (MSA). DA has been applied without the SA initialization phase. AON assignment is the most basic algorithm for traffic assignment, which simply assigns all demands for each OD pair to the shortest path. MSA updates the flow pattern based on a descent direction and a decreasing step size (Mounce and Carey, 2015). The descent direction is extracted from the auxiliary path assignments obtained by the All-or-nothing assignment, and the step size is set as $1/(\lambda+1)$ in the iteration $\lambda$. MSA is also used as initialization for DA to compare with the proposed framework. For all algorithms except for the AON assignment, we run 100 simulation experiments with different random seeds, which control the initial value and the inherent randomness of algorithms. All algorithms are implemented in Java on a laptop with Intel® i7-8750H CPU and 32 GB DDR4 RAM.

To quantify the computational efficiency and performance of different models, Tables 3 and 4 illustrate the results in experiments of all aforementioned algorithms. Their 95% Confidence Intervals (CI) for standard normal distribution are calculated and expressed in the form of sample mean plus/minus margin of error. The best results are marked in bold.

Among all algorithms, AON assignment has the fastest computation time, but the complexity has barely decreased (roughly 24% in both cases). MSA is the second most efficient algorithm, and the resulting complexity is acceptable. However, if the result is used to initialize DA, DA can converge very slowly and get stuck in a local minimum, which corresponds to the result of MSA+DA in
<table><caption>**Table 3**
Performance comparison of models in terms of 95% CI for computation time, complexity, and flight efficiency in case 1.</caption><thead><tr><th rowspan="2">Models</th><th colspan="2">Case 1</th><th rowspan="2">Complexity reduction percentage</th><th rowspan="2">Average travel time (s)</th><th rowspan="2">Average path length (km)</th><th rowspan="2">Energy consumption</th></tr><tr><th>Computation time (s)</th><th>Total complexity</th></tr></thead><tbody><tr><td>Initial state</td><td>-</td><td>680081.00 ± 2000.58</td><td>-</td><td>404.59 ± 0.32</td><td>12.61 ± 0.01</td><td>228.06 ± 0.18</td></tr><tr><td>AON assignment</td><td>&lt; 10<sup>-4</sup></td><td>509223.03</td><td>25.12%</td><td><strong>248.00</strong></td><td><strong>7.16</strong></td><td><strong>130.27</strong></td></tr><tr><td>SA</td><td>5.86 ± 0.26</td><td>192021.02 ± 1613.07</td><td>71.76%±0.25%</td><td>307.27 ± 0.38</td><td>10.11 ± 0.01</td><td>182.72 ± 0.27</td></tr><tr><td>DA</td><td>38.58 ± 0.79</td><td>98957.43 ± 2121.67</td><td>85.45%±0.31%</td><td>320.10 ± 1.83</td><td>10.07 ± 0.06</td><td>181.11 ± 1.15</td></tr><tr><td>MSA</td><td>0.92 ± 0.02</td><td>228767.25 ± 524.08</td><td>66.36%±0.12%</td><td>309.19 ± 0.07</td><td>10.17 ± 0.00</td><td>183.26 ± 0.04</td></tr><tr><td>MSA+DA</td><td>28.55 ± 1.44</td><td>159262.94 ± 8523.64</td><td>76.58%±1.26%</td><td>316.88 ± 1.48</td><td>10.32 ± 0.05</td><td>185.68 ± 0.84</td></tr><tr><td>Proposed model</td><td>13.62 ± 0.52</td><td><strong>65005.20 ± 3552.16</strong></td><td><strong>90.44%±0.53%</strong></td><td>323.83 ± 2.52</td><td>10.27 ± 0.08</td><td>184.68 ± 1.43</td></tr></tbody></table>

<table><caption>**Table 4**
Performance comparison of models in terms of 95% CI for computation time, complexity, and flight efficiency in case 2.</caption><thead><tr><th rowspan="2">Models</th><th colspan="2">Case 2</th><th rowspan="2">Complexity reduction percentage</th><th rowspan="2">Average travel time (s)</th><th rowspan="2">Average path length (km)</th><th rowspan="2">Energy consumption</th></tr><tr><th>Computation time (s)</th><th>Total complexity</th></tr></thead><tbody><tr><td>Initial state</td><td>-</td><td>665702.86 ± 2017.43</td><td>-</td><td>405.60 ± 0.30</td><td>12.66 ± 0.01</td><td>220.04 ± 0.17</td></tr><tr><td>AON assignment</td><td>&lt; 10<sup>-4</sup></td><td>510435.59</td><td>23.32%</td><td><strong>258.21</strong></td><td><strong>7.48</strong></td><td><strong>131.06</strong></td></tr><tr><td>SA</td><td>6.29 ± 0.34</td><td>183251.32 ± 1694.38</td><td>72.47%±0.27%</td><td>314.62 ± 0.29</td><td>10.34 ± 0.01</td><td>179.66 ± 0.27</td></tr><tr><td>DA</td><td>35.31 ± 0.89</td><td>100567.56 ± 2869.64</td><td>84.89%±0.80%</td><td>323.41 ± 1.55</td><td>10.13 ± 0.05</td><td>175.01 ± 0.55</td></tr><tr><td>MSA</td><td>0.98 ± 0.03</td><td>216930.25 ± 686.09</td><td>67.41%±0.15%</td><td>317.10 ± 0.08</td><td>10.45 ± 0.00</td><td>180.89 ± 0.06</td></tr><tr><td>MSA+DA</td><td>27.53 ± 0.87</td><td>154941.17 ± 3293.94</td><td>76.73%±0.50%</td><td>319.13 ± 0.82</td><td>10.48 ± 0.03</td><td>183.58 ± 0.45</td></tr><tr><td>Proposed model</td><td>14.12 ± 0.33</td><td><strong>52450.41 ± 2376.11</strong></td><td><strong>92.12%±0.35%</strong></td><td>328.58 ± 1.84</td><td>10.39 ± 0.06</td><td>179.49 ± 1.11</td></tr></tbody></table>

Tables 3 and 4. DA is the most time-consuming algorithm. However, if DA is initialized by SA, which is the proposed framework, the convergence speed can be significantly improved. DA only takes 7.76 ± 0.43 and 7.83 ± 0.39 s and the CPU time is reduced by 82.02% ± 1.14% and 80.47% ± 1.34%, respectively. Besides, the total complexity is also significantly decreased by 90.44% ± 0.53% and 92.12% ± 0.35% compared to the initial state in the two cases, respectively. In summary, the proposed two-phase algorithm outputs the flow pattern with the lowest complexity among all algorithms in a reasonable time, which supports the primary objective of this study.

Besides, the flight efficiency in the air transport network is also measured to ensure that it is at an acceptable level. In addition to energy consumption defined in Eq. (10), two indicators are introduced, including average travel time $\bar{i}$:

$$ \bar{i} = \sum_{w \in W} \sum_{p \in P_w} \frac{F_p}{D_w} \sum_{e \in P_e} \frac{L_e}{s_e} \qquad (36) $$

and average path length $\bar{l}$:

$$ \bar{l} = \sum_{w \in W} \sum_{p \in P_w} \frac{F_p}{D_w} \sum_{e \in P_e} L_e \qquad (37) $$

where $l_e$ is the length of link $e$.

The flight efficiency metrics in terms of initial flow patterns and the results of different algorithms are also shown in Tables 3 and 4. The best results are marked in bold. As can be seen, AON assignment outperforms other methods in terms of flight efficiency because the flow is entirely allocated to the shortest path of each OD pair. Besides, the flight efficiency results of other algorithms are markedly improved compared with the flight efficiency of the initial flow pattern. The results reveal that the flow pattern obtained by our proposed model also results in acceptable flight efficiency and operational efficiency.

Next, the evolution of objection function values will be compared between the aforementioned algorithms given the same random seed, which is randomly selected from the 100 experiments. The complexity evolution of the proposed approach is firstly illustrated in Fig. 9. In both cases, SA reduces the total complexity by more than 70%. By applying DA to further optimize the problem, the complexities can still be decreased with only a few iterations. It is a positive result in a view that the initial complexities are reduced by 92.45% and 93.29% in two cases, respectively.

To view the details of the algorithm that produces the lowest complexity flow pattern, the evolution of total complexity is separately analyzed in the two phases, see Figs. 10 and 11. The convergence curve of SA is illustrated in Figs. 10(a) and 11(a). The curve is segmented by temperature steps, and the minimum complexities for each temperature are marked with red dots. A rapid decrease can be seen in the first temperature step. Although the objective function drops slightly in the following temperature steps, the complexity of paths is still being balanced and refined, within a few iterations. The complexity is relatively low when the stopping criterion of SA is met. Then, DA starts to further optimize the problem. The evolution of complexity in terms of DA is illustrated in Figs. 10(b) and 11(b). The convergence speed of DA is quite fast: the objective function value decreases significantly
<figure><img src="image_10.png" /><figcaption>Fig. 9. Complexity evolution of the proposed two-phase algorithm.</figcaption></figure>

<figure><img src="image_11.png" /></figure>

<figure><img src="image_12.png" /><figcaption>Fig. 10. Complexity evolution for two phases of the proposed algorithm in the first case.</figcaption></figure>

<figure><img src="image_13.png" /></figure>

<figure><img src="image_14.png" /><figcaption>Fig. 11. Detailed complexity evolution for two phases of the proposed algorithm in the second case.</figcaption></figure>

<figure><img src="image_15.png" /></figure>
<figure><img src="image_16.png" /><figcaption>Fig. 12. Complexity evolution of DA without initialization phase in two cases.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Fig. 13. Complexity evolution of DA initialized by MSA in two cases.</figcaption></figure>

in less than 100 iterations. This fact demonstrates that SA with a probabilistic performance-based neighborhood generation strategy can greatly improve the convergence performance of DA.

The evolution of complexity in terms of DA with random initialization is shown in Fig. 12. Although the objective is reduced to a low level, it is almost twice the complexity of the flow pattern obtained by our proposed algorithm, as shown in Figs. 9(a) and 9(b). In addition, a large number of iterations is required for DA without appropriate initial values. This fact can also be observed in Fig. 13, which depicts the complexity evolution of DA initialized by MSA in two cases. Although it appears that MSA provides low-complexity results with few iterations, DA quickly falls into local optima initialized by this result, which also corresponds to the results in Tables 3 and 4.

## 5.2. Path flow assignment

In order to better interpret the results of the proposed model, the related results consistent with the experiments in 9(a) and 9(b) are discussed in the air transport network representation. To start with, the path flow allocation for the two cases are set out in Figs. 14 and 15. The OD pairs are divided by light gray dotted lines and are partitioned by a bold gray dotted line according to the flow direction. The color indicates the path complexity. A logarithmic scale colormap is used to amplify the difference of values, which is also adopted in the following figures. As shown in Figs. 14(a) and 15(a), the flow pattern assigned randomly brings high complexity on almost all paths. Optimized by SA, the flow is allocated to fewer paths but with better quality in each OD pair, which is illustrated in Figs. 14(b) and 15(b). The complexity of paths associated with less than half OD pairs is sharply decreased to an acceptable level. However, in both cases, for most OD pairs, the problem lies in that two-way flows are mainly routed on the same path, which is not an optimal solution. These paths are congested, while other alternative paths are unused. In addition, large bidirectional flow on the same link brings doubled complexity to connected nodes. The link capacity is also under intense pressure. The phenomenon is especially obvious in the first case with symmetric two-way demand (Fig. 14(b)). Refined by DA in the second phase of the proposed model, the final path flow allocation result is presented in Figs. 14(c) and 15(c). The path complexity is
<figure><img src="image_18.png" /><figcaption>Fig. 14. Evolution of path flow and path complexity of the first case. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

reduced to a low level. Besides, the two-way demand is prevented from being assigned on the same path for each OD pair. This fact can be seen clearly in the first case (in Fig. 14(c)). For each OD pair, the flow in different directions are assigned to a different path, and the path complexity is much lower.

### 5.3. Link flow allocation

In the network representation, the link flow evolution is shown in Figs. 16 and 17. The link color indicates the amount of flow allocated on each link. As can be seen from part (a) of these figures, nearly all links are assigned with high flows, which brings unnecessary congestion and difficulty to manage traffic on waypoints. Optimized by our model, the link flow allocation has been significantly changed, as shown in 16(b) and 17(b). The main flow is allocated to vertical links and some inner links of each network layer. A few flows are found on marginal links of each network layer. As expected, the flow is distributed on different layers of the network, which reduces the overall air traffic complexity.

### 5.4. Node complexity

Under previously mentioned flow allocation result, Figs. 18 and 19 represents the evolution of node complexity in the two cases. The color is related to the node complexity. The initial complexity of nodes are plotted in Figs. 18(a) and 19(a). High complexities with a scale of $10^5$ can be found on the nodes in the central part of network layers and the nodes connecting vertical links. After optimization, the node complexity is very low. Even though vertical links are assigned with large amounts of flow, their connecting nodes have low complexities. In the network, the main complexity can be found at low-level layers, which is consistent with the result of link flow allocation in Figs. 16(b) and 17(b), that the low-level layer is the most congested in the network. Given limited airspace capacity and high demands, the proposed algorithm produces a good quality solution. The insight gained from this study may be of assistance to Air Navigation Service Provider (ANSP) for complexity-optimal strategic planning of UAM operations in the air transport network.

## 6. Conclusion

In this paper, we introduced an efficient methodology to address the air traffic assignment problem in a future urban environment with intensive UAM operations. As a strategic planning method, it models UAM operations as air traffic flows, which operate within fixed-route structures referred to as volume segments. The air transport network is modeled as a graph, with volume segments as
<figure><img src="image_19.png" /><figcaption>Fig. 15. Evolution of path flow and path complexity of the second case. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>Fig. 16. Link flow evolution of the first case in the network representation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

links; droneports, delivery stations, and waypoints as nodes. To make better usage of airspace capacity, vertical links are added to connect the horizontal layers and all links support two-way traffic. Density points are introduced to model the operational density along each link. Cylindrical airspace is designed to capture the complexity in terms of each node. LDS incorporates temporal and spatial information to measure the local disorder and interaction of traffic. To minimize the complexity of UAM operations in the network, a macroscopic air traffic assignment problem is formulated. The objective function considers dynamic flow structures, link congestion, and operational efficiency. We introduce a two-phase algorithm on the basis of SA and DA to optimize the UAM flow pattern. The part of SA rapidly generates candidate solutions by balancing the path complexity in each OD pair, while the part of DA intensifies the search and accelerates the convergence on the basis of the candidate solution.

Computational experiments are conducted on a multi-level two-way air transport network over the Singapore urban airspace given two different traffic demands. In addition, baselines including AON assignment, DA, SA, MSA, and DA initialized by MSA are
<figure><img src="image_21.png" /><figcaption>Fig. 17. Link flow evolution of the second case in the network representation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>Fig. 18. Node complexity evolution of the first case in the network representation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>Fig. 19. Node complexity evolution of the second case in the network representation. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

compared with the proposed model based on the same simulation module. The numerical experiments and comparative studies in terms of algorithmic performance, flight efficiency, path flow assignment, link flow allocation, and node complexity demonstrate that the proposed model can significantly reduce the complexity of air transport networks. The flow is efficiently allocated and the network capacity is fully used.
Due to operational and geographical restrictions, the selection and construction of delivery sites and the air transport network can be limited in Singapore. With the increasing density of UAM operations, these findings are of assistance to ANSP for complexity-optimal strategic planning of UAM operations, which is envisioned in the low-altitude urban area in class G airspace. Moreover, the parameters of the model can be easily adjusted according to the actual situation. Although the methodology defined here is applied to UAM operations in Singapore's urban airspace, it could be applied to other metropolitan areas in the world.

Future research would focus on dynamical traffic assignment with time-varying demand. A multi-objective optimization framework would be introduced to additionally minimize noise, delay, etc. Besides, we plan to investigate efficient algorithms for air traffic assignment with acceleration strategies such as parallel computing. In addition, extension to multi-class traffic assignments and large-scale air transport networks could be an interesting topic for future research.

### CRediT authorship contribution statement

**Zhengyi Wang:** Conceptualization, Methodology, Software, Formal analysis, Visualization, Writing – original draft. **Daniel Delahaye:** Conceptualization, Supervision, Writing – review & editing, Funding acquisition. **Jean-Loup Farges:** Conceptualization, Writing – review & editing, Supervision, Funding acquisition. **Sameer Alam:** Supervision, Writing – review & editing, Conceptualization.

### Acknowledgments

This work is partially supported by the research project CONCORDE of the Defense Innovation Agency (AID) of the French Ministry of Defense (2019650090004707501). The authors would like to acknowledge the help in the review of this paper by the late Serge Roux from ENAC. The authors would also like to thank the editor and anonymous reviewers for their valuable comments.

## Appendix A. Weighted minimum mean square error estimation of linear dynamical system in three dimensions

Within a continuous-time model scope, the LDS has the following form:

$$ \dot{\mathbf{x}} = \mathbf{A}\mathbf{x} + \mathbf{b} \qquad (A.1) $$

In three dimensions, given $N$ density points with positions $\mathbf{x}_i = (x_i, y_i, z_i)^\top$, speeds $\mathbf{v}_i = (v_{x,i}, v_{y,i}, v_{z,i})^\top$, and densities $\omega_i$, $i = 1, \dots, N$, the Weighted Minimum Mean square Error (WMMSE) criterion between the LDS model and the observations is given as:

$$ E = \sum_{i=1}^{N} \omega_i \| \mathbf{v}_i - (\mathbf{A}\mathbf{x}_i + \mathbf{b}) \|^2 \qquad (A.2) $$

In order to transform Eq. (A.2) into matrix form, the following matrices are defined:

$$ \mathbf{X} = \begin{bmatrix} \sqrt{\omega_1}x_1 & \sqrt{\omega_2}x_2 & \cdots & \sqrt{\omega_N}x_N \\ \sqrt{\omega_1}y_1 & \sqrt{\omega_2}y_2 & \cdots & \sqrt{\omega_N}y_N \\ \sqrt{\omega_1}z_1 & \sqrt{\omega_2}z_2 & \cdots & \sqrt{\omega_N}z_N \\ \sqrt{\omega_1} & \sqrt{\omega_2} & \cdots & \sqrt{\omega_N} \end{bmatrix} \qquad (A.3) $$

$$ \mathbf{V} = \begin{bmatrix} \sqrt{\omega_1}v_{x,1} & \sqrt{\omega_2}v_{x,2} & \cdots & \sqrt{\omega_N}v_{x,N} \\ \sqrt{\omega_1}v_{y,1} & \sqrt{\omega_2}v_{y,2} & \cdots & \sqrt{\omega_N}v_{y,N} \\ \sqrt{\omega_1}v_{z,1} & \sqrt{\omega_2}v_{z,2} & \cdots & \sqrt{\omega_N}v_{z,N} \end{bmatrix} \qquad (A.4) $$

$$ \mathbf{A} = \begin{bmatrix} a_{11} & a_{12} & a_{13} \\ a_{21} & a_{22} & a_{23} \\ a_{31} & a_{32} & a_{33} \end{bmatrix} \qquad (A.5) $$

$$ \mathbf{b} = \begin{bmatrix} b_1 & b_2 & b_3 \end{bmatrix}^\top \qquad (A.6) $$

$$ \mathbf{C} = [\mathbf{A} | \mathbf{b}] \qquad (A.7) $$

Then, the error criterion $E$ can be reformulated as

$$ E = \|\mathbf{V} - \mathbf{C}\mathbf{X}\|_F^2 \qquad (A.8) $$

where $\|\cdot\|_F$ represents the Frobenius norm.

The WMMSE aims to find the matrix $\hat{\mathbf{C}}$ that minimize the error criterion:

$$ \hat{\mathbf{C}} = \arg \min_{\mathbf{C} \in \mathbb{R}^{3 \times 4}} E \qquad (A.9) $$

To this end, the gradient of $E^2$ is calculated as:

$$ \nabla_{\mathbf{C}} E = -2 (\mathbf{V} - \mathbf{C}\mathbf{X}) \mathbf{X}^\top \qquad (A.10) $$
* If $\mathrm{Re}(\lambda) > 0$, then $\lim_{t \to +\infty} x(t) = \infty$, which means the component is unstable and diverging,

* If $\mathrm{Re}(\lambda) < 0$, then $\lim_{t \to +\infty} x(t) = 0$, which means the component is stable and converging,

* If $\mathrm{Re}(\lambda) = 0$, the component is conserved. Depending on the imaginary part, it can be further classified into:

- If $\mathrm{Im}(\lambda) \neq 0$, the component is rotating,

- If $\mathrm{Im}(\lambda) = 0$, the component is in translation.

**Proof.** In order to generalize all forms of LDS by a constant-free representation, let $y = \begin{bmatrix} x \\ 1 \end{bmatrix}$, then Eq. (A.1) can be expressed as follows by adding one more dimension to the system:

$$ \dot{y} = \begin{bmatrix} \dot{x} \\ 1 \end{bmatrix} \qquad (B.1) $$

$$ = \begin{bmatrix} \mathbf{A} & \mathbf{b} \\ 0 & 0 \end{bmatrix} \begin{bmatrix} x \\ 1 \end{bmatrix} \qquad (B.2) $$

$$ = \mathbf{B}y \qquad (B.3) $$

A closed-form solution of B.1 for any square matrix B is:

$$ y(t) = e^{\mathbf{B}t} y_0 \qquad (B.4) $$

Denote the dimension of B as N. B is assumed to be diagonalizable with independent eigenvectors $v_1, \dots, v_N$ and the associated eigenvalues $\lambda_1, \dots, \lambda_N$, then the initial state of the system can be represented as:

$$ y_0 = \sum_{i=1}^{N} b_i v_i \qquad (B.5) $$

By applying the following property,

$$ e^{\mathbf{B}} v_i = e^{\lambda_i} v_i \qquad (B.6) $$

the state can be formulated as:

$$ x(t) = e^{\mathbf{B}t} \sum_{i=1}^{N} b_i v_i \qquad (B.7) $$

$$ = \sum_{i=1}^{N} b_i v_i e^{\lambda_i t} \qquad (B.8) $$

Thus, the corresponding eigenvector component of a particular eigenvalue $\lambda_i \in \{\lambda_1, \dots, \lambda_N\}$ is $b_i v_i e^{\lambda_i t}$. Its property is consistent with Proposition 2. $\square$

## Appendix C. Representative qualitative behavior of LDS characterized by the dominant eigenvalue

### Proposition 3. The dominant eigenvalue of A determines the stability of the LDS.

**Proof.** Following the constant-free form of LDS in B.1, let $\hat{\lambda}$ be the dominant eigenvalue of B, then Eq. (B.7) can be written as:

$$ x(t) = e^{\hat{\lambda}t} \sum_{i=1}^{N} b_i v_i e^{(\lambda_i - \hat{\lambda})t} \qquad (C.1) $$

Because the real part of $\lambda_i - \hat{\lambda}$ is negative, we have:

$$ \lim_{t \to +\infty} x(t) \approx b_k v_k e^{\hat{\lambda}t} \qquad (C.2) $$

According to Eq. (C.2), A dominant eigenvalue $\hat{\lambda}$ determines the stability of the whole LDS. $\square$

## Appendix D. Eigenvalue loci for typical three-dimensional traffic situations

In Fig. D.20, four typical scenarios of three-dimensional traffic situations are presented. These cases include translation, rotation, convergence, and divergence. All traffic situations can be seen as the combination of these typical scenarios. The location of the observed traffic is indicated by the black dots and speed vectors are indicated by red arrows. The traffic situation estimated by LDS is denoted by blue arrows. The arrow head indicates the direction and the length of the arrow stands for speed. The eigenvalues associated to the coefficient matrices A are given by the LDS and the corresponding loci are presented in the complex plane.

In the translation case, the speed vectors are parallel and the relative distances between aircraft remain unchanged. The associated eigenvalues are zero. With regard to the rotation case, the aircraft moves around a circle but the relative distances between aircraft
<figure><img src="image_24.png" /><figcaption>Fig. D.20. Eigenvalue loci for four typical traffic situations in three-dimensional airspace. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)</figcaption></figure>

remain constant. All associated eigenvalues have a zero real part. The eigenvalues associated with these two cases are located within the strip in Fig. 3, which corresponds to organized traffic situations.

When the relative distances between aircraft diminish with time, which corresponds to the convergence case, the real part of non-zero eigenvalues is negative, and the LDS is in contraction mode. Finally, in the divergence case, those relative distances increase with time. The real part of non-zero eigenvalues is positive and the LDS is in expansion mode. The aircraft are diverging and their relative distances increase with time. The third and the fourth case belong to unorganized traffic patterns because the relative distance between aircraft varies with time. It is worth mentioning that the third case is more critical because aircraft are converging to each other. In addition, the larger the absolute value of the real part of eigenvalues, the faster the evolution of LDS.

In conclusion, when the relative distances between aircraft change slowly with time, and the relative speeds between aircraft are close to zero, the eigenvalues of coefficient matrix A of the associated LDS have a small real part (Delahaye et al., 2004). This corresponds to organized traffic situation where the traffic has no or few interaction.

References

Bauranov, Aleksandar, Rakas, Jasenka, 2021. Designing airspace for urban air mobility: A review of concepts and approaches. Prog. Aerosp. Sci. 125, 100726.
http://dx.doi.org/10.1016/j.paerosci.2021.100726.

CAAS, 2019. Air navigation (101 — Unmanned aircraft operations) regulations 2019. https://www.caas.gov.sg/legislation-regulations/legislation/air-navigation-act.

Chizho, Larisa N, Okuneva, Tatyana D, Hillawi, Gharbawee Mohammed Anwer, 2021. Sustainable Urban development models. In: Sustainable Development of Modern Digital Economy. Springer, pp. 211–222.

Dafermos, Stella C., 1971. An extended traffic assignment model with applications to two-way traffic. Transp. Sci. 5 (4), 366–389. http://dx.doi.org/10.1287/trsc.5.4.366.

Dafermos, Stella C., 1972. The traffic assignment problem for multiclass-user transportation networks. Transp. Sci. 6 (1), 73–87. http://dx.doi.org/10.1287/trsc.6.1.73.

Dafermos, Stella C., Sparrow, Frederick T., 1969. The traffic assignment problem for a general network. J. Res. Natl. Bureau Stand. B 73 (2), 91–118.

Delahaye, Daniel, Chaimatanan, Supatcha, Mongeau, Marcel, 2019. Simulated annealing: From basics to applications. In: Gendreau, Michel, Potvin, Jean-Yves (Eds.), Handbook of Metaheuristics. In: International Series in Operations Research & Management Science (ISOR), vol. 272, Springer, pp. 1–35.
http://dx.doi.org/10.1007/978-3-319-91086-4_1, URL https://hal-enac.archives-ouvertes.fr/hal-01887543.

Delahaye, Daniel, Paimblanc, Philippe, Puechmorel, Stephane, Histon, Jonathan M., Hansman, Robert J., 2002. A new air traffic complexity metric based on dynamical system modelization. In: Proceedings. the 21st Digital Avionics Systems Conference. Vol. 1. IEEE, p. 4A2. http://dx.doi.org/10.1109/dasc.2002.1067960.

Delahaye, Daniel, Puechmorel, Stéphane, 2000. Air traffic complexity: towards intrinsic metrics. In: Proceedings of the Third USA/Europe Air Traffic Management R & D Seminar.

Delahaye, Daniel, Puechmorel, Stephane, 2010. Air traffic complexity based on dynamical systems. In: 49th IEEE Conference on Decision and Control. CDC, IEEE, http://dx.doi.org/10.1109/cdc.2010.5718004.

Delahaye, Daniel, Puechmorel, Stéphane, 2013. Modeling and optimization of air traffic. Model. Optim. Air Traffic 1–329. http://dx.doi.org/10.1002/9781118743805.

Delahaye, Daniel, Puechmorel, Stephane, Hansman, John, Histon, Jonathan, 2004. Air traffic complexity map based on non linear dynamical systems. Air Traffic Control Q. 12 (4), 367–388. http://dx.doi.org/10.2514/atcq.12.4.367.
Delahaye, Daniel, Sofiane, Oussedik, Puechmorel, Stephane, 2005. Airspace congestion smoothing by multi-objective genetic algorithm. In: Proceedings of the 2005 ACM Symposium on Applied Computing - SAC'05. ACM Press, pp. 907-912. http://dx.doi.org/10.1145/1066677.1066887.

Deschinkel, Karine, Farges, Jean-Loup, Delahaye, Daniel, 2002. Optimizing and assigning price levels for air traffic management. Transp. Res. E Logist. Transp. Rev. 38 (3-4), 221-237. http://dx.doi.org/10.1016/s1366-5545(02)00007-8.

Essén, Hanno, 1993. Average angular velocity. Eur. J. Phys. 14 (5), 201-205. http://dx.doi.org/10.1088/0143-0807/14/5/002.

FAA, 2020. Concepts of Operations V2.0, Foundational Principles, Roles and Responsibilities, Scenarios and Operational Threads. Technical Report, U.S. department of transportation.

García, Adrían, Delahaye, Daniel, Soler, Manuel, 2020. Air Traffic Complexity Map based on Linear Dynamical Systems. Working Paper or Preprint.

Geister, Dagi, 2017. Concept for Urban Airspace Integration DLR U-Space Blueprint. Technical Report, Institute of Flight Guidance.

Histon, Jonathan M., Hansman, R. John, Aigoin, Guillaume, Delahaye, Daniel, Puechmorel, Stephane, 2002. Introducing structural considerations into complexity metrics. Air Traffic Control Q. 10 (2), 115-130. http://dx.doi.org/10.2514/atcq.10.2.115.

Ho-Huu, Vinh, Ganić, Emir, Hartjes, Sander, Babić, Obrad, Curran, Richard, 2019. Air traffic assignment based on daily population mobility to reduce aircraft noise effects and fuel consumption. Transp. Res. D Trans. Environ. 72, 127-147. http://dx.doi.org/10.1016/j.trd.2019.04.007.

Jahn, Olaf, Möhring, Rolf H., Schulz, Andreas S., 2000. Optimal routing of traffic flows with length restrictions in networks with congestion. In: Operations Research Proceedings 1999. Springer Berlin Heidelberg, pp. 437-442. http://dx.doi.org/10.1007/978-3-642-58300-1_68.

Jang, Dae-Sung, Ippolito, Corey A., Sankararaman, Shankar, Stepanyan, Vahram, 2017. Concepts of airspace structures and system analysis for UAS traffic flows for Urban areas. In: AIAA Information Systems-AIAA Infotech @ Aerospace. American Institute of Aeronautics and Astronautics, http://dx.doi.org/10.2514/6.2017-0449.

Juntama, Paveen, Chaimatanan, Supatcha, Alam, Sameer, Delahaye, Daniel, 2020. A distributed metaheuristic approach for complexity reduction in air traffic for strategic 4D trajectory optimization. In: 2020 International Conference on Artificial Intelligence and Data Analytics for Air Transportation. AIDA-AT, IEEE, pp. 1-9. http://dx.doi.org/10.1109/aida-at48540.2020.9049200.

Lascara, Brock, Lacher, Andrew, DeGarmo, Matthew, Maroney, David, Niles, Rick, Vempati, Lakshmi, 2019. Urban Air Mobility Airspace Integration Concepts. Technical Report, The MITRE Corporation.

Metropolis, Nicholas, Rosenbluth, Arianna W., Rosenbluth, Marshall N., Teller, Augusta H., Teller, Edward, 1953. Equation of State Calculations by Fast Computing Machines. Technical Report, Office of Scientific and Technical Information (OSTI), http://dx.doi.org/10.2172/4390578.

Migliore, Michele, Martorana, Vincenzo, Sciortino, Francesco, 1990. An algorithm to find all paths between two nodes in a graph. J. Comput. Phys. 87 (1), 231-236. http://dx.doi.org/10.1016/0021-9991(90)90235-s.

Ming, Cher (Ed.), 2008. Simulated Annealing. InTech, http://dx.doi.org/10.5772/67.

Mounce, Richard, Carey, Malachy, 2015. On the convergence of the method of successive averages for calculating equilibrium in traffic networks. Transp. Sci. 49 (3), 535-542.

Nosedal, Jenaro, Piera, Miquel A., Ruiz, Sergio, Nosedal, Alvaro, 2014. An efficient algorithm for smoothing airspace congestion by fine-tuning take-off times. Transp. Res. C 44, 171-184. http://dx.doi.org/10.1016/j.trc.2014.03.017.

Pang, Bizhao, Dai, Wei, Ra, Thu, Low, Kin Huat, 2020. A concept of airspace configuration and operational rules for UAS in current airspace. In: 2020 AIAA/IEEE 39th Digital Avionics Systems Conference. DASC, IEEE, pp. 1-9.

Schneider, Oliver, Kern, Stefan, Knabe, Franz, Gerdes, Ingrid, Delahaye, Daniel, Vidosavljevic, Andrija, van Leeuwen, Pim, Nieuwenhuisen, Dennis, Sunil, Emmanuel, Hoekstra, Jacco, Ellerbroek, Joost, 2015. Metropolis - Urban Airspace Design. Technical Report, The European Commission, pp. 1-56.

Singireddy, Shiva Ram Reddy, Daim, Tugrul U., 2018. Technology roadmap: Drone delivery - amazon prime air. In: Innovation, Technology, and Knowledge Management. Springer International Publishing, pp. 387-412. http://dx.doi.org/10.1007/978-3-319-68987-6_13.

Sunil, Emmanuel, Hoekstra, JM, Ellerbroek, Joost, Bussink, Frank, Nieuwenhuisen, Dennis, Vidosavljevic, Andrija, Kern, Stefan, 2015. Metropolis: Relating airspace structure and capacity for extreme traffic densities. In: Proceedings of the 11th USA/Europe Air Traffic Management Research and Development Seminar (ATM2015), Lisbon (Portugal), 23-26 June, 2015. FAA/Eurocontrol.

Thippavong, David P, Apaza, Rafael, Barmore, Bryan, Battiste, Vernol, Burian, Barbara, Dao, Quang, Feary, Michael, Go, Susie, Goodrich, Kenneth H, Homola, Jeffrey, et al., 2018. Urban air mobility airspace integration concepts and considerations. In: 2018 Aviation Technology, Integration, and Operations Conference. p. 3676.

Treimuth, Tambet, Delahaye, Daniel, Puechmorel, Stephane, Ngueveu, Sandra Ulrich, 2015. Parallel complexity computation based on dynamical systems. In: 2015 IEEE/AIAA 34th Digital Avionics Systems Conference. DASC, IEEE, http://dx.doi.org/10.1109/dasc.2015.7311340, 1C2-1.

UN, 2018. World urbanization prospects. https://population.un.org/wup/Download. (Accessed 15 2022).

Wang, Zhengyi, Delahaye, Daniel, Farges, Jean-Loup, Alam, Sameer, 2021. Air traffic assignment for intensive urban air mobility operations. J. Aerosp. Inf. Syst. 1-15. http://dx.doi.org/10.2514/1.i010954.

Yixi, Zeng, Huat, Low K., Duong, Vu N., Schultz, Michael, 2020. Droneport placement optimization and capacity prediction. In: 9th International Conference on Research in Air Transportation (ICRAT), Florida (Virtual Format), USA.