Prof. Fedja Netjasov

Prof Karine Deschinkel

Prof Wei Peng

Prof Pierre Maréchal

Serge Roux

Mohammed Sbihi

Hélène Weiss

Catherine Migot

Hasna Habouchi-Torchi

Andrija Vidosavljevic

Nicolas Durand

Gilles Baroin

Celine Demouge

Dinh- Thinh Hoang

Clara Buire

Geoffrey Scozzaro

Pierre Dieumegard

Sylvain Roudière

Paveen Juntama

Andreas Guitart

Sana Ikli

Gabriel Jarry

Philippe Monmousseau

Ying Huo

Xiao Hu

Shangrong Chen

Ruohao Zhang

Large numbers of Urban Air Mobility (UAM) vehicles are expected to operate in urban airspace in the near future. However, this growth will soon exceed the capacities of current airspace and air traffic management systems and put strains on transportation infrastructure, resulting in undesirable outcomes such as traffic congestion, complex traffic situations, delays, etc. As a result, innovative UAM solutions are required to ensure safe and efficient urban transportation. In this thesis, we focus on strategic planning in unmanned aircraft system traffic management, including airspace organization and management and traffic flow management. Specifically, we develop approaches for UAM route network design and air traffic assignment.

The structure of the UAM route network strongly impacts many aspects of the operating traffic flow. The UAM route network is expected to be well-designed for the desired network operations, especially with dense traffic demands. In order to minimize the impact on existing air traffic management systems, this thesis proposes a comprehensive framework to design a UAM route network in low-altitude urban airspace in the presence of obstacles and hazardous airspace. Using a variety of open-source data, this approach is applied as a case study for parcel delivery services using UAVs as UAM vehicles in Singapore's urban airspace. Firstly, The UAM route network is designed as a grid-based network that avoids obstacles and airspace within which the aircraft is prohibited. Then, the link cost is defined in terms of noise impact on populations, airspace safety, and flight efficiency. Unlike most research that only calculates the shortest path, we formulate a k-shortest path problem with diversity to select feasible routes between origin-destination pairs that minimize the route cost. The feasibility of this approach is demonstrated in the experiment, which can also be easily transferred to other scenarios. The impact of different parameter settings for link costs on UAM services is also explored. The feasible routes with low similarity provide more travel options, which can be used as pre-computed routes for air traffic assignment models.

To adapt the increasing demand to the current airspace capacity, we proposed several air traffic assignment models, including two static air traffic assignment models for single-layer and multi-layer two-way UAM route networks and a dynamic air traffic assignment model for multi-layer two-way UAM route networks in order to mitigate the congestion and complexity and organize the structure of air traffic flow. Firstly, at the macroscopic level, UAM operations are modeled as flows by aggregating the individual vehicle dynamics to describe the overall flow features with respect to the dense traffic volume. The flows are structured into UAM route networks, in which the air routes are modeled as corridors or volume segments. Then, the air traffic assignment model is formulated as an optimization problem. The objective functions are modeled to describe the air traffic complexity based on the linear dynamical system and the congestion based on energy consumption and traffic density. The simulation-based framework including i I would like to express my sincere gratitude to my supervisor, Prof. Daniel Delahaye. I first met Daniel during a two-month summer internship in 2017 recommended by Dr. Liang. Daniel offered me an opportunity to complete a Master of Science in ENAC majoring in operation research. During this period, I also completed a PFE internship in the Optim Group under his supervision. Daniel has always been a mentor for me since my postgraduate study. At the beginning of my Ph.D, when I was not clear about the

ii different optimization approaches is proposed to efficiently solve these problems. Computational experiments are performed on case studies of UAM route networks at various scales. The comparisons with regard to conventional traffic assignment algorithms are presented. The results show that the proposed approach is capable of assigning flows in an efficient and effective manner, significantly reducing the complexity and congestion of the UAM operations. The proposed model can be used to assist regulators and air navigation service providers for air traffic assignments in the strategic planning of UAM operations.

Résumé

Un grand nombre de véhicules de mobilité aérienne urbaine, en anglais Urban Air Mobility (UAM), devraient opérer dans l'espace aérien urbain dans un avenir proche. Cette croissance dépassera bientôt les capacités de l'espace aérien et des systèmes de gestion du trafic aérien actuels et mettra à rude épreuve les infrastructures de transport, entraînant des résultats indésirables tels que des encombrements, des situations de trafic complexes, des retards, etc. Par conséquent, des solutions UAM innovantes sont nécessaires pour garantir un transport urbain sûr et efficace. Dans cette thèse, nous nous concentrons sur la planification stratégique du trafic des drones, y compris l'organisation ainsi que la gestion de l'espace aérien et des flux de trafic. Plus précisément, nous développons des approches pour la conception du réseau de routes UAM et l'affectation du trafic aérien associé.

Afin de minimiser l'impact sur les systèmes de gestion du trafic aérien existants, cette thèse propose une méthodologie pour concevoir un réseau de routes UAM dans l'espace aérien urbain à basse altitude en présence d'obstacles et d'espaces aériens dangereux. Sur la base de diverses données publiques, cette approche a été appliquée à l'espace aérien urbain de Singapour comme cas d'étude pour les services de livraison de colis utilisant des UAVs. Dans un premier temps, le réseau de routes UAM est conçu à partir de grilles qui évitent les obstacles et les zones interdites. Ensuite, les coûts d'arcs sont définis en termes d'impact sonore sur les populations, de sécurité de l'espace aérien et d'efficacité du vol. Nous formulons un problème de k-plus courts chemins avec diversité pour sélectionner les routes réalisables qui minimisent le coût associé. Les itinéraires réalisables avec une faible similarité fourniront plus d'options de trajet, qui peuvent être utilisées pour générer le réseau de routes UAM afin de prendre en charge les opérations à haute densité et à flux complexe. En outre, l'impact de différents paramètres pour les coûts d'arcs sur les services UAM est également analysé. Bien que la méthodologie proposée soit appliquée à l'espace aérien urbain de Singapour, elle pourrait facilement être généralisée à d'autres espaces aériens urbains dans le monde.

Pour adapter la demande croissante à la capacité actuelle de l'espace aérien, nous avons proposé plusieurs modèles d'affectation du trafic aérien, notamment deux modèles statiques pour les réseaux de routes UAM à une et à plusieurs couches et un modèle dynamique pour les réseaux de routes UAM birectionnels multicouches, afin d'atténuer la congestion et la complexité et d'organiser la structure du flux de trafic aérien. Premièrement, au niveau macroscopique, les opérations UAM sont modélisées sous forme de flux de trafic en agrégeant la dynamique des véhicules individuels. Les flux sont répartis sur des réseaux de routes UAM, pour lesquels les routes aériennes sont modélisées comme des couloirs ou des segments volumétriques. Ensuite, le modèle d'affectation du trafic aérien est formulé comme un problème d'optimisation. Les critères sont modélisées de manière à décrire la complexité et la congestion du trafic aérien sur la base d'un système iii iv dynamique linéaire. Un cadre basé sur la simulation, comprenant différentes approches d'optimisation, est proposé pour résoudre efficacement ces problèmes. Des expériences informatiques sont réalisées sur des études de cas de réseaux de routes UAM à différentes échelles. Une comparaison avec les algorithmes conventionnels d'affectation du trafic est présentée. Les résultats montrent que l'approche proposée est capable d'affecter les flux de manière efficace et efficiente, en réduisant de manière significative la complexité du réseau de routes UAM 3D. Le modèle proposé peut être utilisé pour aider les régulateurs et les fournisseurs de service de navigation aérienne pour l'affectation du trafic aérien dans la planification stratégique des opérations UAM. morel, 2000). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.12 Eight converging/diverging cases. . . . . . . . . . . . . . . . . . . . . . . . 2. [START_REF] Bruynooghe | Un modèle intègre de distribution et d'afféctation du trafic sur un reseau[END_REF] Complexity heatmap in terms of convergence metric [START_REF] Delahaye | Metropolis -Development & Metrics Definition[END_REF]. 2.14 An iteration example of Dafermos' algorithm on a simple network. . . . . 2. [START_REF]Airspace classification information of Singapore[END_REF] The cooling procedures in terms of different temperature-reducing strategies. For a sudden temperature decrease (hardening), the material will move to a metastable solid state with non-minimal energy (top right).

The distribution of atoms in this state is not symmetric. For a slow decrease in temperature, the material can be cooled sufficiently long and reach a crystal solid state (bottom right). The atoms in this case are organized and symmetric. . . . . . . . . . . . . . . . . . . . . . . . . . . . 2. [START_REF] Caggiani | A metaheuristic dynamic traffic assignment model for OD matrix estimation using aggregate data[END_REF] The procedure of each transition for a given temperature in the SA.

Firstly, a neighboring solution is generated based on the current solution. After objective function evaluation, the acceptance criterion decides whether the new solution is accepted. If accepted, the current solution will be updated. Otherwise, the comeback operator returns to the previous solution and avoids duplication in memory. . . . . . . . . . . . . . . . . .

Urban air mobility

The urbanization process gradually shifts the population from rural to urban areas. Apart from high-volume passenger transport with rail, light rail, or metro lines, a large share of passenger urban and inter-urban transport today relies on an extensive road network in combination with cars or bus services and is often subject to significant congestion at peak periods. In fact, the congestion of ground transportation systems has already emerged in nearly all cities of different sizes since the 1980s. According to a report published by INRIX (2021), in the studied countries with the largest delay in 2021, people in the UK lost an average of 73 hours of congestion during peak commute periods per year, costing the country 8 billion GBP, an average of 595 GBP per person. The root cause is straightforward: the capacity of the existing transportation system did not keep pace with the growth in travel demand (Systematics, 2004). This situation gets worse as the population booms in metropolitan 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 a strain on transportation infrastructure, resulting in undesirable outcomes such as traffic congestion, complex traffic situations, delays, etc. [START_REF] Bauranov | Designing airspace for urban air mobility: A review of concepts and approaches[END_REF].

As a result, the requirement for alternative transportation modes in urban transportation grows. Innovative mobility solutions are required to ensure safe and efficient urban transportation. Urban Air Mobility (UAM) is one of the most potent and effective concepts for achieving high levels of automation in metropolitan areas. In the 1910s, pioneers started to conceive and implement a prototype called the "flying car" [START_REF] Cecil R Roseberry | Pioneer of Flight[END_REF], which is considered the first concept for roadable aircraft. Commercial UAM services can date back to the 1940s. At that time, Los Angeles Airways utilized helicopters to transport mail and passengers in the Los Angeles basin (Harrison, 2017), and New York Airways also used helicopters to transport passengers between Manhattan and the New York area. Due to several accidents, these operations with scheduled helicopter services have ceased after decades. Some other inventors also have made similar attempts in the same period [START_REF] Chabria | Meet George Jetson[END_REF]Patches, 2015). However, with low demand, limited technologies, and high costs, they are not successfully prompted. Until the 2010s, with the development of electric propulsion, enhanced battery capacity, sensing, automation, and miniaturization of UAM vehicles, UAM gradually started to resuscitate and enter the contemporary era. According to Cohen, Shaheen, and Farrar (2021), modern UAM expects a safe, sustainable, affordable, and accessible air transportation system for tasks within or traversing metropolitan airspace, including aerial surveillance, medical evacuations, rescue, filmmaking, news gathering, disaster relief, inspection, scientific research, ground traffic assessment, and weather monitoring (Thipphavong et al., 2018; [START_REF] Hall | The Use of Drones in the Spatial Social Sciences[END_REF]. These UAM operations can be executed with a variety of aircraft with different specifications, such as proportion type, design, technology, capacity, range, and autonomy [START_REF] David P Thipphavong | Urban air mobility airspace integration concepts and considerations[END_REF].

UAM provides on-demand or scheduled operations using manned and Unmanned Aircraft Systems (UAS) for intra-or inter-urban transport [START_REF] David P Thipphavong | Urban air mobility airspace integration concepts and considerations[END_REF]. Due to the development of new technologies, operational needs, and expected benefits in recent decades, there has been renewed interest in the domain of UAM. European union Aviation Safety Agency (EASA) claimed that by 2024-2025, UAM might be a lived reality in Europe (EASA, 2021). [START_REF] Anand | A Scenario-Based Evaluation of Global Urban Air Mobility Demand[END_REF] predicted the annual UAM passenger demands in terms of passenger kilometers, passenger and vehicle trips, and utilization for 542 top global cities considering the low and high scenarios during the years 2035 to 2050. The results over a five-year interval shown in Table 1.1 indicate that there will be a dramatic increase in UAM in the near future.

Unmanned aircraft system traffic management

Specifically, Unmanned Aerial Vehicles (UAVs) play an important role in UAM operations in low-level urban airspace. Also referred to as "drones", UAVs are aircraft without any human on board that are able to fly autonomously or be controlled remotely from 1.2. Unmanned aircraft system traffic management 3 the ground [START_REF] Chamayou | Théorie du drone[END_REF]. A large part of UAM services such as passenger transportation, parcel delivery, and traffic monitoring are expected to be provided by UAVs. There are several advantages of providing flying services by UAVs, including:

1. UAVs provide enhanced connectivity to remote regions. UAVs can make hard-toreach locations no longer out of reach. It can overcome the limitations of traditional methods for human hands-on operation, especially in dangerous situations; International Civil Aviation Organization (ICAO) projected that at the current pace, UAV operations would exceed the number of manned aircraft operations in the near future (ICAO, 2020). According to the newly published Federal Aviation Administration (FAA) aerospace forecast (FAA, 2020a), the total recreation and commercial UAVs already exceed 2 million by 2021. As detailed in Figure 1.1, the projected UAV fleets in three scenarios (low, base, high) are depicted from 2022 to 2026. Despite the two-year impact of COVID-19, UAVs are expected to grow smoothly and healthily in the following years all over the world. If each commercial UAV performs multiple missions per day, the number of daily operations will soon exceed millions (FAA, 2020b).

As of now, most UAM operations, especially remotely piloted UAVs, are required to be segregated in conventional airspace [START_REF] Straubinger | An overview of current research and developments in urban air mobility-Setting the scene for UAM introduction[END_REF]. In turn, the increase in UAM demand and UAM fleet size will result in a significant increase in terms of the volume of operations, control workload, increased traffic complexity, and even safety and security issues such as potential conflicts, which will be beyond the capacity and capability of the existing Air Traffic Management (ATM) system. In order to realize more efficient, scalable, safe, secure, and sustainable UAM operations with ever-growing demands, Unmanned aircraft system Traffic Management (UTM) is proposed as a solution to inter-operate and be consistent with existing ATM systems, with minimum impact on current airspace operations (FAA, 2020b). UTM is a specific subset of ATM, which is separated from ATM but also complementary to the ATM system. UTM aims to manage UAS operations efficiently, safely, and economically through a highly integrated and automated traffic management system in collaboration with airborne and groundbased parties (Thipphavong et al., 2018). In the UTM framework, the UAS operations are managed, organized, and coordinated by actors (operators and stakeholders) through highly automated traffic management systems. Air Navigation Service Provider (ANSP) still maintains the regulatory and operational authorities in terms of traffic and airspace operations. According to the UTM concept of operations (FAA, 2020b), the benefits of the UTM can be listed as follows:

Introduction

• A novel strategy for satisfying service requirements that makes use of commercial services to satisfy consumer needs while putting less pressure on the government in terms of infrastructure and manpower;

• A secure and stable environment where operators can share situational awareness and an operational framework made up of standards, regulations, and common procedures to satisfy business needs while lowering risk and maintaining stability;

• A scalable and adaptable framework that can adjust and evolve as the trade space develops and matures;

• A system that allows industry to oversee activities inside of designated low-altitude UAS flight areas and permits the FAA to maintain its authority over the airspace.

The first UTM conceptual framework was introduced by the National Aeronautics and Space Administration (NASA) (Kopardekar, 2014). Since then, together with FAA, NASA has continued to lead projects to support the development of UTM [START_REF] Prevot | UAS traffic management (UTM) concept of operations to safely enable low altitude flight operations[END_REF][START_REF] Ronald | Unmanned Aircraft System Traffic Management (UTM) Project[END_REF][START_REF] Raju | Making a UTM Ecosystem a Reality[END_REF][START_REF] Benjamin | US Airport Infrastructure Funding and Financing: Issues and Policy Options Pursuant to Section 122 of the 2018[END_REF]. In Europe, Single European Sky ATM Research (SESAR) introduced U-space in 2017 as an enabling framework to make large numbers of UAVs accessible to operate in a safe, efficient, and secure way in all types of airspace, route mission, and environment (Undertaking, 2017). Comprehensive implementation details are proposed in the U-space blueprint, including U-space foundation services, initial services for UAV operations management, advanced services for more complex operations in dense areas, and full services. Furthermore, ICAO created an overarching framework as a global baseline to avoid a lack of harmonization that could have an effect on safety, security, system reliability, environment, etc. (ICAO, 2020).

Similar to the ATM, in the operational aspect, there are some main challenges in the UTM, which are listed as follows from the most preventive to the most reactive level: (Mueller, Kopardekar, and Goodrich, 2017):

• Airspace Organization and Management (AOM);

• Traffic Flow Management (TFM);

• Separation provision;

• Collision avoidance.

The first two elements relate to strategic planning processes, with a relatively longer look-ahead time. The last two elements can be considered as tactical planning, in which measures take effect in a very short time. It is well established from a variety of studies that tactical planning processes for the UTM and decentralized UAM operations are only effective at low traffic volumes and become inadequate as traffic volume increases [START_REF] Holland | Optimizing the next generation collision avoidance system for safe, suitable, and acceptable operational performance[END_REF][START_REF] Sedov | Decentralized selfpropagating ground delay for UTM: Capitalizing on domino effect[END_REF][START_REF] Sedov | Centralized and distributed UTM in layered airspace[END_REF]. In addition, strategic planning can reduce the pressure on tactical planning and improve predictability. Strategic air traffic planning is required to build smooth and efficient air transportation systems by reducing airspace complexity and structuring air traffic. Though potential interactions or conflicts between aircraft may still exist, the separation of operations can be assured through various tactical conflict Introduction resolution methods [START_REF] Wing | New Flight Rules to Enable the Era of Aerial Mobility in the National Airspace System[END_REF]. Moreover, these remaining conflicts and interactions will be much easier to solve when the airspace complexity is reduced, and the traffic flow becomes more organized [START_REF] Bruno | A model for aircraft evaluation to support strategic decisions[END_REF].

Most current industrial applications on the UAM enables aircraft to select their preferred routes and maintain safety with tactical collision avoidance system, such as sense-and-avoid [START_REF] Yu | Sense and avoid technologies with applications to unmanned aircraft systems: Review and prospects[END_REF]. However, in future evolution stages, it has been suggested that the centralized UTM system and strategic planning have great advantages with high-density and high-complexity UAM operations in densely populated metropolitan areas (FAA, 2020a; [START_REF] Prevot | UAS traffic management (UTM) concept of operations to safely enable low altitude flight operations[END_REF][START_REF] Kenneth | Exploratory analysis of the airspace throughput and sensitivities of an urban air mobility system[END_REF].

Therefore, this thesis mainly focuses on strategic planning in the centralized control scheme of the UTM, including the AOM and TFM.

Airspace organization and management

The AOM adapts the capacity to the ever-increasing demand by designing, reconfiguring, and extending the infrastructure of airspace [START_REF] Barnhart | Demand and capacity management in air transportation[END_REF]. A variety of preliminary urban airspace design projects were conducted to integrate UAM with the current ATM system, making the urban airspace system safe, efficient, and predictable, with minimum impact on the existing airspace operations (Singireddy and Daim, 2018; FAA, 2020a; [START_REF] Schneider | Metropolis -Urban Airspace Design[END_REF]Jang et al., 2017a;[START_REF] Lascara | Urban Air Mobility Airspace Integration Concepts[END_REF][START_REF] Geister | Concept for Urban Airspace Integration DLR U-Space Blueprint[END_REF][START_REF] Sunil | Metropolis: Relating airspace structure and capacity for extreme traffic densities[END_REF]. Some representative projects including TU Delft [START_REF] Schneider | Metropolis -Urban Airspace Design[END_REF], FAA (FAA, 2020a), and Amazon (Singireddy and Daim, 2018), envisioned the UAM in the low-altitude urban area in class G airspace, which is uncontrolled below 1200 ft and not equipped with any ATM services [START_REF] Bauranov | Designing airspace for urban air mobility: A review of concepts and approaches[END_REF]. The common grounds in these projects are that urban airspace is divided into several airspace classes and layers. Furthermore, UAVs must operate in authorized airspace classes and fly below a certain altitude to keep away from planes, helicopters, and other manned vehicles. In addition, special route segment structures referred to as UAM corridors and volume segments were developed to integrate UAM operations in urban airspace with a longstanding ATM system. These route segments form a multi-level air transport network with high throughput, efficiency, and safety. In order 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).

Furthermore, to categorize the phases during the evolution of the UAM transportation system from the current state to a highly developed future state, NASA has developed a framework for UAM Maturity Levels (UML) (Goodrich and Theodore, 2021). As specified in Figure 1.2, during the intermediate and mature state of UML, the traffic density of UAM vehicles becomes intense. Especially in the mature state, the traffic density of UAM operations is characterized to have high density and complexity. A highly-integrated and automated route network with a centralized UAM transportation system has been claimed to have great advantages to these challenging scenarios that will arise in the near future (FAA, 2020a; Goodrich and Barmore, 2018).

Significant progress has been made in designing UAM route networks in urban airspace. FAA suggested adopting existing helicopter routes to construct route networks for early-stage UAM operations (FAA, 2020a). Tang et al. (2021) generated a 3D route network for all pairs of vertiports by using the visibility graph to avoid obstructions [START_REF] Tang | Automated flight planning of high-density urban air mobility[END_REF]. Some studies focus on UAM hub-and-spoke network design (Wu and Zhang, 2021; Willey and Salmon, 2021). Hong, Kuby, and Murray (2018) designed a delivery network to support commercial stand-alone drone delivery services in an area with obstacles. However, many studies have only focused on on-demand transportation services. Since on-demand services are often represented by a fully-connected graph (Willey and Salmon, 2021), the network size has to be limited in some studies. In addition, very few studies consider the realistic impact of network design. The lack of integration with real-world data is also a shortcoming of some studies.

Traffic flow management

Although the AOM is proven to support efficient UTM operation (ICAO, 2005), the airspace can usually be costly to be extended on a wider scale on the basis of existing infrastructure (Chaimatanan, Delahaye, and Mongeau, 2014). Another de-congestion strategy is to adapt the demand to the current capacity, which aims to make use of the airspace more efficiently by the TFM. Some research on the TFM has been undertaken. [START_REF] Balakrishnan | A distributed framework for traffic flow management in the presence of unmanned aircraft[END_REF] focused on large-scale TFM problems on predicted demand for the year 2030 with mixed UAS operations and passenger flights. Given airport and airspace capacities as constraints, the problem identifies a 4D trajectory for Introduction each aircraft that maximizes system-wide benefits including revenue and cancellation penalty, and minimizes the costs including operating costs and delay costs. [START_REF] Sedov | Centralized and distributed UTM in layered airspace[END_REF] It is noteworthy that the vast majority of trajectory planning studies have focused on individual vehicles. However, with increasing UAM traffic volume and high demand, few previous studies have performed the TFM from a macroscopic perspective, i.e., the traffic flows involving streams of air vehicles. Given the number of vertiports and timedependent demands between these vertiports, UAM flights are expected to operate as flows from origins to destinations following air routes in low-altitude urban airspace.

As a strategic and macroscopic model, air traffic assignment involves the optimal allocation of flows among alternative routes in transportation systems. The first air traffic assignment model was developed as early in 1954 [START_REF] Allen R Ferguson | The problem of routing aircraft: A mathematical solution[END_REF]. This model addressed the assignment for a given fleet to carry an anticipated traffic load over several routes at the minimum cost. Since then, a variety of traffic assignment models have been developed for solving routes and slot allocation problems [START_REF] Delahaye | Airspace congestion smoothing by stochastic optimization[END_REF] Considerable research has been conducted into solving road traffic congestion to minimize a user or system criterion by distributing demand over different routes. A common approach in the literature is to propose link cost functions, for which the parameters are derived from traffic and road conditions. Another type of cost function is developed on the basis of queuing theory [START_REF] Patriksson | The traffic assignment problem: models and methods[END_REF]. Moreover, a vast majority of research on ground transportation systems assumes that the link cost depends directly on the flows of this link and indirectly on other links in the surrounding area [START_REF] Robert | A comparative description of a capacity-restrained traffic assignment[END_REF][START_REF] Mosher | A capacity-restraint algorithm for assigning flow to a transport network[END_REF][START_REF] Spiess | Technical Note-Conical Volume-Delay Functions[END_REF][START_REF] Leta | Bottleneck and Queuing Analysis[END_REF][START_REF] Neuhold | Volume Delay Functions Based on Stochastic Capacity[END_REF]. Some other factors are also found to be influential in the route assignment process, such

Contributions

as distance [START_REF] Wang | Path-constrained traffic assignment: A trip chain analysis under range anxiety[END_REF][START_REF] Ryu | Two-Stage Bicycle Traffic Assignment Model[END_REF], emission (Patil, 2015; Tidswell and Raith, 2017), safety [START_REF] Ryu | Two-Stage Bicycle Traffic Assignment Model[END_REF][START_REF] Lin | Optimal transport network design for both traffic safety and risk equity considerations[END_REF] and road facilities [START_REF] Yang | Traffic assignment and signal control in saturated road networks[END_REF][START_REF] Kamel | Integrated simulation-based dynamic traffic and transit assignment model for large-scale network[END_REF].

However, to the best of our knowledge, these influential factors neglect the geometric characteristics of the network of urban airspace and consider the structured and disordered traffic patterns equally. Besides, the mixing level of aircraft trajectories also reflects inter-dependency between conflicts and the difficulty of managing the airspace [START_REF] Treimuth | Parallel complexity computation based on dynamical systems[END_REF]. In the literature, several metrics have been developed and applied to measure air traffic complexity, including geometric metrics [START_REF] Delahaye | Air traffic complexity: towards intrinsic metrics[END_REF], proximity metric [START_REF] Delahaye | Air traffic complexity: towards intrinsic metrics[END_REF], clusters metric [START_REF] Histon | Introducing Structural Considerations into Complexity Metrics[END_REF], Grassmannian metric (Delahaye and Puechmorel, 2013b), König metric [START_REF] Essén | Average angular velocity[END_REF][START_REF] Delahaye | A new air traffic complexity metric based on dynamical system modelization[END_REF][START_REF] Juntama | A Distributed Metaheuristic Approach for Complexity Reduction in Air Traffic for Strategic 4D Trajectory Optimization[END_REF]. Although most of these metrics are proven to be effective in basic traffic situations, they are not well-suited for large-scale or complex applications. Additionally, spatial-temporal information is not taken into account. To overcome these limitations, as an intrinsic complexity measurement, through analysis of the airspace geometry and traffic structure [START_REF] Delahaye | Air traffic complexity: towards intrinsic metrics[END_REF], Linear Dynamical System (LDS) is used in this thesis to evaluate air traffic complexity. A dynamical system describes the time dependence of an ensemble of particles in a geometrical space. The system behavior can be predicted for a short time into the future. The LDS has been successfully applied in many studies [START_REF] Delahaye | Air traffic complexity: towards intrinsic metrics[END_REF][START_REF] Delahaye | A new air traffic complexity metric based on dynamical system modelization[END_REF][START_REF] Delahaye | Air traffic complexity map based on non linear dynamical systems[END_REF][START_REF] Delahaye | Air traffic complexity based on dynamical systems[END_REF][START_REF] Treimuth | Parallel complexity computation based on dynamical systems[END_REF][START_REF] Delahaye | Air Traffic Complexity Map Based on Linear Dynamical Systems[END_REF]. With a concise mathematical formulation, due to its relative simplicity, efficiency, and mathematically predictable behavior, the LDS is quite suitable as a metric to estimate the local disorder and interaction of a set of trajectories in a traffic system.

Contributions

Based on the aforementioned discussions, this thesis contributes to the AOM and the TFM in the framework of the future UTM paradigm. More specifically, as shown in Figure 1.3, with increasing UAM traffic demand, several challenges such as air traffic congestion, high airspace complexity, and environmental issues are likely to impact UAM operations in the near future. UAM network design of the AOM and air traffic assignment of the TFM are strategic planning techniques to support high-density and high-complexity flow-based UAM operations.

The main contributions of this thesis to the AOM include: generated between OD pairs that minimize the noise impact on the population and maximize flight safety and efficiency;

Introduction

3. The feasibility of this approach is demonstrated by its application to a parcel delivery service for Singapore's urban airspace, and it can be extended to other urban airspaces as well;

4. The generated UAM route network has the potential to support high-density and complex flow-based UAM operations.

The main contributions of this thesis to the TFM include:

1. Several macroscopic formulations of static and dynamic air traffic assignment problems are introduced to handle future high-density UAM operations in UAM route networks with different configurations;

Thesis outline 11

The pre-defined UAM route network avoids the heavy computation of searching shortest paths and extra cost function evaluation for large-scale route networks;

5. The proposed air traffic assignment models are validated through scenarios of parcel delivery service with high-density UAM traffic in UAM route networks with different scales in Singapore's urban airspace;

6. The proposed framework can assist or provide advisories to UTM authorities and ANSP for various issues in the UTM, including the TFM, UAM traffic flow analysis, Urban airspace complexity measurement, and UAM route network evaluation.

Thesis outline

The thesis is organized as follows. Chapter 2 introduces the state-of-the-art in the context of UAM. This chapter consists of the development and trends of UAM airspace design, the literature review of air traffic complexity metrics, optimization methods and applications of air traffic assignment. In Chapter 3, we describe a framework to design the UAM route network in low-altitude urban airspace. Singapore's urban airspace is used as a case study. Chapter 4 and ?? address the problem of static air traffic assignment for high-density UAM operations in 2D and 3D UAM route networks, respectively. In Chapter 6, we deal with a dynamic air traffic assignment problem with time-varying demand. The UAM route network constructed by the methodology proposed in Chapter 3 is utilized as a case study. Finally, conclusions and perspectives are summarized in Chapter 7.

Chapter 2

State-of-the-art

Airspace design for UAM

As previously mentioned, the fundamental barrier to the growth of urban air transportation is the current ATM system's incapacity to manage urban airspace (Vascik and Hansman, 2017). In addition, several challenges including low-altitude operations, the varying performance of air vehicles, high-density and high-complexity UAM operations, further impede the integration of the UAM with the current ATM system (Vascik, Balakrishnan, and Hansman, 2018). Therefore, airspace design is urgently required for UAM operations.

To promote safety and facilitate the management of airspace, ICAO standardizes airspace classifications by designating airspace into categories A, B, C, D, E, F, and G, where classes A-E are controlled airspace and classes F and G are referred to as uncontrolled airspace (ICAO, 2001). Each airspace class contains a set of restrictions, requirements, and air traffic services provided in the associated airspace. Each nation is allowed to alter the requirements or restrictions for airspace classes. According to the literature (Cho and Yoon, 2018a), several terms that describe the availability of airspace are introduced. The raw availability of urban airspace is represented by free airspace, which is the airspace free of static obstacles. Free airspace can be further divided into unusable airspace and usable airspace. Unusable airspace includes prohibited areas, protected areas, and airspace affected by geofence associated with obstacles. The rest of the free airspace is denoted as usable airspace.

As stated by Bauranov and Rakas (2021), airspace-based operations are used to manage almost all flight operations in controlled airspace, where separations are managed and trajectories are assigned by Air Traffic Controllers (ATCos) within each sector. These processes are transferred from sector to sector, which cannot be applied to UAM operations. An idea for integrating the UAM with the existing ATM system is to increase the airspace capacity and require ATCos to manage all operations within the responsible airspace class. However, the current ATM system would have to be redesigned to accommodate the changes and the workload of ATCos would increase significantly, resulting in higher costs and difficulties in implementation. Therefore, it is more feasible to design new segregated airspace with specific rules and standards for future UAM operations (Jang et al., 2017b). The safety and security of high-density UAM operations can be improved by well-designing new urban airspace structures, which are expected to minimize the complexity and maximize the capacity [START_REF] Sridhar | Airspace complexity and its application in air traffic management[END_REF]. A large and growing body of urban airspace design initiatives and projects are proposed by industry, research institutions, and the government. Interested readers may refer to [START_REF] Bauranov | Designing airspace for urban air mobility: A review of concepts and approaches[END_REF] for a comprehensive review of urban airspace design. Several representative studies are summarized as follows:

FAA developed a concept of operations for UAM to describe the operational environment in the context of ATM and UTM (FAA, 2020a; FAA, 2020b), including airspace, types of operations, regulations, and procedures to support such operations. The urban airspace structure is shown in Figure 2.2. On the basis of UAM ConOps 1.0, the UAM corridor is defined to support the point-to-point operation of UAM aircraft without tactical ATC separation services [START_REF] Johnson | UAM Coordination and Assessment Team (UCAT)[END_REF]. It is a performance-based airspace of defined dimensions in which UAM aircraft follows specific rules, procedures, and performance requirements. UAM corridors can be seen as a primary mechanism to integrate safe and efficient UAM operations in urban airspace. The impact on existing ATM and UTM operations can be reduced. Besides, the needs of public interest stakeholder (e.g., environmental factors, congestion, safety, security) and stakeholder utility (e.g., customer need) can also be addressed. Figure 2.3 illustrates some planar UAM corridors of different flight levels. Further in UAM ConOps 2.0 (FAA, 2020b), Operation Volume Segments (OVS) concept is introduced to model the intended flight path. They are defined as 4D blocks of airspace. Each OVS has specified dynamic information including entry and exit times for the operator's UAM aircraft. An example of the OVS is illustrated in Figure 2.4. With higher navigational performance, UAM vehicles are allowed to operate in smaller volumes [START_REF] Decker | Economic policy choices and tradeoffs for Unmanned aircraft systems Traffic Management (UTM): Insights from Europe and the United States[END_REF]. The UTM environment is defined at or below 400ft above ground level. Other airspace operations occur in the ATM environment. Amazon proposed a layered airspace structure below 500ft. (Air, 2015a; Air, 2015b; A new airspace structure shown in Figure 2.6 has been proposed by SESAR Joint Undertaking (2020) in the concept of the Single European Airspace System. With the exponential growth in the use of drones across Europe, a new UTM system called Uspace has been put in place. The objective is to enable simultaneous drone operations and integration of all airspace users in all types of airspace, especially in urban areas. By Design, U-space is set to be scalable, highly autonomous, and connectible in the combination of emerging technologies. The first of four phases of U-space is currently underway (Undertaking, 2018).

To manage the increasing number of unmanned and manned traffic in urban airspace, 4 potential airspace concepts in the future were introduced in Hoekstra et al. (2015b), as shown in Figure 2.7. In these concepts, tactical CD&R is performed on each individual aircraft.

• Full mix: the traffic is unstructured and only subjected to physical constraints (weather, static obstacles, terrain, etc.). The aircraft are allowed to choose their optimal velocities, altitudes, and trajectories;

• Layers: the airspace is segmented into vertically stacked bands, in which each

State-of-the-art

altitude layer limits horizontal travel to within an allowed heading range. This structure is designed to reduce the probability of conflicts by limiting the relative velocity between aircraft in the same layer. But traveling at the optimal altitude cannot be ensured and efficiency is then limited;

• Zones: rather than vertical segmentation, the horizontal layout of the city is taken into account in the structure design. The airspace is segmented into circular and radial zones. The altitude of the aircraft is selected flexibly;

• Tubes: To structure airspace maximally, bi-directional 4D tubes that provide fixed route structure are implemented in this concept. Tubes are at the same horizontal level and never intersect except at the nodes. Only one aircraft is contained in each tube within a timeslot, and the time-based separation is introduced to maintain safety. Each node is occupied for a time interval if an aircraft pass it. The common points of 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 [START_REF] Bauranov | Designing airspace for urban air mobility: A review of concepts and approaches[END_REF]. 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 multilevel 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 CD&R.

In future high-density UAM operational environments, the most natural way to describe the air traffic is probably the air transport network (Vascik, Hansman, and Dunn, 2018). A transport network represents a structure of a spatial network that describes the movement of vehicles or flow of some commodity [START_REF] Barthélemy | Spatial networks[END_REF]. The construction of the air transport network follows the following general 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 [START_REF] Pang | A concept of airspace configuration and operational rules for UAS in current airspace[END_REF].

2. A trade-off of numbers of layers. The fewer horizontal layers, the more constrained and complex the network, however, the more efficient the computation. 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 increase.

Any transport network is concerned with demand and capacity. These two terms are opposing because capacity describes the network's ability in response to demand. A transport network is designed to support passengers, cargo, and other operations. The node represents waypoints, the link represents the air route segment, and the path represents air routes connecting two nodes. A feasible path set between an OD pair consists of a sequence of paths connecting the OD without cycles. Demands are associated with every pair of OD, which give rise to a traffic pattern.

In the metropolitan area, the number of paths connecting each OD pair may be very large. To reduce the computational burden and potentially decrease the airspace complexity, a natural idea is to reduce the number of candidate paths between each OD pair by searching the k-Shortest Paths (KSP) (Dunn, Grover, and MacGregor, 1994; Scano, Huguet, and Ngueveu, 2015). However, in KSP problems, it is very likely that some of the shortest paths are highly similar. They share a large part of sub-paths, and the resulting paths provide low flexibility and few alternatives for air traffic assignment. To this end, we search k-Shortest Paths with Diversity (KSPD). Since KSPD has been proven to be NP-hard [START_REF] Liu | Finding top-k shortest paths with diversity[END_REF], some algorithms are developed to approximate the KSPD problem [START_REF] Liu | Finding top-k shortest paths with diversity[END_REF][START_REF] Hanaka | Computing diverse shortest paths efficiently: a theoretical and experimental study[END_REF][START_REF] Chondrogiannis | Finding k-shortest paths with limited overlap[END_REF][START_REF] Häcker | Most Diverse Near-Shortest Paths[END_REF]. In addition, various similarity functions have been proposed to measure the diversity between two paths in different perspectives [START_REF] Liu | Finding top-k shortest paths with diversity[END_REF].

Air traffic complexity metrics

Both in newly designed airspace and existing airspace, it is necessary to assess the impact of air traffics on control systems and their interactions.

To these ends, air traffic complexity is widely used in the ATM as it is a critical component of the airspace optimization process [START_REF] Lee | Air traffic complexity: an inputoutput approach[END_REF]. Precursory studies towards air traffic complexity emerged in the 1960s [START_REF] Davis | The influence of selected sector characteristics upon ARTCC controller activities[END_REF][START_REF] Atid | The control load and sector design[END_REF] and have continued to develop since then. According to Meckiff, Chone, and Nicolaon (1998), air traffic complexity consists of three elements: the geometric nature of air traffic, the operational procedures and practices used to handle the traffic, and the characteristics as well as the behavior of the traffic control system involving human operators or automatic processes. These elements can be further refined to control workload and traffic complexity. More concretely, control workload considers these elements and thus relies on a precise definition of operational procedures, for example, the procedures related to sectorization and control system behavior, while traffic complexity focuses on the geometric nature of traffic.

Control workload

In air traffic management, airspace is divided into virtual sections known as sectors. The tasks are assigned to each sector and are decomposed into workloads. Control workload is a measurement of the difficulty of the traffic control system in treating a traffic situation. The traffic control system may be a human operator or an automatic process. In the 20th century, pioneering models existed in the following literature regarding the measurement of control workload associated with the given traffic situation:

• Queue-based model [START_REF] Maugis | Techniques de détermination de la capacité des secteurs de contrôle de l'espace aérien: statistiques et simulations[END_REF]: the aim of this model is to evaluate a maximum acceptable arrival rate in a sector by using queueing theory.

The control sector is modeled as a system that is able to receive aircraft and provide services. Then, the sector can be modeled as a service center containing several servers and an aircraft queue.

However, this model has various well-known limitations. In many future concepts of air traffic management, such as dynamic sector configuration and advanced traffic flow management, the essential is to evaluate and predict the control workload.

For instance, the operational capacity of a control sector is evaluated by the maximum number of aircraft that cross the sector in a given period of time. This evaluation does not consider the aircraft directions, geometric structures, and disorder of traffic. Therefore, in some cases, when the traffic is structured, even if the operational capacity in terms of the total flow is exceeded, the controller is likely to continue receiving traffic. In other cases, when the traffic is disordered, even if the operational capacity in terms of the total flow is not reached, the controller may prohibit additional traffic. .

• Workload model based on the traffic level [START_REF] Kirchner | The human operator in air traffic control systems[END_REF]): The workload is defined as the proportion of control time per hour. The following indicators can be used to characterize the stress of ATCos (David K. Schmidt, 1978):

the average duration of routine control tasks for an aircraft;

the average time to resolve conflicts per aircraft;

the average rate of arrivals in a sector per hour;

the average rate of conflicts in a sector per hour.

It is commonly known that airspace with a low-level structure is likely to lead to more conflict scenarios. The severity of conflicts is derived from the comparison with the threshold of predicted time to the closest point of approach. The percentage of aircraft in conflict can be used to evaluate the number of conflicts.

• Model based on airspace structure [START_REF] Janić | En route sector capacity model[END_REF]: in this model, the workload of a sector mainly depends on its geometric structure:

- -Number of conflicts and distance between aircraft for different ranges Nevertheless, it is difficult to generalize these models of control workload to new sectors, as the models are highly dependent on the involved controllers. Consequently, the basic measurement is not sufficient to represent the difficulty of a certain traffic situation.

Several attempts have been made to improve the modeling of the workload in air traffic control in terms of operational aspects. Some important operational metrics that evaluate the operational differences between airspace structures are summarized (Hoekstra et al., 2015a):

• Operational metrics: they are defined in terms of the following four aspects:

-Safety: the safety metric concentrates on the capability of maintaining safe separation between aircraft, which can be evaluated by the number and severity of loss of separation and conflicts: * Loss of separation: the severity of loss of separation is computed as the maximum intrusions during the considered time period, where the intrusion is defined as the minimum value of the horizontal and vertical intrusions that are normalized with respect to the dimension of their protection zones. The average severity is calculated by the sum of severity averaged by the number of losses of separation. * Conflicts: as one of the important indicators of airspace structures, the conflict has been widely investigated in some pioneering control workload models [START_REF] David | A queueing analysis on the air traffic controller's workload[END_REF]. Further studies were made to measure the control workload by automatic conflict resolution algorithms [START_REF] Lee | Describing Airspace Complexity: Airspace Response to Disturbances[END_REF][START_REF] Lee | Air traffic complexity: an inputoutput approach[END_REF][START_REF] Granger | A traffic complexity approach through cluster analysis[END_REF]. In these researches, the control workload is represented by the number of required modifications on trajectories for conflict resolution.

-Stability: when a conflict is detected, aircraft are maneuvered to avoid conflict, which may generate new conflicts and lower the stability. To measure the stability of the airspace, a domino effect parameter is proposed. The idea of the domino effect parameter is to compare the number of conflicts with and without conflict resolution in the same traffic situation.

-Efficiency: the efficiency metric is designed to measure the utilization of airspace. Efficiency metrics can be categorized into four types: * Route efficiency: route efficiency is linked with fuel consumption. The route efficiency is measured by comparing the actual trajectory flown with the great circle distance; * Delay and delay absorption capability: when the airspace capacity is fully utilized, delays may be expected until the traffic congestion is alleviated.

The delay metric can be expressed as the sum and average of delays of all aircraft. Relative delay absorption is another efficiency metric designed to measure the self-adjusting capability of delay when facing different levels of traffic situations. The trip times from origin to destination can be used to reflect this metric; * Sequencing: the sequencing metric is adopted to evaluate the impact of airspace structure on arrival or departure management procedures. This metric measures the time interval between two successive arrivals of an airport or a droneport in a given period of time.

-Capacity: through the analysis, the capacity metric is related to the gradient of the safety and efficiency metrics with respect to demand. Another method to measure the capacity is to evaluate the extent to which traffic density matches the predefined traffic demand.

These metrics provide novel and valuable insights in terms of operational aspects. Nevertheless, they still depend on the conflict resolution algorithms involved. In addition, most metrics require a large amount of data over a long period of time, which is not suitable for evaluating the complexity of new airspace.

In order to overcome these issues, more intrinsic approaches emerged. One of the representative works is the flow-based model of control workload [START_REF] Delahaye | Airspace congestion smoothing by stochastic optimization[END_REF]. This model is a macroscopic model designed for en-route traffic. The traffic flow circulates on transportation networks, where nodes represent airports or main waypoints and links represent airways. In this work, the control workload is defined as the sum of the following quantitative criteria:

• Conflict workload: the conflict resolution workload is related to the crossing of flows at nodes in the network. It is defined as the sum of workloads at each node within a sector;

• Coordination workload: this workload comes from the transfer of control. The transfer imposes a heavy burden on controllers in charge of different sectors. The burden may be derived from the change in frequency, misunderstanding, and errors that sometimes occur. For instance, in a sectorized transportation network shown in Figure 2.8, the coordination workload is linked to the flow cuts across boundaries of sectors, which include two cases: only one extremity of the link belongs to the sector as [START_REF] Ameli | Cross-comparison of convergence algorithms to solve trip-based dynamic traffic assignment problems[END_REF][START_REF] Ameli | Simulation-based dynamic traffic assignment: Meta-heuristic solution methods with parallel computing[END_REF] and two extremities of the link are not in the sector as [START_REF] Prime | Determining safe access with a best-equipped, best-served model for small unmanned aircraft systems[END_REF][START_REF] Ameli | Simulation-based dynamic traffic assignment: Meta-heuristic solution methods with parallel computing[END_REF]. Note that in the case where the two extremities of the link belong to the sector, as [START_REF] Prime | Revising the airspace model for the safe integration of small unmanned aircraft systems[END_REF][START_REF] Ameli | Cross-comparison of convergence algorithms to solve trip-based dynamic traffic assignment problems[END_REF], no coordination workload is produced.

• Monitoring workload: monitoring in a control sector is conducted on aircraft, not in conflicts or involved in transfers by controllers. The workload consists of realtime checking of the features (position, speed) of aircraft in a sector. This workload can be modeled with the length of the link that intersects with the sector and the average flow speed.

Figure 2.8: An example of a sectorized air transport network. For sector S1, there are conflict workloads at node 2, coordination workloads at links 6-2, 1-2, 3-5, 3-4, and 1-4, and monitoring workloads at all links.

Traffic complexity

Traffic complexity is an intrinsic measurement of the complexity involved with a traffic situation. It is independent of the traffic control system and only dependent on the geometry of trajectories. As is known, maintaining the safe separation between air vehicles and obstacles is always a crucial issue in aviation, which is ensured by CD&R. The process of CD&R is to detect a potential conflict in the future, inform the operator, and in some cases, assist to resolve the conflict (Kuchar and [START_REF] James | A Review of Conflict Detection and Resolution Modeling Methods[END_REF]. When a conflict is detected, aircraft are maneuvered to avoid conflict, which may generate new conflicts. This interdependency between conflicts is dependent on the mixing level between trajectories. Figure 2.9 illustrates 3 traffic situations with the same number of trajectories but different levels of difficulty, predictability, and interdependency. In view of parallel trajectories, there is no conflict, and the sensitivity is low, which is an easy situation. As for concentrated trajectories, potential conflicts are likely to occur, but they are easy to be resolved by creating a roundabout with trajectories by giving all aircraft the same change of direction order (±90 • ). It is regarded as an average situation. The last case involves disordered trajectories. The potential conflicts exist with high-level interdependency between trajectories, which is a difficult situation. In airspace including a large number of air vehicles, interdependency appears much more often. This impact may destabilize the airspace by means of propagating large numbers of conflicts. Only a small amount of maneuvering space is remained, which imposes a heavy computational burden on CD&R algorithms. Besides, the initial conditions and uncertainties in terms of intents such as positions and speeds make it more difficult to conduct precise trajectory prediction, especially for UAVs, which are more flexible and less constrained. In summary, the difficulty of controlling airspace links to its sensitivity to initial conditions and the interdependency of conflicts.

In recent years, intrinsic traffic complexity metrics that are solely linked to trajectory structures have been investigated:

• Geometric metric [START_REF] Delahaye | Air traffic complexity: towards intrinsic metrics[END_REF]: these metrics are based on the real-time geometric information of aircraft in the airspace, such as the position and speed. The following geometric metrics exhibit the different characteristics of traffic complexity:

-Proximity metric: the proximity metric is determined by the positions of aircraft in the airspace. It is used to distinguish the distribution of aircraft and to detect dense traffic structures. An example is given in Figure 2.10. In -Clusters metric: clusters metric is designed to model multiple interactions between aircraft. It is constructed by combining aggregation and convergence levels based on relative distances and speeds;

-Grassmannian metric: it is proposed to calculate a pseudo-measurement of disorder by constructing the Grassmannian matrix of relative speed vectors between aircraft pairs, taking account of relative distances;

-König metric [START_REF] Juntama | A Distributed Metaheuristic Approach for Complexity Reduction in Air Traffic for Strategic 4D Trajectory Optimization[END_REF][START_REF] Delahaye | Air traffic complexity: towards intrinsic metrics[END_REF]: this metric is used for identifying the complexity related to translation or curl by application of König's theorem. The idea is to compute the dispersion of the kinetic moment with regard to the average kinetic moment located at the barycenter.

Geometrical metrics are very suitable for capturing different kinds of complex features. Gathering them together in some way can be much more effective in measuring airspace complexity. Nevertheless, when a traffic situation has a more complex structure, these geometric metrics would fail. In addition, the measurement is made for a given instant, which couldn't extend to a time period.

• Metrics based on dynamical systems (Delahaye et al., 2022; Delahaye and Puechmorel, 2000): Independent of airspace structure and ATC system, these metrics are intrinsic to identify organized and disordered traffic according to traffic dynamics. They are based on linear or nonlinear dynamical systems. They are utilized to quantify the level of disorder and interaction between trajectories in the airspace. These metrics can be applied to specific airspace with a series of aircraft trajectory points including position and speed information.

-Metric based on linear dynamical systems: air traffic complexity metrics based on the LDS has been successfully applied in many studies ( In this metric, an LDS is applied to model a set of trajectories in the airspace with the following equation:

ẋ = Ax + b (2.1)
where x is the state vector of the system or the position vector of aircraft, and ẋ is the speed vector. The eigenvalues of matrix A determine the mode of evolution of the dynamic system. The evolution of the LDS is illustrated in Figure A.1 in Appendix A.5. The vertical strip around the imaginary axis corresponds to organized traffic situations, where the eigenvalues of A e,i have a real part close to zero. In these situations, the relative distances between points change slowly with time, namely, their relative speeds are close to zero and there is nearly no interaction between them.

State-of-the-art

-Metric based on nonlinear dynamical systems: this measurement is designed to increase the accuracy when the number of observations is increased. Unlike linear dynamical systems that represent the aggregated behavior of the underlying traffic, nonlinear dynamical systems exhibit the local behavior of the traffic. The air traffic complexity is measured by Lyapunov's exponent theory, which quantifies the local level of organization of the vector field. However, it is much more computationally consuming. Especially, the regression step has a time complexity of O(N 3 ), where N is the number of aircraft (Delahaye and Puechmorel, 2013c).

Air traffic assignment

To mitigate air traffic complexity and congestion, air traffic assignment models, as strategic planning approaches, have been commonly studied. This study considers air traffic assignment as a general term that encompasses the allocation of trips on routes between OD in air transportation networks. In the following sections, The principles and optimization algorithms of traffic assignment will be explained. Then, its application to air traffic management and UAM will be discussed.

Traffic assignment principles

Traffic assignment problem refers to allocating estimated future OD trips to routes in the existing transportation networks or systems [START_REF] Wang | Dynamic traffic assignment: A review of the methodological advances for environmentally sustainable road transportation applications[END_REF]. Given the complete description of transportation systems and trip movements as inputs, traffic assignment aims to output the estimation of traffic volumes by minimizing transport costs, including travel time, travel distance, etc. Traffic assignment can also be carried out for the following purposes: [START_REF] Patriksson | The traffic assignment problem: models and methods[END_REF]:

• To identify the deficiencies in the existing transportation system by allocating estimated future demand to routes;

• To test the impact of improvement and extensions on the current transportation system through traffic assignment on the modified transportation network;

• To determine the construction priorities for the current transportation system through traffic assignment in a long range;

• To test alternative transportation systems;

• To design hourly traffic volumes and turning movements.

Depending on the characteristic of demand, the traffic assignment can be categorized into Static Traffic Assignment (STA) and Dynamic Traffic Assignment (DTA), where traffic demand is constant over time in the STA, and varies over time in the DTA. On the other hand, from the user's point of view, traffic assignment can also be classified into system optimal assignment and user optimal assignment, where users minimize the global travel cost in system optimal assignment, and users try to minimize their own travel cost in user optimal assignment. Furthermore, traffic assignment may also be defined for deterministic and stochastic cases. The detailed interpretations of these terms are given as follows.

• User optimal assignment: The principle of route selection is that each user selects the route with the minimum cost. When it is impossible for any user to reduce the transport cost by choosing alternative routes, equilibrium is reached. All routes of an OD pair with the flow have the same cost, while other routes not in use have a higher cost. This equilibrium is described by Wardrop's first principle (Wardrop, 1952):

Definition 2.1 (Wardrop's first principle)
The travel costs on all routes actually used are not greater than those which would be experienced on any unused routes.

Wardrop's equilibrium is a Nash equilibrium between network users. When this user equilibrium is reached, the routes currently used have the minimum transport cost. If it is assumed that each user has complete knowledge of the network as well as the estimation of costs in terms of different routes, deterministic user equilibrium is reached.

Another variant is the stochastic user equilibrium. This problem considers OD pairs of several user groups existing in a stochastic network. Each OD route is modeled by a set of links with random non-negative cost values. Each user perceives a transport cost on its route with its own random error linked to its class. The objective is to determine the stochastic equilibrium of the resulting flow in the network. Besides, it is supposed that each user behaves unilaterally depending on the traffic encountered on its route. The behavior can be modeled by a repugnance function. In the stochastic framework, it is possible to find a situation of stochastic user equilibrium that satisfies the generalization of Wardrop's first principle [START_REF] Mirchandani | Generalized traffic equilibrium with probabilistic travel times and perceptions[END_REF]:

Definition 2.

(Stochastic static generalization of Wardrop's first principle) The expected value of the repugnances perceived by users on current routes is not greater than those on any unused routes.

The stochastic formulation is more stable and insensitive to slight variations in the network compared to the deterministic one, but the resolution method requires much more computational costs. It is proved that for heavily congested networks,

State-of-the-art

stochastic and deterministic equilibrium are equivalent in practice. Nevertheless, in the context of air traffic, the network is generally not heavily congested, so the stochastic hypothesis is more suitable. Unfortunately, very few resolution algorithms are available for this case.

• System optimum assignment: The system equilibrium assumes that the network is non-autonomous, in which the users cooperate to minimize the overall transportation cost. It is expressed by Wardrop's second principle (Wardrop, 1952):

Definition 2.3 (Wardrop's second principle)
The average travel cost for all users is the minimum at system equilibrium.

This principle states that not all users but the average take the shortest route.

In the resulting flow, the marginal costs on the routes used to connect an OD pair are not greater than those of the unused routes between the same OD pair. System optimum assignment is not a behaviorally realistic model, but it still can be widely applied in practice. On the opposite, it is realistic in the situation where the authority in charge of traffic is able to control the traffic flows.

• Dynamic Assignment: To obtain a more precise adjustment of demand in terms of capacity and represent the microscopic relationship between them, the temporal evolution of traffic in the network needs to be taken into account. Dynamic assignment deals with traffic demand and link costs that vary with space and time. It consists of predicting the evolution of traffic in a congested network. The dynamic assignment can be classified into two classes:

-Optimal dynamic assignment from the user's perspective: Aims to minimize users' own transport costs.

-Optimal dynamic assignment from the system's perspective: Aims to minimize overall transport costs.

Dynamic assignment adapts in deterministic and stochastic cases (Delahaye and Puechmorel, 2013c):

-Deterministic dynamic assignment: The Deterministic Dynamic user equilibrium can be expressed as a generalization of Wardrop's first principle:

Definition 2.4 (Deterministic dynamic generalization of Wardrop's first principle) The disutilities (including schedule delays and the travel costs) on all routes actually used are not greater than those which would be experienced on any unused routes.

The objective of this model is to determine the departure time and route for each user in the network in order to reach the user equilibrium. The computational complexities of resolution algorithms are very high, especially when optimizing the departure time and routes at the same time. This problem is proved to be NP-hard and multimodal. The proposed methods only provide a partial solution to the problem.

-Stochastic dynamic assignment: The stochastic dynamic user equilibrium can be similarly expressed as a generalization of Wardrop's first principle in terms of disutilities perceived by users:

Definition 2.5 (Stochastic dynamic generalization of Wardrop's first principle) The expected value of disutilities (including schedule delays and the travel costs) perceived by users on all routes actually used is not greater than those which would be experienced on any unused routes.

This model seeks to optimize the flow on different routes to reach equilibrium.

To solve these traffic assignment problems, the basics of optimization and the related optimization methods are given in the following section.

Optimization methods for traffic assignment

Basics of optimization

Mathematical optimization consists of selecting the best decision variables in terms of some criterion from a feasible set of solutions that satisfy some constraints. Optimization can be considered as a decision-support tool to solve real problems. The first step involved in an optimization process is to model the real problem through an efficient mathematical abstraction. A mathematical model is constructed by characterizing state space, objective function, and constraints. These terms are detailed as follows:

• State space: The state space refers to the set of decision variables of the system in which one or more objectives can be optimized.

• Constraint: Constraints are the conditions that the solution must satisfy.

• Objective function: The objective function is the function to minimize or maximize that evaluates the solution. A feasible solution that minimizes or maximizes the objective function is referred to as the optimal solution. According to the number of objective functions that are optimized simultaneously, the optimization problem can be classified into a single-objective problem or multi-objective problem.

In addition, the selection of suitable resolution algorithms is influenced by many factors. For example, one of the important factors is the type of objective function, which can be distinguished as follows (Delahaye and Puechmorel, 2013a):

• Criterion with analytical form: It is the most computation-friendly case. Sometimes, it is even possible to derive the optimality conditions or the analytical form of the optimal solution. However, many real problems cannot be abstracted into a formulation that consists of criteria with an analytical form.

• Cheap criterion evaluated by numerical computation: For criterion without analytical forms such as black-box function, the numerical computation can be conducted to evaluate the criterion. Additional information needs to be provided to resolution algorithms. According to the resource required to evaluate the criterion, the cheap function can be evaluated many times, typically more than thousands of times.

• Costly criterion evaluated by complex simulation process: Some criteria need to be evaluated by a complex simulation process. These criteria can only be evaluated a few times, typically less than hundreds of times. As the evaluation is expensive, it is usually difficult to obtain additional information.

The global optimization algorithms can be classified into:

• Deterministic approaches Deterministic optimization algorithms conduct the exact computation and provide the global optimal solution based on the analytical properties of the problem.

It is ensured that the global optimal solution will be approximated within predefined tolerance after a finite time. Typical deterministic approaches include branch-and-bound (Lawler and Wood, 1966), generalized benders decomposition (Geoffrion, 1972), and outer approximation (Duran and Grossmann, 1986), which are designed for problems such as linear programming, nonlinear programming, and mixed-integer nonlinear programming. An important strength of most deterministic methods lies in the fast convergence speed. Only a few iterations involving the evaluation of the objective function are required to reach the optimal solution. However, deterministic optimization methods highly depend on the nature of the problem. For many real problems such as complex, high dimensional, and NP-hard problems, a large computation time is required to solve these problems, which is not feasible in practice (Garey and [START_REF] Michael | Computers and Intractability[END_REF].

• Heuristic approaches

Heuristic optimization algorithms are designed to be more flexible and efficient than deterministic approaches by sacrificing optimality. Heuristic approaches mainly search for a near-optimal solution by choosing the nearest neighbor in the solution space. The optimality is not guaranteed, but a feasible and good-quality solution can be provided in a reasonable time. Heuristic methods are useful for NP-hard problems or problems where the objective function has no analytical form. However, conventional heuristic methods are highly dependent on the problem to solve, which is difficult to generalize to a new scenario or problem. In addition, as the searching strategy is greedy, heuristic approaches can easily fail in local optima.

On the other hand, some higher-level heuristic approaches are developed to improve the performance of conventional heuristic approaches, which are referred to as metaheuristic approaches. They are more advanced in balancing exploration and exploitation by designing different strategies. Exploration, also referred to as diversification, aims to explore a bigger domain of state space to avoid getting trapped in local optima. Exploitation, also referred to as intensification, searches for a good-quality solution in the promising bounding region of state space. In addition, metaheuristics are problem-independent. They are usually inspired by natural phenomena and mimic the behavior of agents. There is various kind of metaheuristics, from the simplest ones such as local search to the complex ones including swarm intelligence algorithms.

Traffic assignment algorithms

The mathematical modeling of traffic assignment is usually performed by abstracting traffic networks through graphs. According to the components of a graph, traffic assignment problems can be formulated by two types of models, i.e., link-based models or path-based models [START_REF] Long | Link-based system optimum dynamic traffic assignment problems with environmental objectives[END_REF]. For the first type, the decision variables correspond to the link flows, and the constraints enforce flow conservation at each node. The second type computes a set of routes between each OD pair and outputs path flows. Moreover, path constraints may be included in this formulation. The benefit of this formulation is that it can provide accurate path-based flow allocation results, which can facilitate the decision-making process.

Next, we will introduce widely used optimization methods in the context of pathbased traffic assignment studies including STA or DTA to achieve system or user equilibrium.

• All-Or-Nothing (AON) assignment

The AON assignment is one of the simplest assignment approaches. This model assumes that all traffic between a specific OD pair will take the shortest path. No consideration is given to the link capacity as well as congestion effects. All users perceive the route choice in the same way. No user is assigned to other routes if there exists more than one path with the same or approximate cost. This model is generally unrealistic. However, it is reasonable to be applied in uncongested networks where alternative routes are limited and have very different costs. In practical use, the assignment produced by such an algorithm is considered as the initial solution or a component for more advanced assignment models.

• Capacity restraint assignment

Capacity restraint assignment handles the overloaded links in the network. It approximates an equilibrium solution by iteration between all-or-nothing traffic loadings and recalculating link costs based on a congestion function that reflects link capacity. A representative example is an incremental assignment. In addition to the hypothesis that each user will choose a path to minimize the link costs, the incremental assignment also assumes that the link cost varies with the flow. In the assignment process, fixed proportions of traffic volumes are assigned in steps based on an all-or-nothing assignment till all trips are assigned. After each step, the link cost is re-evaluated based on link volumes until all trips are assigned to the transportation network. However, this approach does not provide an equilibrium solution but is only close to user equilibrium when many increments are used. Also, in view of a large-scale network, the incremental assignment becomes unsuitable. Furthermore, the inconsistencies between link volumes and costs can lead to bias in results.

• Frank-Wolfe algorithm

The Frank-Wolfe algorithm is a first-order iterative algorithm designed for constraint convex optimization problems [START_REF] Frank | An algorithm for quadratic programming[END_REF]. The algorithm approximates the objective by linear approximation and conducts a line search to determine the step size by minimizing the objective over the defined line segment. A termination is expected when the step size is sufficiently small. The Frank-Wolfe algorithm was firstly applied to traffic assignment by [START_REF] Bruynooghe | Un modèle intègre de distribution et d'afféctation du trafic sur un reseau[END_REF]. Since then, variants of the Frank-Wolfe algorithm have been proposed to improve the original version (Daneva and Lindberg, 2003; Mitradjieva and Lindberg, 2013).

• Dafermos' Algorithm (DA)

As an efficient sequential decomposition algorithm, the DA has been used to efficiently solve traffic assignment problems. It is designed to determine the flow patterns under system equilibrium in transportation networks (Dafermos and [START_REF] Stella | The traffic assignment problem for a general network[END_REF][START_REF] Stella C Dafermos | An extended traffic assignment model with applications to two-way traffic[END_REF][START_REF] Stella C Dafermos | The Traffic Assignment Problem for Multiclass-User Transportation Networks[END_REF]. The DA requires successive application of the equilibrium operator to decompose the optimization problem into a sequence of simpler ones. This algorithm enumerates the set of OD pairs W. For each OD pair w ∈ W, the most expensive route p w,max and the least expensive route p w,min in terms of the performance indicators are identified in the set of paths P w . The marginal cost C ′ p (f ) is introduced as the performance indicator that measures the variation of the cost function C (f ) in response to a small change of the flow F p on a path p ∈ P w . It is defined as follows:

C ′ p (f ) = ∂C ∂F p (2.2)

Air traffic assignment 35

where f is the set of link flows.

Considering a cost function C (f ) defined as the sum of link cost functions C e (f ):

C (f ) = e∈E C e (f ) (2.3)
where E is the set of links in the transport network.

Taking into account the structure of the cost function, the marginal cost can be expressed as:

C ′ p (f ) = a,b∈E ∂C a (f ) ∂f b δ b,p (2.4)
where

δ b,p = 1 if link b is contained in path p, otherwise δ b,p = 0.
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 (Dafermos, 1971). In the proposition, F i w is the flow on i-th path of OD pair w, R w is the cardinality of the path set between OD pair w.

Proposition 2.1 F is a system equilibrium flow pattern if and only if for any OD pair w ∈ W with feasible routes p 1

w , • • • , p Rw w , the marginal cost of these paths can be sorted to satisfy:

C ′ p 1 w (f ) = • • • = C ′ p m w (f ) ≤ C ′ p m+1 w (f ) ≤ • • • ≤ C ′ p Rw w (f ) (2.5) With F i w > 0, i = 1, • • • , m (2.6) F i w = 0, i = m + 1, • • • , R w (2.7)
Since the marginal costs for all paths with the non-zero flow are equal when the system equilibrium is reached, this algorithm searches the most expensive route p w,max and the least expensive route p w,min in terms of marginal cost for each OD pair. Then, a quantity of flow is transferred from p w,max to p w,min such that the global cost is minimized. This procedure is repeated for all OD pairs iteratively until the global cost cannot be further reduced. Details of the process are given below.

To balance the marginal cost, a certain quantity of flow should be transferred from the most expensive route p w,max to the least expensive route p w,min such that the

State-of-the-art

total complexity is minimized, where

   p w,min = arg min p∈Pw C ′ p (f ) p w,max = arg max p∈Pw,Fp>0 C ′ p (f ) (2.8)
The flow transfer operation is defined as an equilibrium operator E w :

E w : R Rw → R Rw (2.9) F w → E w F w (2.10)
where F w is the path flow pattern between OD pair w.

Starting from an initial flow pattern

F (0)
w , the cost function converges through iterating E w :

F (n) w = E w F (n-1) w , n ∈ N * (2.11)
New path flow pattern Fw is updated using the following rule:

       Fp = F p , ∀p ∈ (P w -{p w,min , p w,max }) Fp w,min = F p w,min + σ w Fpw,max = F pw,max -σ w (2.12)
where σ w is the solution to the following optimization problem: (2.13)

The element σ 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.

The process of an iteration in Dafermos' algorithm is illustrated in Figure 2.14 in terms of a simple network with several links connecting two nodes v i and v f .

• Method of Successive Averages (MSA)

The MSA [START_REF] Robbins | A stochastic approximation method[END_REF]) is widely used in traffic assignment problems (Zhang, Liu, and Waller, 2019; Tran et al., 2021). It is easy to implement and the derivative of the cost function is not required.

In the basic formulation of the MSA, the algorithm is initialized by a feasible flow pattern. At iteration n, n ≥ 0, the path flow pattern of current OD w is denoted as

F (n) w . The cost of all paths in F (n)
w is then computed. The path with the minimum cost is assigned with the demand D w , which is referred to as the All-Or-Nothing (AON) assignment. The resulting path flow is denoted as F(n) w . The path flow pattern is updated by the following strategy:

F (n+1) w = F (n) w + η n ( F(n) w -F (n) w ) (2.14)
where η n is the n-th step size. It satisfies the following constraints:

n η n = ∞ (2.15) n η n 2 ≤ ∞ (2.16)
A typical way is to predetermine the step size as follows:

η n = 1 n + 1 (2.17)
Though MSA is efficient on small networks, it has drawbacks in determining step size and flow-swapping strategy.

• Probabilistic Method (PM) Ameli, Lebacque, and Leclercq (2020a) made a cross-comparison between conventional optimization algorithms for DTA and suggested that most DTA algorithms including the MSA (Robbins and Monro, 1951) and its extensions (Lu, Mahmassani, and Zhou, 2009; Sbayti, Lu, and Mahmassani, 2007) only exhibit good performance in small-scale networks and classical DTA problem formulation. They also proposed the PM that has better performance in a medium or large-scale network.

In addition, the ranking process can be avoided.

The PM is a trip-based algorithm based on the swapping probabilities. At iteration n, a fraction of vehicles on non-least-cost paths swaps to the path with the minimum cost. For each path p in OD w, the swapping probability is given by the following formula:

Pr(swap(p) = 1) = C p -min p∈Pw (C p ) C p (2.18)
where C p is the cost for path p and swap(p) is a binary decision variable for swapping.

• Temporal extension of the network

This method involves applying static assignment techniques while performing a temporal extension on each node and link in a network (Peeta and Mahmassani, 1995; Peeta and Ziliaskopoulos, 2001). The network is duplicated as many times as the number of time periods. To prevent missing important network events, the increment of time is dimensioned based on the intended equilibrium's granularity and the flow speed. This method can only be used with small-scale networks.

• Optimal control

According to Delahaye and Puechmorel (2013c), optimal control-based traffic assignment strategies have been put forth in the setting of small-scale networks, which is now well established from a variety of studies [START_REF] Wie | An application of optimal control theory to dynamic user equilibrium traffic assignment[END_REF][START_REF] Terry L Friesz | Dynamic network traffic assignment considered as a continuous time optimal control problem[END_REF]; Ran, Boyce, and LeBlanc, 1993). In the optimal control formulation, the state variable is modeled as the amount or volume of traffic along each link at the beginning of a time period, and the control variable is the flow entering each link over the same period.

• Metaheuristic algorithms

Due to the aforementioned advantages, metaheuristic algorithms have been widely applied to traffic assignment problems, such as Simulated Annealing (SA) ( Furthermore, as mentioned before, traditional optimization algorithms designed for traffic assignment rely on descent direction and converge by swapping the flow from one path to another. The next decision requires that the previous step is completed and the results are obtained, which prevents the process from being parallelized. These drawbacks can be overcome by metaheuristics (Ameli, Lebacque, and Leclercq, 2020b).

Among these approaches, the SA is capable of addressing NP-hard problems with high-dimensional state space with properties of stochastic convergence. The advantage of SA has been thoroughly discussed in (Delahaye, Chaimatanan, and Mongeau, 2019a). The SA is one of the non-population-based metaheuristic algorithms that intensify the search in the local area of the solution space in order to prevent being trapped in local minima. The SA analogizes the physical annealing in metallurgy that imposes heating and cooling on the material to increase the size of its crystals and reduce their defects by affecting its thermodynamic free energy [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF][START_REF] Černỳ | Thermodynamical approach to the traveling salesman problem: An efficient simulation algorithm[END_REF]. Annealing involves two procedures, which are heating and cooling. It aims to change the physical properties of the material. In the heating process, the material is heated to an annealing temperature that allows for atomic rearrangements. Then, the material is cooled sufficiently long to a low-energy state. The temperature should be decreased slowly and carefully to bring the material to the crystalline state, in which atoms are organized in orderly arrangements. Otherwise, if the temperature is decreased suddenly, the material will move to a metastable state, in which the energy is non-minimal. These processes are illustrated in Figure 2.15.

Figure 2.15: The cooling procedures in terms of different temperature-reducing strategies. For a sudden temperature decrease (hardening), the material will move to a metastable solid state with non-minimal energy (top right). The distribution of atoms in this state is not symmetric. For a slow decrease in temperature, the material can be cooled sufficiently long and reach a crystal solid state (bottom right). The atoms in this case are organized and symmetric.

The cooling procedure in the annealing process is simulated by the Metropolis algorithm [START_REF] Metropolis | Equation of state calculations by fast computing machines[END_REF]. This algorithm is based on the Markov chain Monte Carlo method to generate a sequence of random states from a probability distribution. Given a current state i with energy E i , the acceptance probability to

State-of-the-art

transit from state i to j is:

P a =    1, E j < E i e E i -E j k b T , otherwise (2.19)
where T is the temperature and k b is the Boltzmann constant. It can be observed from equation (2.19) that the temperature determines the acceptance probability of a state with higher energy. The higher the temperature, the higher the probability that atoms move to a higher energy state. As temperature reduces in the cooling process, the acceptance probability decreases until almost all atoms are permitted to move to a lower energy state.

In the SA, the Metropolis-Hastings algorithm can be used to generate solutions in the state space. By analogy, the state space is related to the states of the material, and the objective function corresponds to the energy of the material. T is still the temperature, which is considered a control parameter.

The pseudocode of the SA is given in Algorithm 1 and the procedure of each transition is illustrated in Figure 2.16. The SA begins with an initial temperature T 0 and an initial solution γ 0 , which is selected to make sure that transitions with major degradation in terms of the objective function are accepted.

In that way, the state space can be explored homogeneously and thoroughly. It can also prevent being trapped in local optima. N tr transitions are generated at each temperature assuming the state of equilibrium. In each transition, the neighborhood solution is generated by GenerateNeighbor(•) function and accepted or refused according to Accept(•) function. As the temperature decreases by DecreaseTemperature(•) function, fewer degraded solutions tend to be selected. In addition to the formulation in the pseudocode, the best solution so far can be memorized in each loop.

The SA requires the following information to be provided apriori:

-Initial solution

The choice of initial values may affect the convergence speed of the SA.

-Number of transitions N tr

The number of transitions N tr controls the inner iteration at each temperature. Each transition consists of neighborhood generation and acceptance of the neighborhood solution. It can be proved that the SA converges to the global optimum after an infinite number of transitions. In the Metropolis-Hastings algorithm, a large number of transitions are required at each temperature to achieve thermal equilibrium. N tr is also a trade-off between the optimality and CPU time.

-Strategy to generate neighborhood solutions in each state If accepted, the current solution will be updated. Otherwise, the comeback operator returns to the previous solution and avoids duplication in memory.

The neighborhood generation starts from a feasible solution randomly drawn in the state space to search for another solution in the neighborhood of the current solution. The neighborhood generation strategy is a key choice in terms of the SA. The performance of the SA is directly related to the neighbor generation strategy [START_REF] Henderson | The theory and practice of simulated annealing[END_REF]. A proper neighborhood selection strategy can significantly speed up the convergence of the SA algorithm (Tan, 2008). A variety of research states that the asymptotic convergence of SA to the global optimum is dependent on the selection of neighborhood [START_REF] Hajek | Cooling schedules for optimal annealing[END_REF]; Solla, Sorkin, and White, 1986; Fleischer and Jacobson, 1999). Another important factor in choosing a neighborhood generation strategy lies in the size of the neighborhood. Although there is no available theoretical proof, it is well established that if the size of the neigh-

Algorithm 1 Simulated annealing procedure SA(γ 0 , T 0 , N tr ) γ i ← γ 0 , y i ← f (γ i ), k ← 1 do for n ← 1 to N tr do γ j = GenerateNeighbor(γ i ) y j ← f (γ j ) if accept(y i , y j , T k ) then γ i ← γ j y i ← y j end if end for T k+1 ← DecreaseTemperature(T k ) k ← k + 1 while (Stopping criterion is not satisfied) return γ i end procedure
borhood is small in comparison to the cardinality of the solution space, the SA cannot traverse the solution space in a reasonable amount of time to search for the optimal solution. On the other hand, if the size of the neighborhood is quite large, the SA is forced to sample randomly from a large volume of solution space and prevent it from concentrating on a particular area of the solution space [START_REF] Goldstein | Neighborhood size in the simulated annealing algorithm[END_REF]. In sum, the size of the neighborhood is highly problem-dependent. In addition, the structure of the neighborhood and the solution space itself are also issues that need to be addressed [START_REF] Chardaire | Thermostatistical persistency: A powerful improving concept for simulated annealing algorithms[END_REF][START_REF] Fleischer | Information theory and the finitetime behavior of the simulated annealing algorithm: Experimental results[END_REF].

-Cooling strategy

According to the aforementioned discussion, the temperature should be reduced sufficiently and slowly. There are several approaches in terms of temperature reduction: * Geometric process

In the geometric process (Kirkpatrick, Gelatt, and Vecchi, 1983), the temperature is decreased in the following way:

T k+1 = αT k , α ∈]0, 1[, k ∈ Z + (2.20)
where α is the decay control parameter. If α is too large, the temperature will be reduced too slowly and it may take a long computation time to converge. If α is too small, the temperature will be quickly decreased and SA is more likely to be trapped in the local optimum. Typically, α ∈ [0.8, 0.99] (Delahaye, Chaimatanan, and Mongeau, 2019b). * Linear process

In the linear process [START_REF] Kirkpatrick | Optimization by simulated annealing[END_REF], the change of temperature follows a linear decrease with constant coefficients β:

T k+1 = T k -β, β ∈ R + , k ∈ Z + (2.21)
In this case, it should be ensured that the temperature remains positive.

* Logarithm process

The logarithm cooling process is also an efficient strategy [START_REF] Geman | Stochastic relaxation, Gibbs distributions, and the Bayesian restoration of images[END_REF]:

T k+1 = T 0 log(k + d) , d ∈ R + , k ∈ Z + (2.22)
where T 0 is the initial temperature, d is a coefficient that is generally selected as 1.

-Stopping criterion

Theoretically, the execution of the SA terminates when the value of the objective function for the current solution is sufficiently close to the global optimum. However, due to limited computational resources, the selection of the stopping criterion is a trade-off between computational time and the approximation to the optimum. There are several types of stopping criteria commonly used in the SA (Zielinski and Laur, 2007; Zielinski and Laur, 2008), including: kinds of problems. Besides, these models assume that link costs are independent of the volume of other links, which is not realistic. In addition, Wardrop's principles are initially proposed for road traffic, where the flows are naturally segmented by users. There are no prior justifications that these assignment methods are still suitable in terms of non-segmented flow.

Traffic assignment for ATM and UAM

Based on the aforementioned discussion, we can also observe that the development of algorithms for traffic assignment has slowed down in the last two decades. According to [START_REF] Patriksson | The traffic assignment problem: models and methods[END_REF], only a few of them are applied in the ground transportation system. A possible reason is that the ground transportation system has been planned and developed for decades, in which traffic assignment plays a limited role in enhancement and improvement. Rather, in ATM systems or UTM systems, traffic assignment is still quite necessary to improve efficiency as the airspace structure is still evolving. Especially for UAM, the urban airspace is still being conceived and is yet to be developed. The new operation mode and network structure require the UTM for a large number of UAM aircraft in the future. These characteristics are well suited for the use of air traffic assignment models.

Theoretically, air transport enables greater freedom of route choices than most other transport modes. When a flight plan is assigned to aircraft, some flights may choose the same route (due to a similar cost index based on the shortest route) without considering the global optimum of the system. Such route choices may lead to conflicts as well as frequent congestion on commonly flown routes. To avoid such situations, air traffic assignment is widely used to search for a compromise between aircraft and airspace perspective in a way that not too many aircraft are penalized and no aircraft is heavily penalized.

Air traffic assignment is an essential approach to reducing the congestion of airspace. It is a broad term related to flow control by modifying trajectories in terms of time and/or space and respecting the predefined criteria in the meantime. The first model of air traffic assignment is developed as early in 1954. This model addressed the assignment for a given fleet to carry an anticipated traffic load over several routes at the minimum cost. [START_REF] Allen R Ferguson | The allocation of aircraft to routes-An example of linear programming under uncertain demand[END_REF] Delahaye and Puechmorel (2013a) proposed an air traffic assignment model on the 2D air network sectorized into a specific number of sectors in a highly dense area. The objective is to minimize extra route distance and reduce sector workloads considering the equity constraint. A genetic algorithm is used as a solver. A variety of studies have been carried out in terms of air traffic assignment to reduce environmental impacts, for example, to reduce noise near the airport. Netjasov (2008) and Mirosavljević, Gvozdenović, and Čokorilo (2011) focused on the assignment of runways for take-off and landing in terms of aircraft categories to increase airport capacity and reduce the noise level. Chatelain and Van Vyve (2018) aims to minimize the annoyance to citizens while respecting a fairness scheme in terms of a new α-fairness social welfare function. [START_REF] Ganić | Air traffic assignment to reduce population noise exposure using activity-based approach[END_REF] proposed an activity-based approach considering temporal and spatial variations in population in the vicinity of an airport. Using the branch-and-reduce heuristic algorithm, the aircraft was assigned to routes that fly over the least-populated area. This work has been extended (Ho-Huu et al., 2019). A multi-objective optimization framework is formulated by further considering the fuel assumption and solved by a non-dominated sorting genetic algorithm.

In 2007) solved the UAV assignment problem of optimum paths corresponding to the minimum energy costs using Dijkstra's algorithm. However, these studies tended to focus on the tactical perspective rather than the strategic aspect.

Furthermore, to our best knowledge, there is no application of traffic assignment in the UAM from a macroscopic perspective, i.e., the UAM vehicles are characterized by traffic flow. Existing studies have focused on air traffic management with individual aircraft, similar to the current process of air traffic flow management. These models also target low-volume air traffic and cannot be applied to high-density traffic scenarios. In addition, many studies related to air traffic assignment have not focused on underlying causes, such as intrinsic air traffic complexity.

Conclusions

In this chapter, we first reviewed the airspace design for the UAM. As the main challenge with respect to the increase of urban air transportation is the insufficient capacity of the current ATM system, urban airspace is urgently required to be designed for future UAM operations. Highly related to ATM operations, urban airspace is also required to be integrated with current ATM systems. Several important government-led and industry-led initiatives have been introduced. A framework of UAM airspace design will be presented in Chapter 3.

Then, we addressed air traffic complexity metrics. They are referred to as a measurement of the difficulty of the system in charge for managing the air traffic situation in a specific airspace. Furthermore, air traffic complexity metrics are related to the control workload. Control workload is referred to as a measure of the difficulty, of the system in charge of traffic, in managing a situation in a particular airspace. The traffic control system may be a human operator or an automatic process. Numerous control work models have been reviewed. However, the models of control workload have various well-known limitations. Instead, traffic complexity is an intrinsic measurement of the complexity involved in a traffic situation. It is independent of the traffic control system and only dependent on the geometry of trajectories. Typical metrics include geometric metrics [START_REF] Delahaye | Air traffic complexity: towards intrinsic metrics[END_REF], proximity metric, convergence metric (Delahaye and Puechmorel, 2000), clusters metric [START_REF] Histon | Introducing Structural Considerations into Complexity Metrics[END_REF], Grassmannian metric (Delahaye and Puechmorel, 2013b), König metric [START_REF] Essén | Average angular velocity[END_REF][START_REF] Delahaye | A new air traffic complexity metric based on dynamical system modelization[END_REF][START_REF] Juntama | A Distributed Metaheuristic Approach for Complexity Reduction in Air Traffic for Strategic 4D Trajectory Optimization[END_REF]. These basic geometrical metrics are suitable for capturing different kinds of features in terms of complexity. Gathering them together in some way can be much more powerful. Nevertheless, when a traffic situation has a more complex structure, these geometric metrics would fail. Besides, the evaluation of complexity only represents the current situation and cannot reflect the inter-dependency of future times. To overcome these issues, metrics based on dynamical systems are used to quantify the disorder, interaction, and evolution between trajectories in airspace [START_REF] Delahaye | Air traffic complexity: towards intrinsic metrics[END_REF]. Linear Dynamical System has been proven to be suitable and effective for intrinsic air traffic complexity estimation in many studies ( Lastly, the air traffic assignment was presented. As a macroscopic model for optimally assigning traffic flow operated in the route network, traffic assignment has been widely applied to ground transportation systems. The fundamentals of traffic assignment have been first described. To avoid the time-consuming path enumeration or path set generation in solving path-based traffic assignment models, Chapter 3 will address the computation of candidate paths for each OD in order to construct the UAM route network. Then, different algorithms were proposed to optimize the traffic assignment problems. The application of these algorithms to air traffic assignment has been reviewed.

According to the literature review, the DA presents interesting features for solving traffic assignment problems. However, the DA could converge to local optima if the air traffic complexity is evaluated by the LDS, as the cost function has no analytic expression. Therefore, Chapter 4 introduces some randomness in the DA. Another possibility is to use the metaheuristic together with the DA, which is presented in Chapter ??. Since the DA is computationally expensive in computing the swapping flow and marginal cost, it is only feasible for small-scale networks. To overcome this problem, Chapter 6 introduces a simulated-based framework involving rolling horizon and parallel simulated annealing algorithm to solve high-density traffic assignment problems on large-scale UAM route networks.

Chapter 3

UAM route network design in low-altitude airspace

Problem description

This chapter addresses designing UAM route networks in low-altitude airspace to minimize the noise impact and maximize the efficiency and safety of UAM operations. Various open-source data are used to propose a realistic solution for UAM route network design. The UAM route network provides a set of candidate paths for each OD. As a country with 100% urban population (Chizho, Okuneva, and Hillawi, 2021), Singapore is prosperous in the development and implementation of the UAM. We select Singapore's urban airspace as the case study.

The definitions of the variables and notations used in this chapter are listed in Table 3.1.

Data source and data description

We first use the ALOS World 3D-30m (AW3D30) data (Japan Aerospace Exploration Agency (JAXA), 2022) to model static obstacles in the urban airspace of Singapore. It is a global Digital Surface Model (DSM) with a resolution of approximately 30m. The terrain is described in this data, including natural and built features. The ALOS World 3D data that covers Singapore is visualized in Figure 3 Then, the airspace classification information of Singapore (CAAS, 2022) is used in this study to identify airspace that is potentially hazardous to UAM operations. These airspace types are promulgated as follows (Civil Aviation Authority of Singapore (CAAS), 2015; CAAS, 2019):

• 5km of aerodromes: The airspace within 5km of the aerodrome area.

• Danger areas: The airspace within which activities dangerous to aircraft may exist at specified times.

• Prohibited areas: The airspace within which aircraft is prohibited, usually due to security purposes.

• Restricted areas: The airspace within which aircraft is restricted in accordance with certain specified conditions.

• Protected areas: The reserved airspace within which aircraft is not allowed to operate.

To sum up, UAM aircraft are not allowed in restricted and protected airspace. All types of airspace are shown in Figure 3.2 bounded by discrete points. In addition, the vertiports and delivery sites are indexed by numbers. In addition, the location of vertiports is collected from a study that determined the optimal distribution of vertiports for Singapore considering both costs and society value (Zeng, Low, Schultz, et al., 2020). Some locations of vertiports are modified to avoid hazardous airspace and terrain. Some delivery sites are randomly generated in usable airspace. The location of vertiports and delivery sites are depicted in Figure 3.2.

The Singapore residents by planning area data as of June 2021 (Department of statistics Singapore, 2021) and the boundary of planning area (Government of Singapore, 2019) data are collected to measure the noise impact of the UAM network on populations in different regions. Singapore is divided into 332 subzones, and the total resident population in each subzone is illustrated in Figure 3.3. Note that in some subzones, the population can be nil or negligible.

Methodology

In addition to safety and efficiency, the environmental effect of aircraft movements is one of the essential factors in the acceptance of new UAM route networks, droneports, air routes, or new aircraft types. The impact on the environment can be divided into two types:

• Traditional environmental impacts: These impacts come from traditional air transportation and have been thoroughly assessed, such as emissions, pollution, noise, third-party risk, etc.

• Novel environmental impacts: These impacts are newly discovered. They are still unclear and need further assessment. These impacts include light pollution, privacy concerns, distraction, effects on the biotope, etc.

Among these impacts, the following metrics are proposed:

• Third-party risk metric: Third-party risk is the risk of the third party, for example, the people on the ground, inadvertently exposed to an aircraft accident. The thirdparty risk metric is based on the existing accident data in terms of different aircraft types in different phases.

• Emission and pollution metric: The emission and pollution metrics are also called energy usage metrics because the total energy usage is estimated rather than di-rectly measuring the pollution and emission. The energy usage can be computed based on aircraft weight and flight hours or thrust and traveled distance.

• Noise metric: The more time aircraft are at a lower altitude, the more noise is imposed on the ground. The noise metric is also related to the flying time below a certain altitude.

In the low-altitude urban airspace of a metropolitan area, the problem addressed here aims to design a multi-layer UAM route network that optimizes the airspace usage, flight efficiency, and noise impact on residents, while keeping away from terrain and hazardous airspace. Based on the aforementioned datasets, we first determine the minimum spatial area to design the network and partition the area into grids. The grid points and air routes that have no interaction with the terrain and hazardous airspace are considered as nodes and links, respectively. We search for the paths with k-minimum cost and with diversity between each OD to construct the route network. The network design relies on horizontal 2D grids. Concerning the regulations for UAM operations in low-altitude urban airspace (Cho and Yoon, 2018b), a vertical separation standard of 25ft is ensured between network layers. Then we create network layers on altitudes 150ft, 175ft, 200ft, 225ft, and 250ft. In different layers, we consider the terrain at and above the current altitude as static obstacles. To maintain safety separation between UAVs and obstacles, we apply geofence on terrain to impose geographical boundaries around the perimeter of obstacles as a virtual barrier. We define the minimum keep-out distance as d k =30m in accordance with the specification of most countries (Cho and Yoon, 2018b), as illustrated in Figure 3.4. Consequently, the usable airspace is defined at each altitude.

Construction of the grid-based network

To start with, World Geodetic System (WGS) is converted to Projected Coordinate System (PCS) EPSG:3414 for Singapore. To reduce the search space and the network construction cost, we calculate the Axis-aligned Minimum-Area Bounding Box (AM-ABB) of vertiports and delivery sites. We use the method of rotating calipers (Toussaint, 1983) to calculate the AMABB. The idea is to calculate the convex hull, then rotate by each edge orientation and compute the bounding rectangle. The minimum intersection The obstacles with geofence in the construction area at different altitudes are shown in Figure 3.6. The urban airspace, especially below 40m, contains a large number of geographical constraints, which seriously limits the connections between nodes in the grid-based network. This in turn causes bottlenecks of congestion, especially at intersections, which reduces overall airspace capacity and performance. As a remedy, the grid-based network has a high level of flexibility and connectivity in the airspace at higher altitudes, where few or no static obstacles still exist.

At each altitude, the grid points in the usable airspace and the points that have the same longitude and latitude as vertiports, and delivery sites are selected as potential waypoints. Any pair of potential waypoints with a distance not greater than √ 2d g and not colliding with unusable airspace and obstacles are connected with a horizontal link. In addition, the vertiports and delivery sites, and all waypoints directly above them are connected with vertical links. All these links form the set of link E. All nodes that connect links in E form the set of node V. The resulting network is represented as a graph G = (E, V).

Constraints and costs for candidate paths

In order to support high-density and complexity UAM traffic flow in the grid-based network G, the candidate paths P = {P w |w ∈ W} should be determined for a given set of OD pair W. We search for candidate paths that maximize flight safety and efficiency, minimize the noise impact on populations in the nearby subzones, and satisfy operational constraints in G. The constraints are related to the range of UAV which is a measure of distance, and the endurance of UAV which is a measure of time.

Firstly, the length of path d p should not exceed the range of UAV, denoted as R max :

d p < R max , p ∈ P (3.1)
In addition, the flight time should also be limited to the endurance of the UAV:

e∈p d e s e < T max , p ∈ P (3.2)
where s e is the average speed of UAM traffic flow on link e. Note that maintenance and battery charging or refueling is assumed to be available both at vertiports and delivery sites. If not, the range and endurance of the UAV should be cut by half to consider the return trip.

It is noteworthy that some parts of usable airspace may still be at risk, including aerodromes A a , danger areas A d , and prohibited areas A p . The link traversing these areas is characterized by δ e,m and penalized by a coefficient α m , m ∈ {a, d, p}. If link e traverse hazardous airspace m, δ e,m = 1, otherwise δ e,m = 0. Then, we define the cost for crossing hazardous airspace C h,e for link e ∈ E based on the airspace type:

C h,e = (α a δ e,a + α d δ e,d + α p δ e,p ) (3.3)

UAM route network design in low-altitude airspace

The cost of flight efficiency is calculated based on link length d e and the type of link:

C f,e = (ϕ h δ e,h + ϕ u δ e,u + ϕ d δ e,d )d e (3.4)
where d e is the length of link e, δ e,h = 1 if link e is horizontal, otherwise 0. δ e,u and δ e,d similarly for upward and downward links. ϕ h , ϕ u , ϕ d are efficiency coefficients for horizontal, upward, and downward links, respectively. The coefficient is much larger for vertical links than horizontal links because crossing altitudes will introduce traffic complexity with other levels. The link for climbing is more penalized than the link for descent as it is more energy-consuming. Another focus is to mitigate the noise impact of the designed link to populations in each subzone of Singapore. We use Sound Pressure level (SPL) to measure the pressure level of noise produced by UAVs (Scozzaro, Delahaye, and Vela, 2019). We assume the UAV at height h produces the same reference noise of L h at the ground point directly below, as shown in Figure 3.7. Based on the inverse proportional law, the SPL L r at a distance r from the point directly under UAV can be simplified as:

L r = L h + 10 log 10 h 2 h 2 + r 2 (3.5)
The noise impact of the UAV traffic flow on a link e to a nearby subzone S can be represented as the integral of the SPL from a point distributed on the link to a point distributed in the subzone: Many geographic regions, including census tracts, block groups, and blocks, are largely delineated by streets and are close to the rectangle in shape (U.S. Census Bureau, 1994). To efficiently calculate the multiple integrals in L e,S , the subzone S can be approximated by an AMABB (Mu and Tong, 2020; Mu and Tong, 2021). The grid-based network is rotated by a specific angle according to the reference point of the AMABB so that the Minimum Bounding Box (MBB) becomes the AMABB.

L e,S =
After the rotation, the notation for calculating the noise metric of horizontal links and vertical links is given in Figure 3.8. If the link is horizontal, it is characterized by start point (x 0 , y 0 ), end point (x 1 , y 1 ), and a constant height h. Else if the link is vertical, it is represented by the start point, end point, the minimum height h min , and the maximum height h max . The AMABB of the subzone S is bound by x max , x min , y max , and y min . We give the exact expression of noise impact from a link to the nearby subzones:

If the link is vertical, which corresponds to

L h + 10log 10   h 2 h 2 + (x -u) 2 + (y -(u-x 0 )(y 1 -y 0 ) x 1 -x 0 -y 0 ) 2   dudxdy (3.9)
If the link is horizontal and the x-coordinates of both endpoints are the same, L e,S = max(y 0 ,y 1 ) min(y 0 ,y 1 )

xmax x min ymax y min L h + 10log 10 h 2 h 2 + (x -x 0 ) 2 + (y -v) 2 dvdxdy
(3.10) The numerical integration method is used to solve these integrals. We search for subzones with polygonal shapes in the vicinity of each link modeled as a line segment. The total noise impact from a link e to all nearby subzones is defined as:

C where x ∈ {h, f, n}. The path cost is the sum of the costs of links that belong to this path:

C p = e∈p C e (3.15)

Finding k-shortest path

Denote the acyclic paths satisfying constraints (3.1-3.2) as P, then the problem is to find for each OD w the set of k paths P w satisfying ∀p ∈ P/P w , ∀q ∈ P w , C p ≥ C q .

We use a greedy algorithm introduced by Liu et al. (2017) to approximate the KSPD problem. The pseudocode of the algorithm to generate candidate path sets for all ODs is given in Algorithm 2. The idea is to insert the next shortest path that is sufficiently dissimilar to all shortest paths in the result set. Among similarity functions in the 3.3. Methodology 57 literature [START_REF] Liu | Finding top-k shortest paths with diversity[END_REF], the following metric is used to measure the similarity between two paths p and q based on the path cost:

sim(p, q) = C Ep∩Eq 2C p + C Ep∪Eq 2C q (3.16)
where E p , E q is the set of links on path p and q, respectively. Note that the similarity measurements proposed in the literature are essentially the same and are based on path cost, only the form is different. Other similarity measures can also be chosen according to different requirements. If the similarity between two paths is higher than a threshold, they are considered similar paths and cannot exist at the same time. This algorithm has been shown to have a high approximation ratio to the exact solution in real-world applications (Liu et al., 2017).

Constructing the route network

Finally, the candidate path sets for all ODs generate the route network. A link or a node in the grid-based network is removed from the route network if none of the

UAM route network design in low-altitude airspace

candidate paths pass through it.

Result and discussion

Experiment setting

The size of discretization is set as d g = 2000m to reduce the computation time while providing a suitable amount of node and link options. According to the hazard level of airspace (Civil Aviation Authority of Singapore (CAAS), 2015; CAAS, 2019), α a , α d , and α p is set to be 0.2, 1, and 1, respectively. Considering the distance gap between horizontal and vertical links, we set ϕ h = 0.005, ϕ u = 1 and ϕ d = 0.5 to make their costs of flight efficiency comparable. According to acoustical measurement data of UAV in [START_REF] Schäffer | Drone Noise Emission Characteristics and Noise Effects on Humans-A Systematic Review[END_REF], we roughly set

L h = (80 -2/25(h -150))dB, h ∈ [150ft, 250ft].
More precise data should be obtained through flight experiments. d t is set as a short distance of 2000m to avoid significant effects from other impact factors such as wind. Referring to [START_REF] Troudi | Sizing of the drone delivery fleet considering energy autonomy[END_REF], the range and endurance of UAV are set as R max = 25km and T max = 1h. The average speed of UAM traffic flow on horizontal links and vertical links is set to be 15m/s and 2m/s. For each OD pair, we select the number of shortest paths k = 20. The similarity threshold τ p is chosen to be 0.6. The framework is implemented in Python 3.7 on a laptop with an Intel i7-8750H CPU and 32GB DDR4 RAM.

Construction of the basic grid-based network

The air transport network presented in this section corresponds to the result of Section 3.3.1. The resulting multi-layer grid-based network is shown in Figure 3.9. The vertiports, delivery sites, and waypoints are marked in the figure. The links in different altitude layers are marked in different colors. The resulting network contains 5 altitude layers. To compare the horizontal network in different altitude layers by removing vertical connections, some basic topology features are calculated on the basis of network theory. These topology features include the number of nodes, number of links, Average Degree (AD), Average Clustering Coefficient (ACC) [START_REF] Saramäki | Generalizations of the clustering coefficient to weighted complex networks[END_REF], and Degree Of Assortativity (DOA) (Newman, 2003). The ACC is approximated by repeating the following experiment a specific number of times: randomly choosing a node and randomly selecting two of its neighbors, then checking whether they are connected. The ACC is then calculated as the ratio of the found triangles to the number of trials. Here we repeat 1000 times trials. The DOA is calculated as the standard Pearson correlation coefficient of the degree between the linked nodes. Their values for different altitude layers are summarized in Table 3.2.

As shown in Table 3.2, the grid-based network has more nodes and links at higher altitudes because higher altitude layers have fewer obstacles. The result of the AD indicates that each node has more connections to neighbors as altitude increases. The growth of the ACC with altitude signifies that nodes tend to cluster together as altitude increases. Except for the altitude layer of 150ft which is non-assortative, the DOA for other network layers is positive, which means that high-degree nodes are more likely to attach to high-degree nodes, that is, these networks are assortative.

Figure 3.9: The resulting multi-layer network in Singapore urban airspace.

Link costs

The normalized link cost for crossing hazardous airspace Ch,e , of flight efficiency Cf,e , and of noise impact Cn,e are calculated and represented in Figure 3.10, 3.11, and 3.12, respectively.

According to Figure 3.2, since hazardous airspace occupies a large amount of Singapore's airspace, most of the links in the grid-based network have to traverse hazardous airspace. This fact is consistent with Figure 3.10. Nevertheless, airspace within 5km of aerodromes takes up the largest proportion of hazardous airspace, so most links have a low cost for traversing hazardous airspace. There are some areas where different types of hazardous airspace overlap, in which link costs are high. In terms of link cost of flight efficiency shown in 3.11, for horizontal links, longer links are less efficient and have higher costs. Though vertical links are heavily penalized, it can be seen that the costs of flight efficiency for vertical links are still less than for horizontal links. The main reason for this observation is that the length of vertical links is 25ft, far less than the length of horizontal links, which is 2000m. Also note that the flight efficiency costs for upward links are two times larger than for downward links. In Figure 3.12, links with the high cost of noise impact are computed. Referring to Figure 3.3, they are mainly distributed around densely populated areas.

UAM route network generation

To compare the relative importance of the three types of link cost in Equation (3.13), the different combinations of weighting parameters are examined. The route network is generated with weighting parameters w h , w f , w n [1/3, 1/3, 1/3], [2/3, 1/6, 1/6], [1/6, 2/3, 1/6], and [1/6, 1/6, 2/3]. In these four cases, the first parameter setting treats all link costs as equivalent. Other cases consider that one of the link costs is more important than the other link costs.

As the previous steps can be performed offline, we focus on the computation time of candidate path generation. For these four cases, the candidate path generation costs 11.19s, 11.74s, 13.39s, and 11.31s, respectively. It can be seen that these computation times are short in terms of the relatively large scale of route networks, which can support efficient constructions of the UAM route network. In addition, the average cost for kshortest paths between OD pairs is shown in Figure 3.13. The coordinates indicate the index of vertiports and delivery sites. For OD pairs involving some vertiports, such as the vertiports 0, 1, 2, and 3, the average path cost for different weighting parameters remains low. The reason is that they are more accessible to most delivery sites through the UAM route network. On the contrary, the average path costs for OD pairs that involve vertiports 4, 5, and 6 are always at a high level. In Figure 3.13b, the link cost for crossing hazardous airspace is considered more important. In this case, the peak value of the average path cost is the lowest among the four cases. The reason can be found in Figure 3.2, i.e., most of the hazardous airspace is covered by aerodromes, which receive fewer penalties than others. Figure 3.13c corresponds to the case in which the link cost of flight efficiency is more important. The pattern of average path cost, in this case, is similar to the first case, shown in Figure 3.13a, but the peak value is much higher. This fact indicates that candidate paths between some OD pairs are less efficient, such as OD pairs in the lower left part of Figure 3.13c. In the last case, it can be seen from Figure 3.13d that paths connecting vertiport 4 have high costs. The reason lies in the fact that vertiport 4 is located in a densely populated downtown area. Next, we visualize the route networks for these four cases in Figure 3.14. The color indicates the number of times that each link appears in the candidate paths. At a quick glance, all altitude layers in these four networks are utilized. Vertical links, as the connection of adjacent altitude layers, appear most often in the route network. According to different weighting parameters, the route network is created by minimizing the criteria formulated in Equation (3.15). Slight differences, including link counts and network topology, can be observed in these networks. To sum up, the resulting UAM route network provides many travel alternatives, which can be used by air traffic flow

Conclusion 63

management and air traffic assignment to mitigate the complexity and congestion of airspace. Note that the dead ends that can be observed in Figure 3.9 to 3.12 have been successfully removed in Figure 3.14.

Conclusion

In this chapter, we have considered the problem of UAM route network design. A methodology has been proposed that aims to define a set of alternative routes between numbers of OD pairs in the UAM context. Based on various open-source data, the obstacles and hazardous airspace are identified and modeled. The grid-based network is firstly designed to avoid obstacles and airspace within which the aircraft is prohibited. Then, we define the link cost in terms of noise impact on populations, airspace safety, and flight efficiency. Unlike most research that only calculates the shortest path, relevant constraints and costs are involved in searching for k-shortest paths with diversity. The feasibility of this approach is demonstrated by its application to a parcel delivery scenario for Singapore. The framework is easy to generalize to other urban airspaces all over the world. The proposed methodology can generate a high-quality UAM route network to support high-density and complexity flow-based operations using the sets of paths with low similarity resulting from the solution to the KSPD problem.

A conference paper related to this chapter has been published (Wang et al., 2022c).

Chapter 4

Static air traffic assignment model in 2D UAM route network with segregated layers

Problem description

To adapt the increasing demand to the current airspace capacity, this chapter focuses on strategic UAM operational planning through mitigating congestion and organizing the flow structure from a macroscopic perspective. Methodologies are designed for highdensity UAM operations in a 2D route network with segregated layers.

The definitions of the variables and notations used in this chapter are summarized in Table 4.1.

2D UAM route network modeling

To adapt to high-density demands, we choose a centralized UAM route network for UAM operations (FAA, 2020a; Goodrich and Barmore, 2018).

Graph representation of UAM route network

Consider a geographical area A ⊂ R 2 corresponding to a certain flight level in urban airspace. A UAM route network in A is built with a direct graph G = (V, E) with V being the set of nodes and E being the set of links. The average speed of aircraft on link e is set as s e . The link connecting nodes i and j is indexed as e i,j . The coordinates of nodes are {(x v , y v ), v ∈ V}.

Let W be the set of OD pairs, in which each element w = (O, D) represents the origin node O ∈ V and the destination node D ∈ V. P = w∈W P w denotes the set of feasible paths for all OD pairs, where P w = {p 1 w , ..., p Rw w } represents all feasible paths that connects a specific OD pair w. All possible paths without dead ends and cycles between OD pairs are obtained based on depth-first search (Migliore, Martorana, and Sciortino, 1990). Moreover, we introduce two logical functions to characterize the relationship between link and path. The first type of logical function indicates the path index and the OD pair the link belongs to:

δ e,w,i =    1 if e is contained in the i th path of OD pair w 0 else (4.1)
Regardless of the detail of the OD pair, the logical function of the second type only reflects the relationship between a link and a path:

δ e,p =    1 if e is contained in path p 0 else (4.2)
The demand between OD pairs is denoted by D = {d w , w ∈ W}, which can be represented as a vector.

UAM corridors modeling

Based on the UAM corridor concept introduced in Section 2.1, for link e ∈ E, we design the related UAM corridor U e ∈ A as a rectangle route segment with the same length as link e. To adapt UAM corridors to a high-density operational perspective, two important terms, including density points and complexity area, are introduced.

Density points

To describe the traffic density of UAM operations from a macroscopic perspective, density points are introduced to represent the flow distribution along the UAM corridor. A set of density points ρ e is defined as virtual grids located on each UAM corridor. Each density point g is characterized by coordinates (x g , y g ), heading θ e and speed s e inherited from the link e. A certain number of aircraft is considered at each point, which is also referred to as density and is dependent on the link flow. For every longitudinal position, the lateral density points are assumed to follow a uniform distribution of flow. Let the number of points on each side of a link be N s . To equally share the flow in the lateral direction, the lateral flow sharing coefficient α e is thus defined as:

α e = 1 2N s + 1 (4.3)
Due to the varying length of the UAM corridor, the distance of successive density points is not possible to be the same in the longitudinal direction. An alternative grid partition approach is formulated in Algorithm 3. For all UAM corridors, the size of grids is first assumed to be the same and denoted as L l × L w . The objective is to set the longitudinal distance between successive density points to L l with the exception that one longitudinal distance in the middle of the link belongs to [L l , 2L l [, which corresponds to lines 4 to 11 in Algorithm 3. Though it is possible to distribute this excess space around the margin of the UAM corridor, the choice of concentrating the excess space in the middle of the link is more justified because it is the area where the complexity is likely to be the lowest, and thus the bias is the smallest. To facilitate the computation, the grid partition is firstly conducted on a predefined coordinate system, then translated and rotated to the original direction, which corresponds to Figure 4.1 and lines 3 and 14 in the pseudocode.

Finally, following the uniform distribution of flow, the speed s g for all density points in ρ e is set to the same value s e .

Complexity area

In order to measure the air traffic complexity of the UAM corridor, not only its traffic should be included, but also the surrounding airspace. The complexity evaluation

Algorithm 3 Grid partition on UAM corridor U e

Input: (x i , y i ), (x j , y j ): Coordinate of start and end node of link e N s : Numbers of parallel points on each side L l , L w : Grid length and width Output: P e g : coordinates of density points θ e : common heading for all density points

1: procedure grid_partition(x i , y i , x j , y j , N s , L l , L w ) 2: d ← (x i -x j ) 2 + (y i -y j ) 2 3: θ e ← atan2(y j -y i , x j -x i ) 4: n l ← floor(d/L l ) 5: if mod(n l , 2) == 1 then 6:
x s ← linspace(0, (n l -1)L l /2, (n l + 1)/2) 7:

x e ← linspace(d, d -(n l -1)L l /2, (n l + 1)/2) 8: else 9:
x s ← linspace(0, n l L l /2, n l /2 + 1) 10: return P e g , θ e 16: end procedure airspace for each UAM corridor is called the complexity area. For the UAM corridor U e , e ∈ E, the complexity areas {A e,1 , . . . , A e,Ne } are defined as disks with radius R e , where N e is the number of complexity areas for U e . Each area A e,i includes a set of density points ρ e,i with ρ e,i ∩ ρ e ̸ = ∅. In general, ρ e,i ̸ ⊂ ρ e because the complexity area ρ e,i may include density points of other links e ′ ̸ = e. Figure 4.2 illustrates a case where the complexity areas are disks and the density points are located on grids aligned with the link. The centers of discal complexity areas are also longitudinally distributed between two nodes except for the extra space in the middle of the link. Those centers are referred to as evaluation points Q e = {q e,1 , . . . , q e,Ne } and can be determined by Algorithm 3 under L w = 0, and N s = 0. All points in the i-th complexity area of link e, ∀q ∈ A e,i , ensures that the euclidean distance between q and q e,i is not greater than the radius R e . The set of links whose UAM corridors intersect with the complexity area A e,i is denoted as E q e,i . Figure 4.2 shows a small number of point-to-point UAM corridors. The complexity area of each UAM corridor is presented as a grey disk. From left to right, the complexity areas involve 4, 3, and 4 links, respectively. The density points contained in each complexity area are specifically shown in detail.

x e ← linspace(d, d -n l L l /2, n l /2 + 1) 11: end if 12: y ← linspace(-N s L w , N s L w , 2N s + 1)

Mathematical model of 2D SATA problem

The problem addressed in this chapter can be stated as follows: given a route network within UAM operations, the objective is to find the complexity-optimal flow distribution for UAM corridors such that the congestion and air traffic complexity of the route network is minimized. In future UAM operations, air vehicles will follow the instructions delivered by external control systems, which involve cooperation between users. The total costs incurred by all users can be reduced to the minimum, corresponding to system equilibrium in Definition 2.3. The optimization problem is then formulated as follows:

Decision variables

The flow on the feasible paths of the w-th OD pair is denoted as F w = {F 1 w , ..., F Rw w }, which can be represented as a vector. The decision variables are flows F i w on all paths indexed by i ∈ {1, . . . , R w } for all OD w ∈ W.

Constraints

For w ∈ W, the path flows satisfy the conservation law, i.e., the sum of flows on all paths between each OD pair is equal to the demand on this OD:

d w = Rw i=1 F i w (4.4)
Moreover, the flow on all paths must not be negative:

F i w ≥ 0 ∀w ∈ W, ∀i ∈ {1, . . . , R w } (4.5)
Given the flow of all admissible paths, the flow on a link e ∈ E is given by:

f e = w∈W Rw i=1 δ e,w,i F i w (4.6)
To share the flow for density points, the traffic flow of aircraft within UAM corridor U e is assumed with average speed s e and heading θ e . The density of each density point g ∈ U e is defined as:

ω g = β α e f e s e , ∀g ∈ U e (4.7)
where β is a parameter to convert ω g to a convenient unit and adjust the scale. The density is considered to be proportional to the flow f e and the lateral flow sharing coefficient α e , and inversely proportional to the average speed s e .

Objective function based on linear dynamical systems

To evaluate the air traffic complexity, the LDS is taken into account in the link cost function, in which more intrinsic information will be involved, such as flow, information of density points, and configuration of nearby links around link e. In addition, the LDS is able to quantify the disorder, interaction, and evolution between trajectories in the airspace. As a strategic planning model, when the complexity metric based on the LDS becomes lower, the flow will become more organized, and the conflicts or interactions in real UAM operations will be reduced. The safety and efficiency of an autonomous UAM system can be guaranteed. In order to compute the complexity on the whole graph, the complexity of each complexity area A e,i with i ∈ {1, ..., N e } and e ∈ E are computed using the LDS. The density, position, heading, and speed of density points in each complexity area are used to determine the coefficient matrices and parameters of the corresponding LDS that best fit the density points. Then, the calculated dominant eigenvalue and the flow on the associated links can be used to construct the link cost function.

An LDS in geographical area A can be expressed as a differential equation of motion:

v = Ax + b (4.8)
where x denotes the position vector, v represents the speed vector, A stands for a coefficient matrix, b is a coefficient vector.

Conventional implementations of the LDS simply concatenate the positions and speeds of all points in the complexity area A e,i [START_REF] Delahaye | Air traffic complexity: towards intrinsic metrics[END_REF]. To calculate the error criterion much more efficiently, a Weighted Minimum Mean Square Error (WMMSE) estimator is proposed to adjust the parameters of the LDS:

E e,i = g∈A e,i ω g ∥v g -(A e,i x g + b e,i )∥ 2 (4.9)
where the density ω g , position vector x g = (x g , y g ) ⊤ , heading θ g , speed s g , and speed vector v g = s g (cos(θ g ), sin(θ g )) ⊤ are inherited from the UAM corridor they belong to.

It is worth noting that density points contained in a certain complexity area may be derived from one UAM corridor or multiple UAM corridors.

The complete procedure for calculating coefficient matrices of the LDS is provided in Appendix A.1 and A.2. After determining A e,i and b e,i , the LDS is capable of generating a vector field that fits the observations with the minimum WMMSE. For instance, three detailed coordinate systems at the bottom of Figure 4.2 illustrate the vector fields generated by the LDS based on density points in each complexity area. The arrow length is proportional to the velocity. The color of density points corresponds to values in the colorbar. It can be seen that the goodness of fit in terms of the vector field generated by the LDS is quite satisfying.

The matrix Âe,i that minimizes the WMMSE can be derived. Let λ 1 e,i and λ 2 e,i be the two eigenvalues of Âe,i respectively. The convergence or divergence property of the system can be determined in terms of such eigenvalues. Eigenvalue with a positive real part corresponds to a system in expansion mode. If the real part is negative, the system is in contraction mode. The absolute value of the real part represents the level of contraction or expansion of the system. Besides, the imaginary part of the eigenvalues, Im(λ 1 e,i ) and Im(λ 2 e,i ), represents the curl organization of the system, i.e., the rotation tendency of the vector field. In summary, an LDS can evolve in contraction, expansion, rotation, or their combinations. The basic evolution of the LDS is presented in In the proposed model, the convergence and divergence cases are both considered unorganized structures. The dominant eigenvalue λe,i is then proposed to evaluate the evolution tendency of the system. It is defined as the largest absolute value of the real parts of the two eigenvalues:

λe,i = max(|Re(λ 1 e,i )|, |Re(λ 2 e,i )|) (4.10)
The absolute value operation unifies the convergence and divergence cases. The more the value of the dominant eigenvalue is closer to zero, the more the traffic situation is organized.

To involve the air traffic complexity measurement defined by the LDS and the traffic flow of neighboring links, the cost for link e is designed as:

c e = 1 N e Ne i=1
λe,i

e ′ ∈Eq e,i f e ′ (4.11)
This cost function corresponds to the weighted sum of the dominant eigenvalue of the LDS in each complexity area multiplied by the sum of flow for evaluated links. That is, the complexity is weighted by the flows involved in the complexity area.

Finally, the optimization problem can be formulated as follows:

min f ={fe|e∈E} S (f ) = e∈E 1 N e Ne i=1
λe,i

e ′ ∈Eq e,i f e ′ s.t. f e - w∈W Rw i=1 δ e,w,i F i w = 0, e ∈ E F i w ≥ 0, w ∈ W, i ∈ {1, ..., R w } Rw i=1 F i w = d w , w ∈ W (4.12)

Solution approaches

Dafermos' algorithm (DA) is used as a solution approach. For details, please refer to Section 2.3.2.2. Applying DA to this problem requires overcoming two difficulties:

1. Computation of marginal costs ∂ca ∂f b ;

2. Optimization problem to determine the transfer flow σ w .

In this problem, the cost function depends not only on the flow pattern but also on the air traffic complexity metric, which is the dominant eigenvalue of the LDS in terms of complexity areas. The representation of the cost function does not enable us to find a closed-form expression of the derivative. Numerical differentiation techniques have been adopted to estimate the marginal cost. For a, b ∈ E, the derivative of c a with respect to f b is calculated by finite difference approximations:

∂c a ∂f b = lim h→0 c a (f b + h) -c a (f b ) h (4.13)
where h represents a small variation in f b . It is noteworthy that according to Equation (4.7), f b is proportional to ω g . Since the density ω g ∈ R + makes sense in the weighted form of Equation (4.9), h can be taken as sufficiently small as possible. Based on Dafermos and Sparrow, 1969; Dafermos, 1972, the analytical representation of σ w was easily computed with the assumption of a quadratic form of the cost function. However, it is not the case here. Moreover, for the sake of convergence, the overall global cost should be convex with regard to link flows. In our problem, the dominant eigenvalue is involved in the cost function, and it is dependent on link flow. It could be very difficult to compute the exact mathematical derivation of the cost function, which cannot assure its convexity. Besides, the scale of dominant eigenvalues depends on the traffic situations, which may result in large differences between the marginal costs of each path. The path with the maximum or minimum marginal cost may not change after transferring the flow, making the whole algorithm get stuck in a local optimum. To this end, we propose a two-step strategy to improve the DA. The flowchart of this approach is presented in where arg min(c ′ p , N p w,min ) and arg max(c ′ p , N pw,max ) are defined to return paths with the N p w,min smallest and the N pw,max largest marginal cost under given conditions, respec-tively. random() picks a random path in the candidate set. A possible way to determine the values of N p w,min and N pw,max can be:

N p w,min = min(R w -1, N r ) (4.16)
N pw,max = min(card({p|F p > 0, p ∈ P w }) -1, N r ) (4.17)

where R w -1 ∈ Z + and card({p|F p > 0, p ∈ P w }) -1 ∈ Z + . These two variables control the lower limit of N p w,min and N pw,max . N r is a parameter to limit the maximum number of arguments to choose. By this means, the paths p w,min and p w,max are randomly chosen from sets of paths with maximum and minimum marginal cost rather than choosing paths with exact maximum and minimum marginal cost. It forces to transfer flow between paths, which helps escape from local optima. An iteration of the modified DA is depicted in Figure 4.3 given the random select number N p w,min = N pw,max = 2.

To solve the optimization problem related to flow transfer formulated in Equation (2.13), a Random Optimization (RO) strategy is used. As a numerical optimization method, this strategy is suitable for discontinuous and non-differentiable functions. The pseudocode is shown in Algorithm 4 and is also illustrated in Figure 4.5. This approach proceeds N t epochs in each iteration, by moving to other positions in the search space and sampling with a given distribution D if it is a better position.

The whole resolution algorithm is terminated if any of the following three criteria S min ← +∞ 1. Limit on the number of iterations: when the total number of iterations exceeds the maximum number of iterations N m .

2. Number of iterations without progress: when the number of iterations with no progress on the global cost function exceeds the number of patience N p .

3. Convergence test: when the value of the global cost function is sufficiently small in some sense, lower than a threshold of v s .

Results

Scenario definition

To demonstrate the applicability and effectiveness of the proposed algorithm, two planar UAM route networks in segregated altitude layers are introduced for the experiment. These two networks are denoted as G 1 and G 2 . G 1 is a network modified on the basis of a transportation network for last-mile delivery in Maribyrnong, Melbourne [START_REF] Ewedairo | Estimating transportation network impedance to last-mile delivery[END_REF]. This transportation network is able to cover the suburban and urban areas of most cities in the world. G 2 is a real air transportation network in ATM operations in (Delahaye and Puechmorel, 2013a). G 2 is larger in size and scale than G 1 .

The coordinates of nodes and the index of nodes and links are illustrated in Figure 4.6. The topological characteristics of G 1 and G 2 are summarized as follows. The first subscript indicates the index of the network, and the second subscript indicates the index of the associated attribute.:

• Set of nodes V 1 = {v 1,1 , ..., v 1,17 }; V 2 = {v 2,1 , ..., v 2,24 }, • Set of links E 1 = {e 1,1 , ..., e 1,28 }; E 2 = {e 2,1 , ..., e 2,53 } • Set of OD pairs W 1 = {w 1,1 , w 1,2 } with w 1,1 = (v 1,1 , v 1,17 ), w 1,2 = (v 1,4 , v 1,7 ); W 2 = {w 2,1 , w 2,2 } with w 2,1 = (v 2,16 , v 2,24 ), w 2,2 = (v 2,1 , v 2,22 ).
All feasible paths are calculated for G 1 and G 2 :

• Set of feasible paths P 1 = P w 1,1 ∪ P w 1,2 with P w 1,1 = {p The distribution of travel demands is initialized by means of distributing the demands equally among the allowable paths. For w ∈ W 1 ∪ W 2 ,

F i w = d w R w , ∀i = 1, ..., R w (4.18)
The parameters used in UAM corridors design and resolution algorithm for optimization problem (4.12) are respectively given in two parts of Table 4.2. These two networks share the same parameters.

Results

Performance of algorithms

The value of the objective function in Equation (4.12) against the number of iterations for this scenario is plotted in , respectively, which corresponds to an assignment with sufficiently low complexity. Compared to the initial costs, the global cost in terms of these two networks dropped by nearly 97.5% and 94.7%. It can be seen that there are some plateau periods followed by sharp decreases during the second optimization step. These cases happen when the flow of a certain path becomes zero, then the complexity of the network decreases sharply. Finally, with respect to the value of the cost function obtained by the classical DA in the first step, an improvement over 95% and 80% is obtained by the DA with the RO strategy in the second step for the two networks, respectively. This fact highlights the benefit of the two-step optimization method.

Flow allocation in network representation

The flow assigned on all paths Taken together with G 1 and G 2 , compared to the initial flow assignment in Figure 4.10, their total link cost is dramatically decreased. If possible, the resolution algorithm diverts the path from complex areas. Some unnecessary or long paths may also be avoided. On the basis of the design of the link cost function, the link flow associated with the same node is balanced and optimally distributed to mitigate congestion and traffic complexity. Even if some links are allocated with massive flows, their costs can be quite low. In addition, since the cost of the paths is the sum of the cost of included links, the paths allocated with a large amount of flow usually have short distances. With such traffic assignment results, aircraft will make effective use of networks in urban airspace in order to minimize air traffic complexity and congestion. Besides, the number and the difficulty of potential conflicts will also be strongly mitigated.

{F i w | i ∈ [|1, R w |]} between OD

Conclusions

This Section presents a novel macroscopic air traffic assignment model to alleviate the congestion and traffic complexity of UAM route networks in future high-density urban airspace. Firstly, UAM corridors are designed which can be fitted in UAM route networks in urban airspace. An intrinsic air traffic complexity metric that analyzes the airspace geometry and traffic structure is defined based on the LDS. From a macroscopic perspective, an optimization problem under system equilibrium is formulated as an air traffic assignment model. The objective function involves both the airspace complexity metric and the congestion factors. The DA is a commonly used sequential decomposition algorithm to solve this optimization problem. On that basis and in consideration of the complicated structure of the link cost function, an efficient two-step strategy is proposed to solve the associated optimization problem. To demonstrate the applicability and effectiveness of the proposed algorithm, case studies are carried out on UAM route networks in urban airspace. The global costs of these two networks are respectively reduced by 97.5% and 94.7% compared to the initial flow allocation. Complex areas and unnecessary and long paths can be avoided. The flow is balanced and optimally distributed.

In summary, this research will shed light on mitigating the traffic complexity and congestion of urban airspace based on the UAM route network representation from a macroscopic perspective. The proposed model is able to assist regulators and ANSP for air traffic assignment in the strategic planning of UAM operations. A journal paper related to this chapter has been published [START_REF] Wang | Air Traffic Assignment for Intensive Urban Air Mobility Operations[END_REF].

Nevertheless, the segregation of 2D network layers, one-way UAM corridors, and the absence of bounds in terms of energy consumption assumed in this chapter have several limitations. The next chapters will address and improve the modeling and optimization encountered in more realistic UAM operations. Chapter 5

Static air traffic assignment model in 3D UAM route network with connected layers

Problem description

We develop in this chapter a complexity optimal air traffic assignment model to organize intensive UAM traffic in a multi-layer two-way UAM route network. Unlike segregated 2D networks, the horizontal layers are connected by vertical links in the multi-layer structure. It makes better use of limited urban airspace by providing more travel options, which has a more significant potential to reduce airspace complexity and congestion by air traffic assignment. The model is macroscopic and aggregates individual vehicle dynamics to describe the overall stream features in the network, such as flow structure, congestion, and efficiency. It is one of the first attempts to conduct air traffic assignments for 3D UAM operations from a macroscopic perspective.

3D UAM route network modeling

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

Graph representation of UAM route network

A UAM route network G consists of nodes V and links E. The nodes represent the waypoints, droneports, and delivery stations. The position of a node v is given as (x v , y v , z v ). A link e ∈ E represents the airway that joins a pair of nodes. The links E include horizontal links E h , upward links E u , and downward links E d . All links are considered to be associated with another link in the opposite direction. 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.

Inspired by (FAA, 2020a; Geister, 2017), volume segments are used to model links (FAA, 2020a; Geister, 2017). It is a performance-based 3D route segment consisting of 86 SATA in 3D UAM route network with connected layers The proportion of flow assigned to i-th path for OD pair w corridors that support one-way or two-way point-to-point UAM operations. All UAM vehicles within the volume segment follow specific rules, procedures, and performance requirements. Figure 1a illustrates an example of 3D volume segments for two horizontal links in opposite directions. The vertical view is shown in Figure 1b. For vertical links in opposite directions, a pair of upward and downward links are grouped through the same structure. The couples of 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. All grid segments have the same length L h , except that the grid length in the middle is in [L h , 2L h [. This is because the length varies for each link, and the grid cannot be divided equally. The height of the grid segment, length L h , and width L w should be determined according to the size of the protected zone of aerial vehicles.

A set of density points ρ e is 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 Figure 1b. Each density point g ∈ ρ e is characterized by position (x g , y g , z g ), heading θ g if e ∈ E h , speed s g , and density ω g . The position is determined based on the position of the 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 ω g , which is defined as follows:

ω g = f e L g N s s e (5.1)
where f e is the link flow. Note that with respect to Equation 4.7:

• The lateral flow sharing coefficient here is 1/N s . Unlike the previous formulation in Section 4.3 that for each link there is a central lane and N s lanes in the same direction, the UAM corridor in this formulation simply specifies N s lanes for each link.

• The conversion parameter β is characterized as a grid length L g . Generally L g = L h and it corresponds to a constant value β. However, for the density point in the middle of the link, L g ̸ = L h .

Feasible paths generation

The proposed method relies on an intelligent choice of the set of paths P w for each OD pair w. For each OD pair w ∈ W, the feasible paths P w in a bidirectional network are generated without cycles based on depth-first search (Migliore, Martorana, and Sciortino, 1990). Furthermore, the range and endurance of any given flight are constrained by Equation (3.1) and (3.2) formulated in Section 3.3.2.

Demand modeling

Delivery services are launched from the droneport to delivery stations, and the drone will also return to the droneport for recharging and maintenance purposes. Based on this assumption, both a droneport and a delivery station can be considered as the origin or destination. Each pair of OD w ∈ W is given with an associated demand d w . This traffic demand will be distributed among paths in the pre-computed path set 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 not necessarily balanced.

Mathematical model

Given a complete description of the proposed or existing UAM route network, air traffic assignment aims to allocate the demand of a given set of OD pairs minimizing the complexity of the air transportation system while keeping flight efficiency and operational efficiency within acceptable limits.

Decision variables

The decision variables are the proportions of flow γ = (γ w ) w∈W assigned to paths for each OD pair w ∈ W, where γ w = (γ 1 w , ..., γ Rw w ).

Constraints

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

γ i w ≥ 0, i ∈ {1, . . . , R w }, w ∈ W (5.2) Rw i=1 γ i w = 1, w ∈ W (5.3)
The flow of i-th path F i w is derived from the demand and the proportion of path flow:

F i w = d w γ i w , i = 1, ..., R w , w ∈ W, (5.4)
In addition, the air traffic flow should be meaningful for the assignment. If the flow on a path is positive, it should be greater than a path flow threshold F t :

F i w ∈ {0} ∪ [F t , d w ], i = 1, ..., R w , w ∈ W (5.5)
For consistency, the following assumption holds:

F t < min w d w (5.6)
Then, link flow f e can be calculated from path flow in the same way for 2D networks as equation 4.6:

f e = w∈W Rw i=1 δ e,w,i F i w , e ∈ E, i = 1, ..., R w , w ∈ W, (5.7)
To model the operational efficiency, air traffic flows on all links are assumed to have the same characteristics. In a multi-level network, there are three flight phases: climb, descend, and cruise. Therefore, an energy consumption coefficient α e is associated with each link according to the flight phase. Then, α e is introduced to measure the operational efficiency on link e in terms of distance and flight phases: (5.8) where φ h , φ u , and φ d are efficiency coefficients per distance unit depending on whether the link is horizontal, upward, or downward.

α e =      L e φ h , e ∈ E h L e φ u , e ∈ E u L e φ d , e ∈ E d
Based on α e , the energy consumption E for flights through the UAM route network is defined as follows. It measures link congestion by involving link flow f e .

E = e∈E α e f e ( 5.9)
An efficiency threshold β Ê is set to ensure the network efficiency remains within an acceptable range:

E ≤ β Ê (5.10)
where β is a factor to describe the proportional increase and β > 1. Ê is the minimum energy consumption, that is, the minimum value of Equation ( 9) when the demand d w is completely assigned to the path with the highest efficiency for all OD pairs w ∈ W:

Ê = w∈W d w min p∈Pw   e∈p α e  
(5.11)

Objective function based on linear dynamical systems

The objective function aims to minimize the complexity of flow patterns in a UAM route network. Because volume segments organize the traffic in parallel, in 3D UAM route networks, the air traffic complexity mainly occurs at intersections. The interactions between the different amounts of path flow from different directions can be complicated. Therefore, the objective function is defined as the sum of node complexities:

C (f ) = v∈V C v (f ) (5.12)
where f is the complete flow pattern {f e |e ∈ E}.

The procedure for computing the node complexity c v is based on the LDS and can be divided into three steps:

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

The principle used to estimate coefficient matrices of the LDS is the same as Section 4.3.3. However, note that in this problem, the LDS is three-dimensional.

To improve the computation efficiency, we propose a WMMSE estimator to derive the coefficient matrix of a three-dimensional LDS. WMMSE solves the following optimization problem:

min A,b g∈ρ ω g ∥ ẋg -(Ax g + b)∥ 2 (5.13)
where ρ are divergent density points ρ d v or convergent density points ρ c v . ω g is the density at point g and x g = (x g , y g , z g ) ⊤ is the three-dimensional position vector of point g. The speed vector of point g depends on the phase of the associated link:

ẋg =      s g (cos(θ g ), sin(θ g ), 0) ⊤ , g attached to e ∈ E h (0, 0, s g ) ⊤ , g attached to e ∈ E u (0, 0, -s g ) ⊤ ,
g attached to e ∈ E d (5.14) where θ g is the heading at point g and only exists in the cruise phase, s g is the speed of flow on the associated link. The details for calculating the coefficient matrix A that minimizes the WMMSE are given in Appendix A.1 and A.2.

The proposition A.1 in Appendix A.3 enables us to partially characterize the qualitative behavior of the system by the eigenvalues of matrix A. Based on the aforementioned discussion and Proposition A.2, the dominant eigenvalue is then introduced as a representative of all eigenvalues to describe the stability of the whole LDS:

Definition 5.1 Given λ 1 , • • • , λ N as the eigenvalues of a matrix, if λ ∈ {λ 1 , • • • , λ N }, and
Re( λ) = max(Re(λ 1 ), • • • , Re(λ N )), then λ is the dominant eigenvalue.
The proof of Proposition A.2 is given as well in A.4.

To summarize, Figure 2 presents the framework for calculating the dominant eigenvalues λ c v and λ d v 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 ∈ V, the efficiency factor α e and flow f e of associated link e ∈ E m v is also taken into account, where m ∈ {c, d} is the operational mode, E c v denotes the set of links with head v, and E d v denotes the set of links with tail v. The node complexity is defined as follows:

C v (f ) = m∈{c,d} |Re(λ m v )| e∈E m v α e f e ( 5.15)
Finally, the optimization problem can be formulated as follows:

min γ C (f ) = v∈V C v (f ) s.t. E ≤ β Ê Rw i=1 γ i w = 1, w ∈ W F i w = d w γ i w , w ∈ W, i = 1, ..., R w w∈W Rw i=1 δ e,w,i F i w = f e , e ∈ E F i w ∈ {0} ∪ [F t , d w ], w ∈ W, i = 1, ..., R w γ i w ≥ 0, i = 1, • • • , R w , w ∈ W (P1)
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 [START_REF] Jahn | Optimal Routing of Traffic Flows with Length Restrictions in Networks with Congestion[END_REF]. We then seek efficient heuristic algorithms to produce good-quality solutions in reasonable computation time for this complex, high-dimensional, and NP-hard problem.

A two-phase algorithm for air traffic assignment problem

Though DA has been successfully applied to small-scale networks in Chapter 4, there are some limitations to applying DA to this problem. Due to the marginal cost calculation, determining the flow and the path of the flow transfer is computationally expensive. This makes it inefficient to apply DA alone to larger-scale networks. In 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.

Therefore, we propose a two-phase optimization algorithm for efficiently solving the problem (P1). The SA is used in the first phase to explore a good-quality solution as initial values for the next phase. Based on the DA, the second phase aims to further optimize the problem by exploiting the search space based on the initial solution generated by SA.

The resulting optimization architecture is shown in Figure 3. It is made up of three main modules: initialization, simulation, and optimization module. The simulation module is first built by constructing the UAM route 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 γ. Based on γ, the simulation module generates the associated flow pattern in the UAM route 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 γ.

Simulated annealing

Due to the high-dimensional continuous search space in the problem (P1), we use the SA algorithm to provide initial values for DA in order to improve its convergence performance.

Figure 4 illustrates the mechanism for each transition of SA. Based on the current state γ i , the neighbor state γ j is generated. In this research, we propose a probabilistic performance-based neighborhood generation strategy. The goal is to prevent DA from falling into local optima 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 the complexity of the nodes on a path:

C p (f ) = v∈p C v (f ), p ∈ P (5.16)
We aim to generate an initial flow pattern such that paths in each OD pair have similar complexities. To this end, for each OD pair w, a random proportion of flow δ 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, δ is set to a small value. The selection probability of the outflow path in each OD pair w ∈ W is given by:

P out w,i =          Cp(f )- min p∈Pw ,Fp>0 Cp(f ) p∈Pw ,Fp>0 Cp(f )- min p∈Pw ,Fp>0 Cp(f ) , if F p > 0 0, Otherwise (5.17)
and the selection probability of the inflow path is given as:

P in w,i = 1/ max(C p (f ), ε) -min p∈Pw (1/ max(C p (f ), ε)) p∈Pw 1/ max(C p (f ), ε) -min p∈Pw (1/ max(C p (f ), ε)) (5.18)
The Metropolis principle, described in Section 2.3.2.2, is used to determine whether a neighbor state can be accepted. If yes, the state is updated by its neighbor in terms of the current OD pair. If the neighbor is rejected, only the variables of the current OD pair need to be recovered. When the stopping criterion of SA is met, the final state is used as the initial flow pattern for DA in the optimization module.

Algorithm 5 Probabilistic performance-based neighborhood generation strategy

Input: γ: Current state Output: γ: Generated neighborhood of γ 1: procedure GenerateNeighbor(γ)

2:
for w ∈ W do 3:

do 4: Calculate P out w,i , P in w,i , i = 1, • • • , R w 5:
r 1 , r 2 ← uniform(0, 1)

6: i 1 , i 2 ← 0 7: do 8: r 1 ← r 1 -P out w,i 1 9: i 1 ← i 1 + 1 10:
while (r 1 > 0)

11:

do 12: r 2 ← r 2 -P in w,i 2 13: i 2 ← i 2 + 1 14: while (r 2 > 0) 15: δ ← min uniform(0, 5F t /d w ), γ i 1 w 16: γi w ← γ i w , i ∈ {1, • • • , R w } \ {i 1 , i 2 } 17: γi 1 w ← γ i 1 w -δ 18: γi 2 w ← γ i 2 w + δ 19:
while (constraint in problem (P1) is not satisfied)

Dafermos' algorithm

Given a link e = (v, v ′ ), the associated link cost function C e defined for the DA is given as follows:

C e (f ) = |Re(λ d v )| + |Re(λ c v ′ )| α e f e ( 5.19)
Considering the link cost function in Equation 19, applying the DA requires overcoming the same difficulties as those presented in Section 4.4. Similarly, the marginal cost is computed by the numerical differentiation approach described in Equation (4.13). Concerning the optimization of flow transfer in Equation (2.13), for the efficiency of the algorithm, the constraint interval of σ w is uniformly discretized into a set of size N d . The element σ 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.

96

SATA in 3D UAM route network with connected layers

Results

Scenario definition

In this chapter, the case study focuses on the urban airspace of 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 UAM route network is constructed as an example network covering most Singapore areas for delivery purposes.

The map view of the resulting network is plotted in Figure 5. The graph representation of the transport network is shown in Figure 6a, 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 15m 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. 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. Figures 7a and7b plot the demand imposed on all OD pairs in these two cases. (5.20) where rand p generates a random number in [0, 1] for path p.

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 AON assignment, SA, DA, MSA, and DA initialized by MSA. The detail of these algorithms can be found in Section 2.3.2.2. 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 equipped with Intel ® i7-8750H CPU and 32GB DDR4 RAM.

To quantify the computational efficiency and performance of different models, several performance indicators are computed. The main indicators are computation time and total complexity defined in Equation [START_REF] Bruno | A model for aircraft evaluation to support strategic decisions[END_REF]. Besides, the flight efficiency in the UAM route network is also measured to ensure that it is at an acceptable level. In addition to energy consumption defined in Equation ( 9 where L e is the length of link e, |W| is the cardinality of W.

Finally, the energy consumption, as defined in Equation ( 9), is also reported. Table 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 the sample mean plus/minus margin of error. The best results are marked in bold. Among all algorithms, the AON assignment has the fastest computation time, but the complexity has barely decreased, which is roughly 24% in both cases. The 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. The DA is the most timeconsuming algorithm. However, if the DA is initialized by SA, which is the proposed As can be seen, the AON assignment outperforms other methods in terms of flight efficiency and energy consumption 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 objective 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 Figure 8. 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 the 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 Figure 9 and 10. The convergence curve of SA is illustrated in Figures 9a and10a. 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 Figures 9b and10b The evolution of complexity in terms of DA with random initialization is shown in Figure 11. 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 Figure 8a and 8b. In addition, a large number of iterations is required for DA without appropriate initial values. This fact can also be observed in Figure 12, 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 Table 3 and4

Path flow assignment

In order to better interpret the results of the proposed model, the related results consistent with the experiments in 8a and 8b are discussed in the UAM route network representation. To start with, the path flow allocation for the two cases is set out in Figures 13 and14. 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 Figures 13a and14a, the flow pattern assigned randomly brings high complexity on almost all paths. The flow optimized by SA is allocated to fewer paths but with better quality in each OD pair, which is illustrated in Figures 13b and14b 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 (Figure 13b). Refined by DA in the second phase of the proposed model, the final path flow allocation result is presented in Figures 13c and14c. The path complexity is 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 Figure 13c). For each OD pair, the flow in different directions is assigned to a different path, and the path complexity is much lower.

Link flow allocation

In the network representation, the link flow evolution is shown in Figures 15 and16. 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 in managing traffic on waypoints. Optimized by our model, the link flow allocation has been significantly changed, as shown in 15b and 16b. The main flow is allocated to vertical links and some inner links of each network layer. A few flows are found on the 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.

Node complexity

Under the previously mentioned flow allocation result, Figures 17 and18 represent the evolution of node complexity in the two cases. The color is related to the node complexity. The initial complexity of nodes is plotted in Figures 17a and18a. 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 Figures 15b and16b, 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.

Conclusions

In this chapter, we introduced an efficient methodology to address the 3D SATA 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 introduced to model the operational density along each link. Cylindrical airspace is designed to capture the complexity at each node. The 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 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 UAM route network over the Singapore urban airspace given two different traffic demands. In addition, baseline methods including AON assignment, DA, SA, MSA, and DA initialized by MSA, are 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 UAM route 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 UAM route network are limited in Singapore. With the increasing density of UAM operations, these findings may assist the ANSP in 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.

A journal paper related to this chapter has been published (Wang et al., 2022b).

Chapter 6

Dynamic air traffic assignment model for high-density UAM operations

Problem description

This Chapter mainly focuses on the dynamic aspect of air traffic assignment.The problem is formulated as a Dynamic Air Traffic Assignment (DATA) model for highdensity UAM operations, which aims to allocate the UAM traffic flows among air routes in the planning horizon to organize UAM traffic flows and reduce airspace complexity within a centralized control scheme while meeting the demand and respecting some criteria.

The definitions of the variables and notations used in this Chapter are listed in Table 6.1.

UAM route network modeling

The planning horizon is discretized into K equal time intervals. The time-varying demand D w,k refers to the assigned flow on OD pair w at time interval k, which can be obtained in advance from Providers of Services for UAM (PSU) and Supplemental Data Service Provider (SDSP) services (FAA, 2020a).

The UAM route network is modeled as a directed graph in a similar way as section 2. In addition, the capacity of traffic flow on links is defined as f c e . The volume segment U e described in section 2.1 is used to model each link e in the route network as 3D blocks of airspace. Each volume segment is composed of several air lanes, which are further divided into grids representing the associated flow information. The center points of grids on a volume segment U e are denoted as density points g e , which is illustrated in Figure 1. Each point g ∈ g e is characterized by static properties including coordinate (x g , y g , z g ), the direction of velocity θ g if it belongs to a horizontal link, and dynamic properties in time interval k including average flow speed s g,k and density ρ g,k . Thereinto, the θ g and s g,k correspond with the traffic flow traversing the grid.

In this study, the amount of air traffic flow distributed per grid segment g during time interval k is defined as the traffic density ρ g,k :

ρ g,k = Fg,k s g,k (6.1)
where Fg,k is the accumulative path flow on point g and is given later in Equation (6.11).

Mathematical model

Decision variables and constraints

The decision variables in this problem are the proportion of path flow in OD pairs during the planning horizon. Unless otherwise noted, path flow refers to the flow that passes through the origin of a path. Rather than link-based models, the benefit of this formulation is that it can provide accurate path-based flow allocation results, which can facilitate the decision-making process. Such decision variables can be formulated as follows:

γ = {γ w,k |w ∈ W, k = 1, . . . , K} (6.2)
where γ w,k = {γ p,k |p ∈ P w }, w ∈ W, k = 1, . . . , K (6.3)

In addition, the decision variables are non-negative and satisfy the conservation law:

γ p,k ≥ 0, p ∈ P w , w ∈ W, k = 1, . . . , K (6.4) p∈Pw γ p,k = 1, w ∈ W, k = 1, . . . , K (6.5)
In addition, the non-zero path flow should be greater than a threshold to make it meaningful for the assignment:

F p,k ∈ {0} ∪ [F t , D w,k ], p ∈ P w , w ∈ W, k = 1, . . . , K (6.6)
where the path flow can be derived as:

F p,k = D w,k γ p,k , p ∈ P w , w ∈ W, k = 1, . . . , K (6.7)
Unlike static formulation, the flow may not traverse a link on a path all the time during a given time interval. To model the propagation of traffic flow, we first introduce the travel time at free-flow conditions from the origin of a feasible path p to a point g on this path:

τ p,g = e∈pg L e s e + d(v f g , g) s eg , g ∈ p, p ∈ P w , w ∈ W (6.8)
where p g is the truncated route segments from the origin of path p to point g, v f g is the final node of p g , L e is the length of link e, s e is the free flow speed on link e, and e g is the link containing point g. where T d is the common duration for all time intervals. The average path flow for a path p on a certain point g during the k-th time interval can be calculated as the weighted average of flows in the previous N -th and (N + 1)-th time interval uniformly distributed in each corridor, where the weight is linked to the time gap τ ε :

F g,p,k = 1 N s (1 - τ ε T d )F p,k-N + τ ε T d F p,k-N -1 , p ∈ P w , w ∈ W, k = 1, . . . , K (6.10)
This rationale can be interpreted in Figure 6.1b, a case τ p,g < T d . The path flow in time intervals [t 0 , τ p,g ], [t 1 , τ p,g ], and [t 2 , τ p,g ] are flows from previous time steps that cannot be finished. Thus the average path flow on a point g in each time interval has to be calculated across two timesteps in that case.

The accumulative path flows on a point g except for nodes in a k-th time interval can be derived as: Due to the limited traffic flow capacity f c e for each link, the link inflow and outflow are constrained so that the average flow speed s e,k can be specified as a reasonable value independent of the traffic density:

Fg,k = w∈W p∈Pw F g,p,k δ eg,p , g ∈ G \ V, k = 1, . . . ,
f in e,k ≤ f c e , e ∈ E, k = 1, . . . , K (6.14) f out e,k ≤ f c e , e ∈ E, k = 1, . . . , K (6.15)
where f c e is the traffic flow capacity for link e.

Objectives

Firstly, we introduce the congestion cost G e,k for each link e in time interval k. In addition to evaluating the traffic density of incoming and outgoing flow, the energy consumption of the air traffic flow on each link is taken into account. G e,k is defined as follows:

G e,k = R∈{in,out} aα e f R e,k + bf R e,k 2 , e ∈ E, k = 1, • • • , K (6.16)
where a, b are coefficients, and α e is the energy consumption factor of link e defined by equation [START_REF] Barnhart | Demand and capacity management in air transportation[END_REF].

The linear and quadratic form of link flow represents the energy of displacement and the energy to maintain the free flow speed, respectively. The purpose of coefficients a and b is to pre-weighting these terms. According to [START_REF] Delahaye | Optimisation de la sectorisation de l'espace aérien par algorithmes génétiques[END_REF], a feasible selection of these coefficients is a = 1/|ᾱ fk |, b = 1/(4 f 2 k ), where ᾱ = e∈E α e /|E|, fk is the average flow of all links in time interval k:

fk = e∈E f in e,k + f out e,k 2|E| (6.17)
However, as mentioned previously, the congestion only reveals partial information about the complexity of UAM traffic. For example, Figure A.3 in Appendix A.5 gives four traffic scenarios with the same traffic density but with very different complexity in managing the traffic. In order to assess the interdependency between congestion and the difficulty of managing air traffic, the objective function takes into account the air traffic complexity cost and congestion cost.

Before modeling the complexity cost, we first define cylindrical airspace A v around each node v. One way to determine the size of A v is to include the first N r rows of the grid closest to node v in the horizontal and vertical direction. The density points involved in the cylinder are the observations to measure the air traffic complexity. Figure 6.2 provides an example of the cylindrical airspace around a node, which is connected by several one-way or two-way links. The density points involved in the cylinder are used to measure the air traffic complexity around the node.

The complexity is computed at each node in cylindrical airspace A v in the same way as section 3.3, except that it is time-dependent.

To fit the observations of the aircraft in each cylindrical airspace A v around node v with the LDS, we propose a WMMSE criterion. Traffic density ρ g,k is added in this formulation to indicate the importance of traffic flow on each grid g to the LDS over the time interval k. The WMMSE is given by: where x g = (x g , y g , z g ) ⊤ is the position vector of point g. The speed vector of a point g depends on the type of associated link e g :

ẋg,k =      s g,k (cos(θ g ), sin(θ g ), 0) ⊤ , e g ∈ E h (0, 0, s g,k ) ⊤ , e g ∈ E u (0, 0, -s g,k ) ⊤ , e g ∈ E d (6.19)
where s g,k is the average flow speed on e g in time interval k, and θ g is the heading of flow.

The negative real part of the eigenvalues λ of the coefficient matrix Â, which is the solution of equation ( 13), corresponds to a contraction mode of the LDS. It introduces a high level of control workloads as trajectory points are converging [START_REF] Delahaye | Air traffic complexity map based on non linear dynamical systems[END_REF]. In addition, the larger the eigenvalue, the more complex the traffic conditions.

We then define the complexity cost X v,k for each node v over a time period k as the sum of the absolute value of negative real parts of λ of the LDS in the cylindrical airspace A v :

X v,k = 3 i=1 | min(Re(λ i ), 0)| (6.20)
The objective function C is defined as the weighted sum of the air traffic complexity and congestion of UAM operations:

C = K k=1 (1 -ϕ) v∈V X v,k + ϕ e∈E G e,k (6.21)
where ϕ ∈ [0, 1] is the weighting parameter to control the cost of each group. Since these two costs measure the disorder of air traffic from different perspectives, they can be seen as independent of each other. In this study, the parameter ϕ is selected to balance the magnitude of the total complexity cost and congestion cost to make them equally important.

Simulation-based rolling horizon framework

Conventional traffic assignment algorithms assume the cost function to be convex and separable with respect to link flow, which is not always true in real operations. For example, in some complex traffic situations, such as the intersection of several links with dense traffic flow, the cost function can be modeled to evaluate the air traffic complexity in a non-separable way, such as equation (6.20). In addition, the cost function may also be non-convex in some cases, even without analytical form.

C p,k = (1 -ϕ) v∈p X v,k + ϕ e∈p G e,k , p ∈ P w , w ∈ W, k = 1, • • • , K (6.22)
Then, the AON assignment is performed on paths in P w :

γ p,k =    1, p = arg min p∈Pw C p,k 0, otherwise (6.23)
This initialization strategy is adapted to the objective function. It is able to provide an initial path flow pattern with relatively low air traffic complexity and congestion, which will improve the efficiency of the optimization algorithm.

Two main modules are contained in this framework: the simulation module and the optimization module. The simulation module includes the modeling of a 3D route network under UAM representation. Together with the network dynamics, it supports calculating the characteristics of the route network and the objective function C . In the optimization module, a rolling horizon approach is introduced to address the DATA problem in the planning horizon by decomposing the original problem into several subproblems. This approach results in substantial savings in computational time while preserving the relationship between the different time intervals. Each sub-problem k can be formulated as a modified SATA problem, with the objective function as

C k = (1 -ϕ) v∈V X v,k + ϕ e∈E G e,k , k = 1, . . . , K (6.24)
and constraints defined by equations (6.4)-(6.15). In order to evaluate C k , previous time intervals {k -1, • • • , k -N max } should be considered, where N max = ⌈max g∈p,p∈Pw,w∈W τ p,g /T d ⌉. Indeed, due to the residual demand, except for the first time interval, the flow pattern in each time interval is determined based on the results of N max previous time intervals. Figure 6.4 gives the scheme of the rolling horizon, corresponding to the case in which the maximum time gap is less than the time interval T d , that is, N max = 1. In this study, the time length of each rolling horizon and the time shift are both set to T d . The start time and the end time of each rolling horizon are aligned with each time interval. Sliding along the time axis, the rolling horizon begins at t 0 , and the optimization algorithm is applied in the time interval [t 0 , t 1 ]. Next, the rolling horizon is shifted by T d , and the current interval for optimization becomes [t 1 , t 2 ]. The receding process is repeated until the final rolling horizon [t K-1 , t K ] is reached. Since the maximum time gap in this example is less than the duration of the time interval, in Figure 6 The optimization module, as shown in Figure 6.3, directly controls and updates the decision variables γ in a nested loop over OD pairs and sub-problems. This formulation strengthens the connection between sub-problems by collaboratively adjusting the decision variables involved in each rolling horizon. A simulated annealing algorithm using parallel computing and a novel neighborhood generation strategy is proposed as the optimization algorithm in the optimization module, which will be introduced in the next section. Thus, based on the current decision variables, the simulation module generates the flow pattern for calculating the value of objective function C , which guides the optimization module to approximate the optimal DATA solution with low air traffic congestion and complexity. Finally, each outer loop performs a convergence test. If the stopping criterion is reached, the whole framework will be terminated.

The main advantages of this framework can be listed as follows:

• The route network has already been designed in Chapter 3. A feasible path set is created in advance for each OD. Each feasible path set includes K-shortest paths that minimize operational cost, efficiency cost, and safety cost of the route network. This avoids the heavy computation of searching shortest paths and extra cost function evaluation for large-scale route networks.

• The rolling horizon approach in the optimization module can reduce the problem size and computational burden of cost function evaluation. The decision variables involved in each rolling horizon are collaboratively updated to strengthen the connection between sub-problems. Moreover, it also enables a refined optimization by multiple threads and placed in a completed task queue (lines 5-7). A thread pool is used in this process to avoid latency. The congestion costs are not involved in parallel computing as they are straightforward to compute. The change in the objective function value of neighboring solutions is computed based on the completed task queue (lines 8). The Metropolis criterion, which can be referred to in Section 2.3.2.2, is used as the acceptance criterion. If the generated neighboring solution is accepted, the current solution and costs will be updated in terms of the impacted paths (lines 9-13). Also note that, after the cost evaluation and update for each time interval, the complexity cost for the next time interval is required to be completely computed to update with the modifications of the flow pattern in the last time interval. Algorithm 6 Neighborhood generation method for DATA problem.

σ ← uniform(0, γ p 1 ,k ), γ p 1 ,k ≥ Ft D w,k γ p 1 ,k , γ p 1 ,k < Ft D w,k 7: γp 1 ,k ← γ p 1 ,k -σ 8: γp 2 ,k ← γ p 2 ,k + σ 9:
while (Constraint (6.15) is not satisfied) 10:

return {p 1 , p 2 }, γw,k 11: end procedure PSA is stopped if any of the following stopping criteria are met:

1. ( Ĉ -Ĉc )/ Ĉ < 0.1%, if the best objective value Ĉc of current temperature is improved less than 0.1% than previous one Ĉc ;

2. T c < 10 -4 T init , if the current temperature is less than 10 -4 of initial temperature.

Experiments

Scenario definition

We create a parcel delivery scenario conducted by intensive UAM operations within a centralized control scheme in Singapore's urban airspace. The associated UAM route network is shown in Figure 6.5b. It was designed by the methodology presented in Chapter 3. Some basic topology features of the UAM route network are summarized in Table 6.2 for each layer and the network, including the number of nodes, the number of links, Average Degree (AD), Average Clustering Coefficient (ACC), and Degree Of Assortativity (DOA). In addition, in this network, there are 7 vertiports and 12 delivery 120 Dynamic air traffic assignment model for high-density UAM operations Algorithm 7 Incremental parallel cost calculation and update.

1: procedure CostEvaluation(γ w,k , X, G) 2: T Q ← ∅ ▷ Task queue 3: p, γw,k ← NeighborhoodGeneration(γ w,k ) 4: T Q ← v∈ p X v,k (γ w,k ) 5: do 6:
T Q comp ← Invoke(T Q) ▷ Concurrent execution by multiple threads of tasks submitted to T Q, until all tasks have been completed

7: while (T Q ̸ = ∅) 8: ∆C w,k ← (1 -ϕ) v∈ p (X v,k -X v,k (γ w,k )) + ϕ e∈ p (G e,k -G e,k (γ w,k )) 9: if Accept(∆C w,k ) then 10: γ w,k ← γw,k 11: X v,k ← X v,k (γ w,k ), v ∈ p 12: G e,k ← G e,k (γ w,k ), e ∈ p 13: end if 14:
return X, G, γ w,k 15: end procedure sites distributed among the city, with a total of 84 OD pairs. Their connections and geographic locations are visualized in Figure 6.5a. According to hourly demands given in advance, a certain type of UAVs make deliveries from vertiports to delivery sites in this UAM route network. In the presence of dense obstructions at low altitudes, this UAM route network has limited route options, making it much more difficult to optimally assign the UAM traffic flow. The parameter settings of the proposed simulation-based framework are summarized in Table 6.3. The time-varying traffic demands (Figure 6.6) are focused on six peak hours of the day, from 09:00-10:00 to 14:00-15:00. For future high-density UAM operations, these demands are simulated based on the expected hourly number of parcel delivery UAVs in metropolitan areas predicted by Doole, Ellerbroek, and Hoekstra (2020) and the To evaluate the effectiveness and performance of the proposed PSA, we also conduct a comparison study with representative conventional DTA algorithms for this DATA problem. We compare PSA with commonly used traffic assignment algorithms, including DA, MSA, and PM. These algorithms are also integrated into the simulation-based framework in Figure 6.3. The details of these comparative algorithms have been presented in section 2.3.2. The stopping criteria for these algorithms are in accordance with PSA. Note that, due to the limitations of traditional traffic assignment algorithms for this DATA problem in terms of network size, traffic volume, and mathematical formulation, we only select representative algorithms for comparison with our proposed algorithm to demonstrate the performance of PSA on UAM route networks with highdensity traffic volume, high complexity, and congestion. All algorithms are implemented in Java on a laptop equipped with Intel ® i7-12700H CPU and 16GB DDR5 RAM. The processor has 20 threads and 14 cores in total, including 6 performance cores for heavy tasks and 8 efficient cores for background tasks.

Performance analysis of PSA algorithm

Parallelization performance

To evaluate the computational time improvement of the proposed acceleration strategy in terms of SA, we conduct the experiment on this scenario by increasing the number of threads from 1 to 20. It is noteworthy that the single-thread case is similar to a sequential process. For each number of threads, we record 10,000 times the CPU time for Figure 6.7 shows the boxplots of the computational time of a transition in PSA for different numbers of threads. It can be noted that the performance is not significantly improved when more than 8 threads are used. This limit can be related to two factors: the number of cores in the processor and the average number of nodes per path. The results indicate that the proposed parallel strategy can result in a speedup of nearly three times over the sequential approach by using a sufficient number of threads. The CPU time per transition of PSA using 20 threads can be reduced to (9729 ± 2.5) × 10 -5 s in 95% confidence interval. It is quite efficient since each transition includes the cost evaluation and update for all OD pairs and all time intervals, which is K|W| = 504 times in this studied case. Averaged over each operation, Algorithm 7 only takes less than 0.2ms.

Comparison with conventional traffic assignment algorithms

In addition to the objective function, the flight efficiency in the UAM route network is also measured to ensure that it is at an acceptable level. Two indicators are introduced, including average travel time t: For all algorithms, We performed 100 experiments with different random seeds, where the random seeds control the initial values and the inherent randomness of the algorithm. Table 6.4 illustrates the experimental results of all aforementioned algorithms. 95% Confidence Interval (CI) for the standard normal distribution is used to report the results, and the format is the sample mean plus/minus margin of error. The best result in each column is underlined. In addition, Figure 6.8 visualizes the reduction of objective function values obtained by these algorithms. As the performance of DA and PM is close, we additionally compare the objective function values produced by these two algorithms. As for representative conventional DTA algorithms including DA and PM, there is no significant reduction in the value of objective function compared to PSA. This result may be explained by the fact that several assumptions for conventional DTA algorithms are not fully satisfied in the DATA formulation. In particular, the objective function does not have an analytical form and is not guaranteed to be convex, which may make traditional DTA algorithms easily fall into local optima. Another important point lies in the scale of the problem. With high-density traffic volume and large network size, traditional DTA algorithms may not be applicable to this scenario. In terms of the computation time, PM is much lower than other algorithms. DA takes the longest computation time, but the objective function value of DA is better than PM. DA also produces the flow pattern with the best flight efficiency.

In view of the proposed model PSA, the value of the objective function is significantly reduced by (32.20 ± 0.29)%. In the objective function, the complexity cost has been reduced by almost half. There is a limited reduction in congestion of nearly 10%. A possible reason could be the limited path selection in the network and the low congestion level of the initial flow pattern. The computation time of PSA is also satisfying, considering the large size of the network, the high-density traffic volume, and the heavy computation involved in the LDS. The differences in average travel time and path length are quite small compared to the best results, which demonstrate that the flight efficiency of flow patterns obtained by PSA is not degraded. In summary, the proposed PSA algorithm is able to output the flow pattern with a significantly decreasing air traffic complexity and low congestion in a reasonable time without degrading the flight efficiency, which supports the primary objective of this study.

To further analyze the optimization results at the operational level, the path flow allocation, link flow allocation, air traffic complexity, and congestion in the UAM route network will be presented in the following sections.

Path flow allocation

To start with, the initial path flow allocations in each time interval are illustrated in Figure 6.9. Due to the high amount of paths in each OD pair, the initial path flow allocations are segmented into 4 subfigures. As can be observed in Figure 6.9, the demand is completely allocated to a specific path in each OD pair. Although the initial path flow pattern following the AON strategy has a low congestion level, there is still considerable potential to reduce air traffic complexity by better distributing traffic on various candidate paths. Additionally, in many OD pairs, the demand is assigned to the same path across all time intervals. The impact of the residual path flow in previous intervals on the current time interval requires the consideration of their interdependency. For each OD pair, it could be a preferable solution to have different flow allocation results for different time intervals. The path flow distribution optimized by PSA is represented in Figure 6.10. A significant difference compared to Figure 6.9 can be observed in that the demands for each OD and time interval can be assigned to multiple paths. In addition, most OD pairs have different path flow patterns across time intervals. This fact demonstrates that after being optimized by PSA, the candidate paths can be sufficiently and optimally utilized to mitigate congestion and air traffic complexity. The optimization results can be accurately and clearly provided to ATC service providers and aircraft operators for implementation in the form of the number of flows on each path in each time interval. The configuration is the same as Figure 6.9.

Link flow allocation

Another important illustration lies in the link flow allocation. Figure 6.11 presents the initial link flow allocation in each time interval in the network representation. According to all subfigures, the traffic flow is mainly distributed on the three lower layers and the vertical links connecting them. The two upper layers are not fully used, especially the top layer. In accordance with the initial path flow distribution result in Figure 6.9 that the demand is assigned to the same path in many OD pairs across all time intervals, the link flow distribution is similar for different time intervals.

The link flow distribution optimized by PSA is shown in Figure 6.12. Compared with the initial link flow pattern in Figure 6.11, the link flows are distributed in all layers. It is also worth noting that, in network layers with few links and low connectivity, such as the top layer and bottom layer, the links are assigned with relatively less flow. An explanation could be that UAM operations in these layers are more likely to introduce congestion and air traffic complexity. On the contrary, the link flows are densely concentrated in the middle layers and the vertical links connecting them. Since there are more route choices in these layers, a large number of link flows can be enabled to operate without bringing high congestion and air traffic complexity. Furthermore, given with time-varying demand for each time interval, the link flow patterns in these subfigures are slightly different with regard to the amount of link flow and link usage.

Air traffic complexity and congestion

To evaluate the air traffic complexity and congestion, we illustrate the initial and optimized complexity cost and congestion cost in the network representation in Figure 6.13 and Figure 6.14, respectively. The color of each node represents the complexity cost, and the width of each link represents the congestion cost. As can be seen in Figure 6.13, the initial flow pattern involves high complexity costs. Each layer has a large number of nodes with a high level of complexity. The congestion mainly occurs in vertical links, especially the vertical links connecting the bottom three layers. This fact also demonstrates that initial flow patterns are mainly assigned in the lower layers. After optimization by PSA, it can be seen from Figure 6.14 that the complexity cost has been obviously reduced to a low level, especially for nodes in the three middle layers. Even for the top layer and bottom layer that have limited route choices, PSA is still capable of mitigating the air traffic complexity. The complexity cost of almost all nodes has been decreased to below 2000. With respect to the congestion cost, unlike Figure per link from the initial flow pattern to the optimized flow pattern in Figure 6.15 and Figure 6.16, respectively. These figures are intended to give a better view of the level of cost reduction compared to Figure 6.13 and 6.14. The colorbar indicates the changes in cost. According to Figure 6.15, most of the changes were made to reduce the complexity for nodes at different time intervals. Especially for some high-complexity nodes, the decrease is significant. In Figure 6.16, although the air traffic congestion is mitigated on a large part of links in different time intervals, there are still some links with increased congestion. This fact is also consistent with Figure 6.14. Rather than concentrating on several links, the congestion is optimally distributed to all links in the network in order to reduce the overall congestion in the UAM route network.

Conclusions

To tackle the time-varying demand in terms of dense UAM operations in the near future, we formulate a macroscopic DATA model within a centralized strategic planning scheme. The focus of this model is to allocate and organize traffic flows on routes in the UAM route network within the planning horizon in order to alleviate congestion and reduce air traffic complexity. Firstly, the UAM route network is represented as a directed graph, with air routes modeled as volume segments; and vertiports, delivery stations, and waypoints modeled as nodes. The DATA problem is then formulated as an optimization problem that considers constrained capacity, traffic flow propagation, and time-varying demand. The objective function involves congestion cost based on energy consumption and traffic density, as well as complexity cost based on the LDS. A simulation-based rolling horizon framework is proposed to decompose the DATA problem into sub-problems in each time interval. To overcome the limitations that lie in traditional traffic assignment algorithms, we propose the PSA as a meta-heuristic method using parallel computing and a novel neighborhood generation strategy to solve this problem efficiently. Experiments are conducted on a previously designed large-scale UAM route network for a high-density parcel delivery scenario in Singapore's urban airspace. The results in terms of parallelization performance, algorithmic comparison, path flow allocation, link flow allocation, and complexity and congestion costs in network representation demonstrate that optimized by the proposed model, the congestion and air traffic complexity can be efficiently mitigated to a satisfying level. The paths are sufficiently utilized, and the flows are well allocated.

In addition to the theoretical significance of this research, the proposed framework can also assist or provide advisories to ATC responsibilities and ANSP for a series of problems, such as unmanned air traffic flow management, air traffic assignment, UAM traffic flow analysis, Urban airspace complexity measurement, and UAM route network evaluation.

A journal paper related to this Chapter has been submitted (Wang et al., 2022a).

Chapter 7

Conclusions and perspectives

Conclusion

In high-density UAM operations in the near future, reducing congestion and structural constraints are key challenges. Pioneering urban airspace design projects expect the air vehicles to fit into structured UAM corridor networks. However, most existing air transport networks are not capable of handling the increasing traffic demand, which is likely to cause congestion, traffic complexity, and safety issues. Additionally, a major challenge is integrating UAM with the ATM system in urban airspace, especially in class G airspace.

Our first contribution was to develop an approach to design the UAM route network in low-altitude urban airspace. The proposed approach is based on open-source data and can be applied to urban airspace around the world. We also contributed in the context of air traffic assignment, which is able to remedy the deficiencies in the existing transportation system by assigning estimated future demands on the existing UAM systems. We tackled the static air traffic assignment for high-density UAM operations in 2D and 3D UAM route networks. In addition, we proposed a dynamic air traffic assignment framework to handle large traffic demands in dense UAM operations. A complete strategic planning strategy can be achieved by integrating the UAM route network design with the traffic assignment model. The details of these contributions are given as follows.

Firstly, we addressed the problem of designing UAM route networks in low-altitude airspace to minimize noise impact and maximize the efficiency and safety of UAM operations. On the basis of the open-source data, the UAM network is designed as a grid-based multi-layer route network that supports two-way traffic. The topology features of the route network have been analyzed. To provide alternative travel options for UAM traffic flow, we searched for feasible routes between OD pairs by solving the KSPD problem that minimizes link costs in terms of noise impact, safety and efficiency. The resulting feasible routes can potentially reduce airspace complexity and can be used for air traffic assignments. In addition, the impact of different parameter settings for link costs on UAM services has been explored. Singapore's urban airspace has been selected for a case study to demonstrate the feasibility of this approach.

Next, we dealt with the air traffic assignment problem. In order to adapt the increasing demand to the current airspace capacity, three novel macroscopic traffic assignment models have been proposed to mitigate the air traffic complexity and congestion by or-

134

Conclusions and perspectives

ganizing the structure of air traffic flows. Firstly, at the macroscopic level, the UAM traffic has been modeled to operate in a UAM route network with respect to the dense volume. The UAM corridor or volume segment has been designed and fitted into a graph representation.

The first air traffic assignment model is referred to as the SATA. Based on 2D LDS, the model has been formulated for 2D segregated UAM route networks. A two-step resolution method based on Dafermos' algorithm has been introduced to efficiently solve this optimization problem. A case study has been carried out on a two-layer segregated UAM route network with intensive UAM operations. The results demonstrate that the proposed model can successfully mitigate urban airspace congestion and organize the UAM traffic into a low-complexity flow pattern.

We then proposed a SATA framework for 3D two-way UAM route networks. Based upon the 3D LDS, a novel complexity metric has been defined as the objective function, which takes into account dynamic flow structure, congestion, and operational efficiency. A two-phase approach combining SA and DA has been introduced to efficiently solve this problem. To validate the proposed model, a case study of a UAM route network in Singapore's urban airspace with two different demands has been conducted. Comparative studies were 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, significantly reducing the complexity of the 3D UAM route network.

In the last part of the thesis, we introduced a DATA model to deal with the dynamic air traffic assignment problem in order to manage large time-dependent traffic demand in dense UAM operations. At each time step, the air traffic complexity and congestion can be evaluated in terms of intrinsic characteristics of UAM traffic flow, link flow distribution, and energy consumption. An advanced macroscopic formulation of dynamic air traffic assignment has been introduced. Then, given time-varying traffic demand, we introduced a simulation-based rolling horizon approach to determine the optimal path flow pattern by decomposing the original problem into several sub-problems and dividing the planning horizon into several time intervals. A modified SATA problem can be formulated at each time interval. We developed a PSA algorithm as the optimization method, with a novel neighborhood generation strategy specifically designed for this problem. The validation of the proposed model through a scenario of parcel delivery service in Singapore's urban airspace has been provided. The UAM route network in the case study was designed according to the methodology described in Chapter 3. Comparison studies were conducted with conventional DTA optimization methods. The proposed approach can find good solutions to decrease the complexity and congestion in the UAM route network with high-density UAM traffic demands. PSA was suitable for finding a near-optimal solution in a short computational time that could be used for real-time deployment.

In summary, we conclude that the UAM route network design approach and air traffic 7.2. Perspectives 135 assignment model have a high potential to handle dense UAM operations in UAM route networks in the near future, resulting in an increase in airspace capacity and a decrease in airspace complexity. The framework may be useful for ANSP in strategic planning for UAM operations and urban airspace design, which supports a robust and efficient UTM system.

Perspectives

Several research directions can be envisioned in future work:

1. In terms of UAM route network design, an interesting direction is to estimate the density of operations and analyze the impact of bottlenecks using different network structures.

2. Another important extension of air traffic assignment could consist in considering mixed traffic, for instance, automated air taxis sharing the airspace with parcel delivery UAVs. For dynamic air traffic assignment, the uncertainty on delay from the origin to a point on the path can also be considered in the mathematical formulation.

3. In air traffic assignment, while reducing the air traffic complexity and congestion, other factors can also be incorporated into the objective function, such as equity, privacy, and other air traffic complexity metrics or indicators. In addition, in this thesis, we dealt with multi-objective optimization by considering the weighted sum or the multiplication of each objective according to the preference or importance of the concerned objective. A better approach for decision-makers could be selecting their preferred solution from the Pareto front according to the different needs of stakeholders.

3. 1 63 4. 1 83 5. 1 123 6. 8 125 6. 9 131 A. 1 159 A. 2

 16318311238125913111592 ALOS World 3D DSM data covering Singapore. . . . . . . . . . . . . . . . 3.2 Existing area limits for UAM flights in Singapore's airspace. . . . . . . . . 3.3 Total resident population of Singapore by subzone. . . . . . . . . . . . . . xi xii LIST OF FIGURES 3.4 Geofence with a keep-out distance of 30m. . . . . . . . . . . . . . . . . . . 51 3.5 The calculation of network construction area based on convex hull and the AMABB. The waypoints are determined using grid-based discretization. . 52 3.6 The distribution of obstacles with geofence in a network construction area at altitudes 150ft, 175ft, 200ft, 225ft, and 250ft. The color indicates the altitude of obstacles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53 3.7 Measurement of SPL from UAV to a location. . . . . . . . . . . . . . . . . 54 3.8 Notations for different types of links and subzones represented by the AMABB. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56 3.9 The resulting multi-layer network in Singapore urban airspace. . . . . . . 59 3.10 Normalized link cost for crossing hazardous airspace in the grid-based network. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 3.11 Normalized link cost of flight efficiency in the grid-based network. . . . . . 60 3.12 Normalized link cost of noise impact in the grid-based network. . . . . . . 61 3.13 Average path cost for different OD pairs with different weighting parameters. 62 3.14 The number of times the links appear in the k-shortest paths between OD pairs with different weighting parameters. . . . . . . . . . . . . . . . . . . Illustration of algorithm 3. . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 4.2 Illustration of density points on UAM corridors, disk complexity areas, and vector fields of the LDS. . . . . . . . . . . . . . . . . . . . . . . . . . 69 4.3 The process of an iteration of the modified DA. . . . . . . . . . . . . . . . 73 4.4 The two-step optimization algorithm for 2D SATA problem. . . . . . . . . 74 4.5 Illustration of random optimization strategy. A neighbor solution is generated, then its objective function value is compared with the previous one to determine if it can be accepted. . . . . . . . . . . . . . . . . . . . . 75 4.6 Representation of the two-layer network, in which the first layer contains 17 nodes, 28 links, and 2 OD pairs, and the second layer contains 24 nodes, 53 links, and 2 OD pairs. . . . . . . . . . . . . . . . . . . . . . . . 77 4.7 Evolution of global cost optimized by the proposed two-step optimization approach for two segregated networks. . . . . . . . . . . . . . . . . . . . . 78 4.8 Optimized allocation of path flows for two segregated networks. . . . . . .79 4.9 Optimized allocation of link flows for two segregated networks. . . . . . . 80 4.10 Representation of link costs and link flows under initial flow allocation in two segregated networks. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 81 4.11 Traffic assignment result in terms of flow allocation, and link cost for two segregated networks. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83 4.12 Traffic assignment result for links with the nonzero flow or nonzero link cost for two segregated networks. . . . . . . . . . . . . . . . . . . . . . . . An example of two-way UAM volume segments . . . . . . . . . . . . . . . 87 LIST OF FIGURES xiii 5.2 A example of cylindrical airspace, density points, and eigenvalue calculation of the associated LDS. . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 The optimization framework of the two-phase proposed algorithm on a simulation environment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.4 Probabilistic performance-based neighborhood generation and evaluation of objective function for each transition in SA. . . . . . . . . . . . . . . . 5.5 A map-based visualization of the UAM route network in Singapore's urban airspace. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.6 Graph representation of UAM route network used in this study. . . . . . . 5.7 The distribution of demand for OD pairs in two cases. . . . . . . . . . . . 5.8 Complexity evolution of the proposed two-phase algorithm. . . . . . . . . 5.9 Complexity evolution for two phases of the proposed algorithm in the first case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.10 Detailed complexity evolution for two phases of the proposed algorithm in the second case. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.11 Complexity evolution of DA without initialization phase in two cases. . . 5.12 Complexity evolution of DA initialized by MSA in two cases. . . . . . . . 5.13 Evolution of path flow and path complexity of the first case. . . . . . . . . 5.14 Evolution of path flow and path complexity of the second case. . . . . . . 5.15 Link flow evolution of the first case in the network representation. . . . . 5.16 Link flow evolution of the second case in the network representation. . . . 5.17 Node complexity evolution of the first case in the network representation. 5.18 Node complexity evolution of the second case in the network representation.106 6.1 Illustration of the travel time and the real-time path flow at point g on a path p connecting OD pair w in an example network. . . . . . . . . . . . . 6.2 Example of cylindrical airspace around a node to measure the complexity of air traffic. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Flowchart of a simulation-based framework for solving the DATA optimization problem. The simulation module includes the UAM representation of the route network, and the optimization module contains the optimization algorithm in a nested loop over OD pairs and rolling horizons.115 6.4 General strategy in the optimization module to approximate the optimal flow pattern of DATA using rolling horizon approach. This example includes six time intervals and the maximum time gap is less than the duration of the time interval. Inner loops are conducted for each OD pair in each time interval. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.5 3D UAM route network in Singapore's urban airspace. . . . . . . . . . . . 6.6 Simulated hourly UAM traffic demands for each OD pair from 09:00-10:00 to 13:00-14:00. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiv LIST OF FIGURES 6.7 Run time for a transition in PSA with the different number of threads. The median is indicated in green. The lower bound and upper bound of the box represent the first and third quartiles, respectively. The whiskers above and below the box include the values between the 5th and 95th percentiles. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . The reduction of objective function value after optimization of different algorithms compared to the initial state and the performance comparison between PM and DA. The boxplot is represented in the same way as Figure 6.7. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Initial path flows distribution in different time intervals. The x-axis represents the index of paths and the y-axis represents the index of time intervals from 09:00-10:00 to 14:00-15:00. The dashed white lines separate the OD pair. The color of each path at each time interval indicates the associated path flow. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 126 6.10 Distribution of path flows optimized by PSA in different time intervals. The configuration is the same as Figure 6.9. . . . . . . . . . . . . . . . . . 127 6.11 The initial link flow distribution in the UAM route network in different time intervals. The color represents the link flow in a logarithmic scale that covers a large range of values. . . . . . . . . . . . . . . . . . . . . . . 128 6.12 The resulting link flow distribution in the UAM route network in the different time intervals, optimized by PSA. The configuration is the same as Figure 6.11. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 6.13 The initial air traffic complexity and congestion are represented respectively on nodes and links of the UAM route network in different time intervals. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 129 6.14 The resulting air traffic complexity and congestion optimized by the proposed model are represented respectively on nodes and links of the UAM route network in different time intervals. . . . . . . . . . . . . . . . . . . . 130 6.15 Change in complexity cost in terms of each node compared with initial flow pattern and the flow pattern optimized by PSA. . . . . . . . . . . . 130 6.16 Change in congestion cost in terms of each link compared with initial flow pattern and the flow pattern optimized by PSA. . . . . . . . . . . . . . . Evolution of the LDS in terms of eigenvalues. . . . . . . . . . . . . . . . . Eigenvalue loci for four typical 2D traffic situations. The location of the observed traffic is indicated by the black dots. The blue arrows indicate the speed vectors. The arrowhead indicates the direction and the length of the arrow indicates the speed. . . . . . . . . . . . . . . . . . . . . . . . 160 LIST OF FIGURES xv A.3 Eigenvalue loci for four typical 3D traffic situations. The black dots represent the observed trajectory points of aircraft, the red arrows indicate their speed vectors, and the blue arrows indicate the speed vectors estimated by the LDS at equidistantly-partitioned points. The arrowhead indicates the direction and the length of the arrow stands for speed. . . . 160 List of Tables 1.1 Estimated UAM demands for low and high scenario (Anand et al., 2021).

3. 1

 1 Nomenclature specified to Chapter 3 . . . . . . . . . . . . . . . . . . . . . 3.2 Basic topology features of different altitude layers and the whole network 4.1 Nomenclature specified in Chapter 4 . . . . . . . . . . . . . . . . . . . . . 4.2 Parameters setting for route network modeling and optimization process in the experiment. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.1 Notations and definitions specified in Chapter ?? . . . . . . . . . . . . . . 5.2 Parameters setting . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.3 Performance comparison of models in terms of 95% CI for computation time, complexity, and flight efficiency in case 1. . . . . . . . . . . . . . . . 5.4 Performance comparison of models in terms of 95% CI for computation time, complexity, and flight efficiency in case 2. . . . . . . . . . . . . . . . 6.1 Nomenclature specified to Chapter 6. . . . . . . . . . . . . . . . . . . . . . 6.2 Basic topology features of different altitude layers and the UAM route network. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 6.3 Parameters setting in the experiment. . . . . . . . . . . . . . . . . . . . . 6.4 Performance comparison of PSA and representative conventional DTA algorithms in terms of computation time, objective function, and flight efficiency for a 95% CI. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xvii Introduction 1.

Figure 1 . 1 :

 11 Figure 1.1: Predicted growth of UAV in terms of commercial and recreation fleet.

Figure 1 . 2 :

 12 Figure 1.2: UAM maturity level (Goodrich and Theodore, 2021).

  conduct the TFM by managing conflicts among UAVs in high-density very low-level uncontrolled airspace. The conflicts are mitigated by assigning UAVs to different layers in designed multi-layer airspace. Chin et al. (2021) developed a strategic TFM technique for UAM operations with dynamic traffic demands by assigning each UAV to enter an airport or a sector during each time interval. The objective function takes delay and fairness into account. An extension of this work incorporates the demand with flights that have low file-ahead times (Chin et al., 2022).

Figure 1 . 3 :

 13 Figure 1.3: Current UAM challenges and solution.

2 . 3 . 4 .

 234 Efficient UAM representations are proposed to model air transport networks with flow-based UAM operations; Measurements of air traffic complexity are developed based on LDS involving the intrinsic characteristic of UAM traffic dynamics and flow distribution. The congestion is evaluated by traffic density and energy consumption; Several optimization methods are proposed to efficiently optimize different air traffic assignment problems within the simulation-based framework. A single simulation environment for each individual process is created to increase efficiency.

  Figure 2.1 represents an example of airspace classes in the US.

Figure 2 . 1 :

 21 Figure 2.1: ICAO Airspace classification (FAA, 2022).

Figure 2 . 2 :

 22 Figure 2.2: Urban airspace structure proposed in ConOps, FAA (FAA, 2020a).

Figure 2 . 3 :Figure 2 . 4 :

 2324 Figure 2.3: UAM corridor integrated in different airspace classes (FAA, 2020a).

Figure 2 . 5 :

 25 Figure 2.5: UAM airspace structure proposed by Amazon (Singireddy and Daim, 2018).

-

Figure 2 . 6 :

 26 Figure 2.6: New airspace structure proposed by SESAR Joint Undertaking (2020).

Figure 2 . 7 :

 27 Figure 2.7: Four potential airspace concepts.

Figure 2 . 9 :

 29 Figure 2.9: Three cases of traffic.

Figure 2 .

 2 10a, the aircraft are distributed homogeneously. On the contrary, the aircraft in Figure 2.10b aggregate and form a cluster. Although the proximity ◢ (a) Aircraft well distributed. (b) Aggregated aircraft.

Figure 2 . 10 :

 210 Figure 2.10: Two cases of aircraft in the sector.

Figure 2 . 11 :-Figure 2 . 12 :

 211212 Figure 2.11: Complexity heatmap in terms of proximity metric (Delahaye and Puechmorel, 2000).

Figure 2 . 13 :

 213 Figure 2.13: Complexity heatmap in terms of convergence metric (Delahaye et al., 2015).

  Delahaye and Puechmorel, 2000; Delahaye et al., 2002; Delahaye et al., 2004; Delahaye and Puechmorel, 2010; Treimuth et al., 2015; Juntama et al., 2022). With a strong theoretical framework, LDS evaluates the air traffic complexity by describing the evolution of a given traffic situation. It is able to efficiently quantify the level of disorder and interaction between a large number of aircraft.

  min σw C ( f ) s.t. fe = f e + (δ e,p w,min -δ e,pw,max )σ w , e ∈ E σ w ∈ [-F p w,min , F pw,max ]

Figure 2 . 14 :

 214 Figure 2.14: An iteration example of Dafermos' algorithm on a simple network.

Figure 2 . 16 :

 216 Figure 2.16: The procedure of each transition for a given temperature in the SA. Firstly, a neighboring solution is generated based on the current solution. After objective function evaluation, the acceptance criterion decides whether the new solution is accepted.If accepted, the current solution will be updated. Otherwise, the comeback operator returns to the previous solution and avoids duplication in memory.

  Delahaye and Puechmorel, 2000; Delahaye et al., 2002; Delahaye et al., 2004; Delahaye and Puechmorel, 2010; Treimuth et al., 2015; Delahaye et al., 2022). The models proposed in Chapter 4, ??, and 6 measure the air traffic complexity by this metric.

  .1.

Figure 3 . 1 :

 31 Figure 3.1: ALOS World 3D DSM data covering Singapore.

Figure 3 . 2 :

 32 Figure 3.2: Existing area limits for UAM flights in Singapore's airspace.

Figure 3 . 3 :

 33 Figure 3.3: Total resident population of Singapore by subzone.

Figure 3 . 4 :

 34 Figure 3.4: Geofence with a keep-out distance of 30m.

Figure 3 . 5 :

 35 Figure 3.5: The calculation of network construction area based on convex hull and the AMABB. The waypoints are determined using grid-based discretization.

Figure 3 . 6 :

 36 Figure 3.6: The distribution of obstacles with geofence in a network construction area at altitudes 150ft, 175ft, 200ft, 225ft, and 250ft. The color indicates the altitude of obstacles.

Figure 3 . 7 :

 37 Figure 3.7: Measurement of SPL from UAV to a location.

e

  f e (b) S f S (a)L dist(a,b) dadb (3.6) where dist(a, b) is the euclidean distance between points a and b, f S (a) is the Probability Density Function (PDF) of a point a in polygon S, f e (b) is the PDF of a point b in line segment e. To calculate L e,S , assuming that individuals are uniformly distributed in a polygon, and there is a uniform traffic flow of UAV on each link, then f S (a) = 1/|S| and f e (b) = 1/d e , where |S| is the area of polygon and d e is the length of line segment e. L e,S can be simplified as:

Figure 3 L h + 10 log 10 h 2 h 2 +

 32 (x -x 0 ) 2 + (y -y 0 ) 2 dhdxdy (3.8)If the link is horizontal (corresponds to Figure3.8a) and the x-coordinates of both endpoints are different, L e,S = max(x 0 ,x 1 ) min(x 0 ,x 1 )

Figure 3 . 8 :

 38 Figure 3.8: Notations for different types of links and subzones represented by the AM-ABB.

Algorithm 2 7 :

 27 Greedy algorithm for KSPD problem Input: Grid-based network: G, OD pairs: W, similarity function: sim(•, •), number of shortest path: k, and similarity threshold: τ p Output: k-shortest path for OD pairs: Ψ 1: procedure KSPD 2: for w ∈ W do 3: P ← Set of acyclic paths that connects OD pair w and satisfies the constraint (3.1-3.2) 4: p * ← arg min p∈ P C p 5: P ← P/{p * } 6: Ψ w ← {p * } while |Ψ w | < k and | P| > 0 do 8: p * ← arg min p∈ P C p 9: P ← P/{p * } 10: if ∀p ∈ Ψ w , sim(p, p * ) < τ p then 11: Ψ w ← Ψ w ∪ {p * }

Figure 3 . 10 :

 310 Figure 3.10: Normalized link cost for crossing hazardous airspace in the grid-based network.

Figure 3 . 11 :

 311 Figure 3.11: Normalized link cost of flight efficiency in the grid-based network.

Figure 3 . 12 :

 312 Figure 3.12: Normalized link cost of noise impact in the grid-based network.

  (a) Weighting parameters [1/3, 1/3, 1/3]. (b) Weighting parameters [2/3, 1/6, 1/6]. (c) Weighting parameters [1/6, 2/3, 1/6]. (d) Weighting parameters [1/6, 1/6, 2/3].

Figure 3 . 13 :

 313 Figure 3.13: Average path cost for different OD pairs with different weighting parameters.

  (a) Weighting parameters [1/3, 1/3, 1/3]. (b) Weighting parameters [2/3, 1/6, 1/6]. (c) Weighting parameters [1/6, 2/3, 1/6]. (d) Weighting parameters [1/6, 1/6, 2/3].

Figure 3 . 14 :

 314 Figure 3.14: The number of times the links appear in the k-shortest paths between OD pairs with different weighting parameters.

Figure 4 . 1 :

 41 Figure 4.1: Illustration of algorithm 3.

Figure 4 . 2 :

 42 Figure 4.2: Illustration of density points on UAM corridors, disk complexity areas, and vector fields of the LDS.

Figure A. 1

 1 in Appendix A.5.

Figure 4 . 3 :

 43 Figure 4.3: The process of an iteration of the modified DA.

Figure 4 . 4 .

 44 The first step corresponds to the classical DA and is related to operator E w,1 . The brute-force approach is used to solve the optimization problem related to flow transfer formulated in Equation (2.13). All values of σ w in the interval are evaluated, and the value that makes S minimum is chosen. The convergence speed is very fast until a local optimum is reached. If the global cost function stops decreasing, the second step associated with operator E w,2 starts by following equations (4.14) -(4.17). The second step determines the paths p w,min , p w,max ∈ P w in an alternative way: p w,min = random arg min p∈Pw (c ′ p , N p w,min ) (4.14) p w,max = random arg max p∈Pw,Fp>0 (c ′ p , N pw,max ) (4.15)

Figure 4 . 4 :

 44 Figure 4.4: The two-step optimization algorithm for 2D SATA problem.

3 : 4 :Figure 4 . 5 :

 3445 Figure 4.5: Illustration of random optimization strategy. A neighbor solution is generated, then its objective function value is compared with the previous one to determine if it can be accepted.

Figure 4 . 6 :

 46 Figure 4.6: Representation of the two-layer network, in which the first layer contains nodes, 28 links, and 2 OD pairs, and the second layer contains 24 nodes, 53 links, and OD pairs.

Figure 4 . 7 .

 47 The initial global cost function value of these two networks are 367,225 and 245,940, respectively. It shows that the two curves decrease significantly within 10 iterations in the first step of optimization, indicating the high convergence rate of DA. The second step of optimization begins when the global cost function stops decreasing, marked by a red arrow in Figure 4.7. The overall trend of the global cost function continues to decrease until the stop criterion is met. The small fluctuations are linked to the randomness of the RO algorithm. For both cases, the algorithm stops because the number of iterations without progress reaches N p . The minimum value of the global cost during the optimization process of the two networks reaches 9,068.5 and 12,963

  (a) UAM route network G 1 . (b) UAM route network G 2 .

Figure 4 . 7 :

 47 Figure 4.7: Evolution of global cost optimized by the proposed two-step optimization approach for two segregated networks.

Figure 4 .

 4 8b. The final result shows that few paths are assigned with a large volume of flows, and the flows on other paths are zero or close to zero. In the first network, only four paths are assigned with flows. These four paths are respectively p 6 w 1,1 = e 1,1 -e 1,4 -e 1,9 -e 1,6 -e 1,21 , p 9 w 1,1 = e 1,1 -e 1,4 -e 1,10 -e 1,12 -e 1,21 , p 2 w 1,2 = e 1,7 -e 1,9 -e 1,6 -e 1,19 and p 6 w 1,2 = e 1,7 -e 1,11 -e 1,17 -e 1,25 -e 1,26 -e 1,19 . Their average length is 189.3 km, and the average flight time is 0.96 h. By contrast, the average length of all 66 paths is 216.2 km, and the average flight time is 1.08 h. In the second network, five paths are distributed with significant flow, they are p 4 w 2,1 = e 2,40 -e 2,41 -e 2,43 -e 2,47 -e 2,53 , p 23 w 2,2 = e 2,1 -e 2,4 -e 2,7 -e 2,15 -e 2,32 -e 2,44 , p 38 w 2,2 = e 2,1 -e 2,5 -e 2,13 -e 2,17 -e 2,32 -e 2,44 , p 52 w 2,2 = e 2,1 -e 2,5 -e 2,15 -e 2,31 -e 2,41 -e 2,44 , p 53 w 2,2 = e 2,1 -e 2,5 -e 2,15 -e 2,32 -e 2,44 . Their average length is 406 km, and the average flight time is 1.88 h. While the average length of all 61 paths is 454 km, and the average flight time is 2.27h. As for both G 1 and G 2 , the result indicates that these paths are efficient in terms of distance and flight time.

  (a) UAM route network G 1 . (b) UAM route network G 2 .

Figure 4 . 8 :

 48 Figure 4.8: Optimized allocation of path flows for two segregated networks.

Figure 4 .

 4 12 is further provided for clearer visualization, in which the links and nodes are labeled as the symbols used in the context, and only the links with the nonzero flow or nonzero link cost are kept. For example, e 1,18 in Figure4.11a is not allocated with the flow but has a nonzero link cost. The link cost is measured by different colors extracted from the colorbar, and the width of the link represents the link flow. For comparison, the link costs under initial flow allocation are depicted in Figure4.10. Note that, the colorbar scales in these figures are different. We can see clearly that the cost for links initially assigned with the flow is extremely high, especially in dense areas. In Figure (a) UAM route network G 1 . (b) UAM route network G 2 .

Figure 4 . 9 :

 49 Figure 4.9: Optimized allocation of link flows for two segregated networks.

  (a) UAM route network G 1 (b) UAM route network G 2

Figure 4 . 10 :

 410 Figure 4.10: Representation of link costs and link flows under initial flow allocation in two segregated networks.

  (a) UAM route network G 1 (b) UAM route network G 2

Figure 4 . 11 :

 411 Figure 4.11: Traffic assignment result in terms of flow allocation, and link cost for two segregated networks.

  (a) UAM route network G 1 (b) UAM route network G 2

Figure 4 . 12 :

 412 Figure 4.12: Traffic assignment result for links with the nonzero flow or nonzero link cost for two segregated networks.

  (a) 3D representation of two-way UAM volume segments (b) The top view of Figure 1a

Figure 5 . 1 :

 51 Figure 5.1: An example of two-way UAM volume segments

1 )

 1 Measure the node complexity by traffics within the surrounding airspace 90 SATA in 3D UAM route network with connected layers To capture the traffic situation around each node v ∈ V, we define cylindrical airspace 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 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 backward directions relative to 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 ρ c v which are inwards towards the node v, and divergent density points ρ d v which are outward from the node v, where ρ c v ∩ ρ e ̸ = ∅ and ρ d v ∩ ρ e ̸ = ∅ if node v is connected by link e. It should also be noted that ρ c v ̸ ⊂ ρ e and ρ d v ̸ ⊂ ρ e since a cylinder A v may also include density points ρ e ′ located on other links e ′ ̸ = e. Figure 2 illustrates an example of 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.

Figure 5 . 2 :

 52 Figure 5.2: A example of cylindrical airspace, density points, and eigenvalue calculation of the associated LDS.

  Figure A.1 in Appendix A.5 shows the evolution of the 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 A.5.

Figure 5 . 3 :

 53 Figure 5.3: The optimization framework of the two-phase proposed algorithm on a simulation environment.

Figure 5 . 4 :

 54 Figure 5.4: Probabilistic performance-based neighborhood generation and evaluation of objective function for each transition in SA.

Figure 5 . 5 :

 55 Figure 5.5: A map-based visualization of the UAM route network in Singapore's urban airspace.

Figure

  Figure 6b 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 used to compute the node complexity. The size of the volume segment is determined in the presence of uncertainty (positioning accuracy, safety margin, emergency response, etc.). According to the regulation of unmanned aircraft (CAAS, 2019), the parameters of the simulation

Figure 5 . 6 :

 56 Figure 5.6: Graph representation of UAM route network used in this study.

Figure 5 . 7 :

 57 Figure 5.7: The distribution of demand for OD pairs in two cases.

  ), two indicators are introduced, including average travel time t:

  . The convergence speed of DA is quite fast: the objective function value decreases significantly in less than 100 iterations. This fact demonstrates that SA (a) Case 1.(b) Case 2.

Figure 5 . 8 :

 58 Figure 5.8: Complexity evolution of the proposed two-phase algorithm.

Figure 5 . 9 :

 59 Figure 5.9: Complexity evolution for two phases of the proposed algorithm in the first case.

  .

Figure 5 . 10 :

 510 Figure 5.10: Detailed complexity evolution for two phases of the proposed algorithm in the second case.

Figure 5 . 11 :

 511 Figure 5.11: Complexity evolution of DA without initialization phase in two cases.

  . The complexity of paths associated (a) Case 1. (b) Case 2.

Figure 5 . 12 :

 512 Figure 5.12: Complexity evolution of DA initialized by MSA in two cases.

  (a) Initial path flow pattern. (b) Path flow pattern generated by SA. (c) Final path flow pattern.

Figure 5 . 13 :

 513 Figure 5.13: Evolution of path flow and path complexity of the first case.

  (a) Initial path flow pattern. (b) Path flow pattern generated by SA. (c) Final path flow pattern.

Figure 5 . 14 :

 514 Figure 5.14: Evolution of path flow and path complexity of the second case.

  (a) Initial link flow. (b) Final link flow.

Figure 5 . 15 :

 515 Figure 5.15: Link flow evolution of the first case in the network representation.

Figure 5 . 16 :

 516 Figure 5.16: Link flow evolution of the second case in the network representation.

  (a) Initial node complexity. (b) Node complexity after assignment.

Figure 5 . 17 :

 517 Figure 5.17: Node complexity evolution of the first case in the network representation.

  (a) Initial node complexity. (b) Node complexity after assignment.

Figure 5 . 18 :

 518 Figure 5.18: Node complexity evolution of the second case in the network representation.

Figure 6 .

 6 1a gives an example for calculating the travel time from the origin of a path to a point g under network representation.Since τ p,g > 0, then ∃N ∈ N and τ ε ∈ [0, T d [ such that τ p,g = N T d + τ ε , (6.9)

  (a) Network representation of traversal time from the origin O to a point g. (b) The real-time path flow through point g given the path flow for two time periods F p,0 and F p,1 when the travel time τ p,g is less than the duration of time interval T d .

Figure 6 . 1 :

 61 Figure 6.1: Illustration of the travel time and the real-time path flow at point g on a path p connecting OD pair w in an example network.

Figure 6 . 2 :

 62 Figure 6.2: Example of cylindrical airspace around a node to measure the complexity of air traffic.

Figure 6 . 3 :

 63 Figure 6.3: Flowchart of a simulation-based framework for solving the DATA optimization problem. The simulation module includes the UAM representation of the route network, and the optimization module contains the optimization algorithm in a nested loop over OD pairs and rolling horizons.

. 4 ,

 4 only the current and last time intervals are involved in calculating C k in time intervals k = 1, • • • , 5. The path flow is updated in each OD according to the cost function or performance indicators.

Figure 6 . 4 :

 64 Figure 6.4: General strategy in the optimization module to approximate the optimal flow pattern of DATA using rolling horizon approach. This example includes six time intervals and the maximum time gap is less than the duration of the time interval. Inner loops are conducted for each OD pair in each time interval.

1 : 3 :

 13 procedure NeighborhoodGeneration(γ w,k ) γw,k ← γ w,k 4: p 1 ← path randomly selected from P with non-zero flow 5: p 2 ← path randomly selected from P w different from p 1 6:

  (a) Connection of OD pairs between delivery sites and vertiports. (b) Network representation.

Figure 6 . 5 :

 65 Figure 6.5: 3D UAM route network in Singapore's urban airspace.

Figure 6 . 6 :

 66 Figure 6.6: Simulated hourly UAM traffic demands for each OD pair from 09:00-10:00 to 13:00-14:00.

6. 4 . Experiments 123 a

 4123 transition in different temperatures throughout the optimization process of SA.

Figure 6 . 7 :

 67 Figure 6.7: Run time for a transition in PSA with the different number of threads. The median is indicated in green. The lower bound and upper bound of the box represent the first and third quartiles, respectively. The whiskers above and below the box include the values between the 5th and 95th percentiles.

Figure 6 . 8 :

 68 Figure 6.8: The reduction of objective function value after optimization of different algorithms compared to the initial state and the performance comparison between PM and DA. The boxplot is represented in the same way as Figure 6.7.

Figure 6 . 9 :

 69 Figure 6.9: Initial path flows distribution in different time intervals. The x-axis represents the index of paths and the y-axis represents the index of time intervals from 09:00-10:00 to 14:00-15:00. The dashed white lines separate the OD pair. The color of each path at each time interval indicates the associated path flow.

Figure 6 . 10 :

 610 Figure 6.10: Distribution of path flows optimized by PSA in different time intervals. The configuration is the same as Figure 6.9.

Figure 6 . 11 :

 611 Figure 6.11: The initial link flow distribution in the UAM route network in different time intervals. The color represents the link flow in a logarithmic scale that covers a large range of values.

Figure 6 . 12 :

 612 Figure 6.12: The resulting link flow distribution in the UAM route network in the different time intervals, optimized by PSA. The configuration is the same as Figure 6.11.

Figure 6 . 13 :

 613 Figure 6.13: The initial air traffic complexity and congestion are represented respectively on nodes and links of the UAM route network in different time intervals.

Figure 6 . 15 :

 615 Figure 6.15: Change in complexity cost in terms of each node compared with initial flow pattern and the flow pattern optimized by PSA.

Figure 6 . 16 :

 616 Figure 6.16: Change in congestion cost in terms of each link compared with initial flow pattern and the flow pattern optimized by PSA.

Figure A. 2 :

 2 Figure A.2: Eigenvalue loci for four typical 2D traffic situations. The location of the observed traffic is indicated by the black dots. The blue arrows indicate the speed vectors. The arrowhead indicates the direction and the length of the arrow indicates the speed.

Figure A. 3 :

 3 Figure A.3: Eigenvalue loci for four typical 3D traffic situations. The black dots represent the observed trajectory points of aircraft, the red arrows indicate their speed vectors, and the blue arrows indicate the speed vectors estimated by the LDS at equidistantlypartitioned points. The arrowhead indicates the direction and the length of the arrow stands for speed.

Table 1 .

 1 1: Estimated UAM demands for low and high scenario[START_REF] Anand | A Scenario-Based Evaluation of Global Urban Air Mobility Demand[END_REF].

	(a) Low scenario
	Years	2035 2040 2045 2050
	Annual UAM passenger trips (billion)	58.34 82.61 99.41 127.8
	Annual UAM passenger kilometers (billion) 716.4 1013 1211 1554
	Annual UAM utilization (billion hours)	6.043 8.447 10.09 12
	Annual UAM Vehicle Trips (Billion)	23.33 33.04 39.76 51.09
	(b) High scenario
	Years	2035 2040 2045 2050
	Annual UAM passenger trips (billion)	227.3 279.1 334.5 405.5
	Annual UAM passenger kilometers (billion) 2762 3408 4098 4990
	Annual UAM utilization (billion hours)	11.51 14.21 17.07 20.79
	Annual UAM Vehicle Trips (Billion)	90.94 111.6 133.8 162.2

  2. Scheduling and routing become more flexible. Since UAVs have different specifications, they can perform various tasks as required. UAVs have a more comprehensive range of movement and require much less navigation than manned aircraft;3. UAVs may contribute to congestion reduction. For example, in urban transportation, large-scale applications of UAVs for urban delivery missions could potentially help reduce traffic congestion;

4. UAV services are cost-saving. No longer exclusive to the military, UAVs are now a common tool for civilian use. It is reasonably priced and very economical in terms of maintenance and fuel costs.

  ; Farges and Delahaye, 2001; Deschinkel, Farges, and Delahaye, 2002; Delahaye, Sofiane, and Puechmorel, 2005; Nosedal et al., 2014), Air Traffic Flow Management (ATFM) (Strub and Bayen, 2006; Bertsimas, Lulli, and Odoni, 2011; Zhang, Cai, and Zhu, 2012; Zhang et al., 2015), en-route network management (Delahaye and Puechmorel, 2013a; Haouari, Aissaoui, and Mansour, 2009), etc. Some other studies have been carried out for different purposes: to increase the airspace capacity, reduce the noise level (Netjasov, 2008; Ganić et al., 2018; Chatelain and Van Vyve, 2018; Ho-Huu et al., 2019), and reduce emission (Economou et al., 2007; Mirosavljević, Gvozdenović, and Čokorilo, 2011; Ho-Huu et al., 2019). Air traffic assignment in the application of the UTM is mainly focused on UAV task assignment (Jiang, Zhou, and Ye, 2017; Zhou et al., 2018; Cheng et al., 2019; Liu et al., 2019).

  Number of entry points;

	-Route configurations;
	-Route intersections;
	-Arriving traffic distribution and mix;
	-Flight levels.
	• Dynamic density: (Laudeman et al., 1998; Sridhar, Sheth, and Grabbe, 1998):
	Dynamic density metric is based on the flow characteristic in the airspace. As
	correlated factors, the on-duty controller activity is captured. Generally, a lin-
	ear regression model is designed to approximate the experienced workload. The
	following variables were selected as inputs of the dynamic density function for a
	sector:
	-Traffic density;
	-Number of aircraft with different filter conditions such as changes in direction,
	speed, altitude, etc.

  further enhanced this model under uncertain demands. Since then, a variety of traffic assignment models have been developed for solving routes and slot allocation problems (Delahaye and Odoni, 1997; Farges and Delahaye, 2001; Deschinkel, Farges, and Delahaye, 2002; Delahaye, Sofiane, and Puechmorel, 2005; Nosedal et al., 2014), ATFM (Strub and Bayen, 2006; Bertsimas, Lulli, and Odoni, 2011; Zhang, Cai, and Zhu, 2012; Zhang et al., 2015), en-route network management (Delahaye and Puechmorel, 2013a; Haouari, Aissaoui, and Mansour, 2009), etc. For example,

  the context of UTM, most research has focused on UAV task assignment (Jiang, Zhou, and Ye, 2017; Zhou et al., 2018; Cheng et al., 2019; Liu et al., 2019). In addition, a few works are based on the UAM route network. Chatelain and Van Vyve (2018) investigates the problem of user-demand-based UAV assignment over a highly dense geographic area. Economou et al. (

Table 3 .

 3 1: Nomenclature specified to Chapter 3 , δ e,u , δ e,d = 1 if the link e is horizontal, upward, and downward, resp.; o/w 0

	Term	Description
	A a , A d , A p	Aerodromes, danger areas, and prohibited areas
	d e	Length of link e
	ϕ h , ϕ u , ϕ d	Efficiency coefficients for horizontal, upward, and downward links
	L r	Sound Pressure Level (SPL) at a distance of r from the source
	L e,S	Noise impact from a link e to a nearby subzone S
	C h,e	Cost for crossing hazardous airspace for link e
	C f,e	Cost of flight efficiency for link e
	C n,e	Total noise impact from link e to all nearby subzones
	C p	Cost of path p
	α a , α d , α p	Penalization parameter for crossing aerodromes, danger areas, and prohibited areas, resp.
	δ e,n	= 1 if link e traverses hazardous airspace n, o/w 0
	δ e,h

48 UAM route network design in low-altitude airspace

Table 3 .

 3 2: Basic topology features of different altitude layers and the whole network

	Network	# of nodes # of links AD	ACC	DOA
	Network layer at 150ft 75	216	5.7600 0.3138 -0.0004
	Network layer at 175ft 102	336	6.5882 0.3875 0.3613
	Network layer at 200ft 111	444	8.0000 0.4709 0.2139
	Network layer at 225ft 115	524	9.1130 0.5081 0.2775
	Network layer at 250ft 115	546	9.4957 0.5154 0.3383
	The whole network	529	2218	8.3856 0.3789 0.3789

Table 4 .

 4

		1: Nomenclature specified in Chapter 4
	Notation Description
	A	a certain layer of urban airspace
	A e,n	n-th complexity area of link e
	c ′ p	marginal cost of path p
	e i,j	link connecting nodes i and j
	E w	equilibrium operator on OD pair w
	f e	flow on link e
	F i w	flow on i-th path of OD pair w
	Q e	evaluation points of link e
	s g	speed on grid point g
	U e	UAM corridor of link e
	x v	x coordinate of node v
	y v	y coordinate of node v
	λ	dominant eigenvalue of a linear dynamical system
	ρ e	density points on link e
	ω g	density on grid point g

UAM route network with segregated layers

13 :

 13 P ← Matrix with 2 rows containing all combinations of append(x s , x e ) and y

	14:	P e g ←	cos(θ

e ) -sin(θ e ) sin(θ e ) cos(θ e ) P + x i y i ▷ column-wise addition 15:

  Within a high-density UAM operational scenario, the demands of each OD pair (aircraft/h) for G 1 and G 2 are set as:• Set of demands D 1 = {d w 1,1 ,d w 1,2 } with d w 1,1 = 1925 and d w 1,2 = 2255; D 2 = {d w 2,1 , d w 2,2 } with d w 2,1 = 1995 and d w 2,2 = 2160

1 w 1,1 , ..., p 55 w 1,1 }, P w 1,2 = {p 1 w 1,2 , ..., p 11 w 1,2 }; P 2 = P w 2,1 ∪P w 2,2 with P w 2,1 = {p 1 w 2,1 , ..., p 7 w 2,1 }, P w 2,2 = {p 1 w 2,2 , ..., p 54 w 2,2 }

Table 4 . 2 :

 42 Parameters setting for route network modeling and optimization process in the experiment.

	Process	Param. Value	Description
		L l	2 (km)	Length of grid
		L w	0.02 (km)	Width of grid
	Route network	s e	200 (km/h) Speed for link e
	modeling	N s	2	Number of lateral points on each side
		β	10 (km -1 )	Unit conversion coefficient
		R e	5 (km)	Radius of complexity area for link e
		h	0.01	Small variation for numerical derivation
		N t	20	Sampling epoch in RO strategy
	Optimization algorithm	N r D N m	2 U (a, b) 1000	Number of alternative suboptimal paths Sampling distribution in RO strategy Total number of iteration
		N p	200	Number of iteration without progress
		v s	100	Target value of criteria
	Note: e ∈ E1 ∪ E2

  pairs of these two networks w ∈ W 1 and w ∈ W 2 is illustrated in Figure 4.8. The paths {p 1 , ..., p 11 w 1,2 } are indexed from 1 to 66 as x-axis in Figure 4.8a, and {p 1 w 2,1 , ..., p 7 w 2,1 , p 1 w 2,2 , ..., p 54 w 2,2 } are indexed from 1 to 61 as x-axis in

	4.5. Results	79
	p 1 w 1,2
		w 1,1 , ..., p 55 w 1,1 ,

Table 5 .

 5 1: Notations and definitions specified in Chapter ??

	Notations Explanations
	Sets
	A v	Cylindrical airspace around node v
	ρ c v , ρ d v	Convergent, divergent density points in A v of node v
	ρ e	Set of density points of link e
	γ	Set of path flow proportions
	f	Set of link flows
	Variables
	F i w	Flow on i-th path of OD pair w
	F p	Flow on path p
	f e	Flow on link e
	ω g	Density on density point g
	γ i w

Table 5 .

 5

				2: Parameters setting
	Module	Parameter Value		Description
		L h	0.3 (km)	Length of grid on horizontal volume seg-
					ment
		L v	2 × 10 -3 (km) Length of grid on vertical volume seg-
					ment
		L w	0.03 (km)	Width of grid
		s h	120 (km/h)	Speed for horizontal link e h
		R max	20 (km)	Maximum range of UAV
	Simulation module	T max s v N s	2000 (s) 5 (km/h) 3	Maximum endurance of UAV Speed for vertical link e v Number of lanes for each link
		β	3		Efficiency threshold parameter
		R e	5 (km)	Radius of complexity area for link e
		F t	10		Meaningful flow threshold e
		φ h	2 × 10 -4	Efficiency coefficient for horizontal links
		φ u	10 -3		Efficiency coefficient for upward links
		φ d	2 × 10 -5	Efficiency coefficient for downward links
		T 0	1000		Initial temperature
		N tr	100		Number of iterations at each tempera-
	Initialization				ture
	module	α	0.9		Cooling rate
		Stopping criterion	Cp-Cc Cp	≤ 1%	The best objective value of current tem-perature is improved less than 1%
	Optimization	ε	0.01		Small value of variation for numerical
	module				derivation
	satisfied, and the random number generator is kept the same for these two cases:
				F p =	d w rand p rand p
				p∈Pw

Table 5 .

 5 3: Performance comparison of models in terms of 95% CI for computation time, complexity, and flight efficiency in case 1.

	Models	Case 1 Computation	Total complexity	Percentage reduction	Average travel	Average path	Energy
		time (s)		in complexity	time (s)	length (km)	consumption
	Initial state -	680081.00 ± 2000.58	-	404.59 ± 0.32	12.61 ± 0.01	228.06 ± 0.18
	AON	< 10 -4	509223.03	25.12%	248.00	7.16	130.27
	SA	5.86 ± 0.26	192021.02 ± 1613.07	71.76% ± 0.25%	307.27 ± 0.38	10.11 ± 0.01	182.72 ± 0.27
	DA	38.58 ± 0.79	98957.43 ± 2121.67	85.45% ± 0.31%	320.10 ± 1.83	10.07 ± 0.06	181.11 ± 1.15
	MSA	0.92 ± 0.02	228767.25 ± 524.08	66.36% ± 0.12%	309.19 ± 0.07	10.17 ± 0.00	183.26 ± 0.04
	MSA+DA	28.55 ± 1.44	159262.94 ± 8523.64	76.58% ± 1.26%	316.88 ± 1.48	10.32 ± 0.05	185.68 ± 0.84
	Our model	13.62 ± 0.52	65005.20 ± 3552.16 90.44% ± 0.53%	323.83 ± 2.52	10.27 ± 0.08	184.68 ± 1.43

Table 5 .

 5 4: Performance comparison of models in terms of 95% CI for computation time, complexity, and flight efficiency in case 2. convergence speed can be significantly improved. The DA only takes 7.76 ± 0.43 and 7.83 ± 0.39 seconds, 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.

	Models	Case 2 Computation	Total complexity	Percentage reduction	Average travel	Average path	Energy
		time (s)		in complexity	time (s)	length (km)	consumption
	Initial state -	665702.86 ± 2017.43	-	405.60 ± 0.30	12.66 ± 0.01	220.04 ± 0.17
	AON	< 10 -4	510435.59	23.32%	258.21	7.48	131.06
	SA	6.29 ± 0.34	183251.32 ± 1694.38	72.47% ± 0.27%	314.62 ± 0.29	10.34 ± 0.01	179.66 ± 0.27
	DA	35.31 ± 0.89	100567.56 ± 2869.64	84.89% ± 0.80%	323.41 ± 1.55	10.13 ± 0.05	175.01 ± 0.55
	MSA	0.98 ± 0.03	216930.25 ± 686.09	67.41% ± 0.15%	317.10 ± 0.08	10.45 ± 0.00	180.89 ± 0.06
	MSA+DA	27.53 ± 0.87	154941.17 ± 3293.94	76.73% ± 0.50%	319.13 ± 0.82	10.48 ± 0.03	183.58 ± 0.45
	Our model	14.12 ± 0.33	52450.41 ± 2376.11 92.12% ± 0.35%	328.58 ± 1.84	10.39 ± 0.06	179.49 ± 1.11
	framework, the

Table 6 .

 6 1: Nomenclature specified to Chapter 6. Free-flow travel time from the origin of a feasible path p to a point g ∈ p F g,p,kAverage flow that traverses a grid point g on path p of OD pair w in time interval k

	Term	Description
	Air route network
	U e	UAM volume segment of link e
	A v	Cylindrical airspace around node v
	g ∈ g e	A point in the set of grid points of link e
	Z v	Intersection zone of node v
	γ	Set of proportions for all paths flows in all OD pairs during the planning horizon
	γ w,k	Set of proportions for all path flows in OD pair w in time interval k
	γ p,k	Proportion for path flow in OD of p in time interval k
	τ p,g
	Fg,k	Accumulative path flow on a grid point g in time interval k
	f in e,k , f out e,k f c e	Inflow and outflow of link e in time interval k Traffic flow capacity for link e
	D w,k	Demand of OD pair w in time interval k
	ρ g,k	Traffic density of point g in time interval k
	s e,k	Average speed of UAM traffic flow on link e at time interval k
	α e	Energy consumption factor of link e
		109

air traffic assignment model for high-density UAM operations

Table 6 . 2 :

 62 Basic topology features of different altitude layers and the UAM route network.

	Network	Number of nodes Number of links AD ACC DOA
	Network layer at 150ft	51	96	4.57 0.35 -0.19
	Network layer at 175ft	62	145	5.27 0.41 -0.12
	Network layer at 200ft	85	216	5.20 0.32 0.31
	Network layer at 225ft	89	199	4.47 0.24 0.23
	Network layer at 250ft	45	89	3.96 0.18 0.33
	The whole route network 332	830	5.00 0.23 0.23

Table 6 . 3 :

 63 Parameters setting in the experiment.

	Module	Parameter	Value	Description
		L h	0.15 (km)	Grid length of horizontal volume segment
		L v	1e-3 (km)	Length of grid on vertical volume segments
		L w	0.03 (km)	Grid width
		s h	54 (km/h)	Speed for horizontal link e h
		s v	15 (km/h)	Speed for vertical link e v
		N s	3	Number of lateral points on each side
		N r	3	Number of rows of the grid included by cylin-
	Simulation			drical airspace
	module	F t	20 (veh/h)	Meaningful flow threshold e
		f c e	1.5e4 (veh/h) Traffic flow capacity for link e
		φ h	0.2	Efficiency coefficient for horizontal links
		φ u	1	Efficiency coefficient for upward links
		φ d	0.02	Efficiency coefficient for downward links
		ϕ	0.9	Weighting parameter in objective function
		T init	100	Initial temperature
	Optimization module	N tr α ε	100 0.9 0.01	Number of iterations at each temperature Cooling rate Small variation for numerical derivation in
				the DA

Dynamic air traffic assignment model for high-density UAM operations distribution

  of hourly demand between several ODs for an urban network during a day (Deng, Zeng, and Mei, 2019).

Table 6 .

 6 4: Performance comparison of PSA and representative conventional DTA algorithms in terms of computation time, objective function, and flight efficiency for a 95% CI.

	Models Computation	Value of objective	Complexity cost	Congestion cost	Average travel	Average path
		time (s)	function			time (s)	length (km)
	Initial -	99081.42 ± 368.48	585380.95 ± 3627.11	45048.14 ± 27.92	1201.94 ± 0.41	17.84 ± 0.01
	PSA	143.47 ± 3.74 67157.45 ± 146.87	302743.55 ± 1574.53	40981.22 ± 36.90	1233.61 ± 0.48	18.28 ± 0.01
	PM	12.315 ± 4.30 95926.01 ± 122.50	557388.35 ± 1260.06	44652.42 ± 65.67	1231.31 ± 0.57	18.23 ± 0.01
	DA	502.15 ± 6.68 94361.08 ± 236.84	541264.75 ± 2368.71	44705.12 ± 26.23	1198.18 ± 0.38	17.78 ± 0.01

A methodology is proposed to design grid-based UAM route networks in lowaltitude urban airspace in the presence of obstacles and hazardous areas using a variety of open-source data;

2. In the designed UAM route network, the feasible routes with low similarity are

* The number of temperature transitions reaches a specific limit; * The temperature reaches a specific limit; * The CPU time reaches a specific limit; * The improvement in terms of the best objective value with respect to the previous temperature is less than a percentage threshold; * No acceptance or improvement in the neighborhood after a specific number of transitions; * No acceptance or improvement in the neighborhood after a specific CPU time; * The combination of the aforementioned criteria.In summary, transportation network theory proposes different modeling technologies for most real situations. However, they have several limitations. The objective function in some optimization problem formulation can be non-convex, nonlinear, and even cannot be represented explicitly as a function as decision variables. The aforementioned resolution algorithms designed for traffic assignment are not capable of solving these

Acknowledgements

This thesis is partially supported by the research project CONCORDE of the Defense Innovation Agency (AID) of the French Ministry of Defense (2019650090004707501). Many thanks to my friends for always supporting me! My Deepest thanks to my girlfriend and soulmate Shan, for

Abbreviations

ACC

Average Clustering Coefficient

AD

Average Degree with fewer uncertainties in the search space.

• Other than most traffic assignment algorithms based on several simulation environments, only one simulation environment will be created for each individual process, which improves the calculation efficiency.

Parallel Simulated Annealing (PSA) algorithm

As the DATA is an optimization problem with continuous variables and a nondifferentiable objective, from the operational point of view, a near-optimal solution is required, rather than a globally optimal solution, which exists theoretically but is difficult to achieve in the limited computational time. In this study, we redesign the classical SA algorithm and introduce a PSA algorithm with a novel neighborhood generation strategy and parallelized cost function evaluations to efficiently optimize the simulation-based DATA problem.

As indicated in section 2.3.2.2, the asymptotic convergence of the SA to the global optimum is highly dependent on the selection of neighborhood. For the DATA problem, we propose a novel strategy to generate the neighborhood solution. The pseudocode is presented in Algorithm 6. To generate the neighboring solution that moves effectively and efficiently, we focus on transferring the flow from one path to another in each OD pair in each transition. The outflow and inflow paths are chosen randomly, whereas the outflow path is chosen only from the paths with the non-zero flow. To avoid redundant operations on paths with less traffic, we determine the transferring flow σ by comparing the flow of outflow path γ p i ,k with the threshold F t /D w,k (line 6). If γ p i ,k is larger, the flow to be transferred will be drawn from a uniform distribution between zero and γ p i ,k . Otherwise, the outflow path will transfer all its flow. In this way, only constraints (6.14) and (6.15) are required to be ensured, and all other constraints are satisfied.

Given a neighboring solution, we need to update the objective function to determine whether to accept this solution. Due to the separability of the objective function (6.24) of rolling horizons, it is noteworthy that the objective function may be updated through the simulation module without involving all decision variables, by calculating the incremental difference of the objective function value. To further speed up the optimization algorithm, the objective function is evaluated in parallel. Unlike problem-independent parallelization methods that involve frequent evaluations of the objective function (Onbaşoğlu and Özdamar, 2001), we design a problem-dependent parallelization approach by dividing ∆C w,k among multiple threads. As the evaluation of the objective function is computationally expensive, especially in solving the WMMSE and computing the eigenvalues, the objective function is only calculated if necessary and updated incrementally. The pseudocode of this parallel computing strategy is formulated in Algorithm 7. After completely computing the objective function (6.21), for each OD pair w and rolling horizon k, we first generate the neighboring solution and submit the evaluation of associated complexity costs to the task queue (lines 2-4). The task queue is executed concurrently Appendix A

Linear dynamical systems A.1 Definition and notation of linear dynamical system

Within a continuous-time model scope, the LDS is given by the following form:

Given N observations with positions x i , speeds v i , and densities ω i , i = 1, ..., N , the Weighted Minimum Mean square Error (WMMSE) between the LDS and observations is formulated as:

Next, two cases including the 2D case and the 3D case are introduced, and the notations are given separately. For both cases, the coefficient matrix and vector of the LDS are grouped in a single matrix C:

In two dimensions, x i = (x i , y i ) ⊤ and v i = s i (cos(θ i ), sin(θ i )) ⊤ , i = 1, ..., N . In order to transform Equation (A.2) into matrix form, the following matrices are defined:

In three dimensions,

Linear dynamical systems to transform equation (A.2) into matrix form, the following matrices are defined:

A.2 Weighted minimum mean square error estimation of linear dynamical system

Denote the dimension of a linear dynamical system as D LDS . The error criterion E can be reformulated as

where ∥∥ F represents the Frobenius norm.

The WMMSE aims to find the matrix Ĉ that minimizes the error criterion:

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

∇ C E = 0 allow us to calculate Ĉ. If XX ⊤ is invertible, namely the columns of X are linearly independent, the LDS problem has a unique solution:

where X + is the pseudo-inverse of X:

The matrix X may not be full-rank or the LDS could be ill-conditioned in some rare A.2. Weighted minimum mean square error estimation of linear dynamical system 155 cases. To avoid numerical problems associated with the determinant when inversing the matrix XX ⊤ , Singular Value Decomposition (SVD) was proposed to solve this problem. X can be decomposed as:

where L ∈ R D LDS +1×D LDS +1 and R ∈ R N ×N are unitary matrices. L and R are not unique, which can be composed respectively by the eigenvectors of XX ⊤ and X ⊤ X. Σ ∈ R D LDS +1×N is a rectangular diagonal matrix with non-negative values on the diagonal.

LDS +1 is a symmetric matrix, and the square roots of its eigenvalues are the singular values of X. The singular values of X can be sorted in descending order:

where r is the rank of X and r ≤ min(D LDS + 1, N ). Then, Σ is then defined as:

where the elements of the last N -r columns and D LDS + 1 -r rows are 0.

According to Equation (A.17), X + can be formulated as :

can be formulated as:

Linear dynamical systems

where the elements of the last D LDS + 1 -r columns and N -r rows are 0. Ĉ is thus given by:

Finally, Â can be extracted from Ĉ according to Equation (A.3).

A.3 Partial qualitative behavior of LDS characterized by an eigenvalue Proposition A.1 A particular eigenvalue λ of the coefficient matrix A reflects the associated component of the state of a continuous-time LDS:

• If Re(λ) > 0, then lim t→+∞ x(t) = ∞, which means the component is unstable and diverging,

• If Re(λ) < 0, then lim t→+∞ x(t) = 0, which means the component is stable and converging,

• If Re(λ) = 0, the component is conserved. Depending on the imaginary part, it can be further classified into:

Proof. In order to generalize all forms of LDS by a constant-free representation, let A closed-form solution of equation (A.22) for any square matrix B is:

Denote the dimension of B as N . B is assumed to be digonalizable with independent eigenvectors v 1 , . . . , v N and the associated eigenvalues λ 1 , . . . , λ N , then the initial state of the system can be represented as:

By applying the following property,

the state can be formulated as:

Thus, the corresponding eigenvector component of a particular eigenvalue

Re(λ i ), this property implies Proposition A.1.

A.4 Representative qualitative behavior of LDS characterized by the dominant eigenvalue

Proposition A.2

The dominant eigenvalue of A determines the stability of the LDS.

Linear dynamical systems

Proof. Following the constant-free form of LDS in equation (A.22), let λ be the dominant eigenvalue of B, then equation (A.28) can be written as:

Because the real part of λ i -λ is negative, we have:

According to equation (A.31), A dominant eigenvalue λ determines the stability of the whole LDS.

A.5 Eigenvalue loci for typical three-dimensional traffic situations

In Figure A.2 and A.3, four typical scenarios of 2D and 3D traffic situations are presented, respectively. These cases include translation, rotation, convergence, and divergence. All traffic situations can be decomposed into one or several typical scenarios. The eigenvalues associated with 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. Concerning the rotation case, the aircraft moves around a circle but the relative distances between aircraft remain constant. All associated eigenvalues have a zero real part. The eigenvalues associated with these two cases are located within the strip in Figure A.1, 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 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 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 on 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 the matrix the coefficient matrix A of the associated LDS have a small real part [START_REF] Delahaye | Air traffic complexity map based on non linear dynamical systems[END_REF]. This corresponds to an organized traffic situation where the traffic has no or few interactions.

A.5. Eigenvalue loci for typical three-dimensional traffic situations