Joint Optimization of Cost and Scheduling for Urban Air Mobility Operation Based on Safety Concerns and Time-Varying Demand

Abstract

As the value and importance of urban air mobility (UAM) are being recognized, there is growing attention towards UAM. To ensure that urban air traffic can serve passengers to the greatest extent while ensuring safety and generating revenue, there is an urgent need for a transportation scheduling plan based on safety considerations. The region of Beijing–Tianjin–Hebei was selected as the case study in this research. A real-time demand transportation scheduling model for a single day was constructed, with the total service population and total cost as objective functions, and safety intervals, eVTOL performance, and passenger maximum waiting time as constraints. A Joint Optimization of Cost and Scheduling Particle Swarm Optimization (JOCS-PSO) algorithm was utilized to obtain the optimal solution. The optimal solution obtained in this study can serve 138,610,575 passengers during eVTOLs’ entire lifecycle (15 years) with a total cost of CNY 368.57 hundred million, with the cost of CNY 265.9 per passenger. Although it is higher than the driving cost, it saves 1–1.5 h and thus has high cost effectiveness during rush hours.

1. Introduction

Urban air mobility (UAM) has significant potential for reducing ground traffic congestion and generating revenue. Studies indicate that during rush hours, UAM can alleviate at least 45% of pedestrian flow [1]. Additionally, NASA projects UAM to be profitable by 2028, with an estimated 740 million passenger trips via aerial buses by 2030 [2]. Moreover, it is anticipated to handle 500 million last mile delivery transport services by 2030, potentially generating billions of dollars [2]. In this context, UAM is advocated vigorously in America [3], China [4], and Europe [5].

For the maximum efficacy of UAM, it is crucial to implement a well-organized dispatch system based on UAM operational models for electric Vertical Takeoff and Landing (eVTOL). The operational model of UAM entails real-time responsiveness and sensitivity to safety concerns. Therefore, establishing a scheduling plan centered around safety and demand is fundamental for maximizing the efficacy of UAM.

Several studies have been conducted to solve this problem. Z. Wu et al. [6] proposed a charging sequencing and scheduling model with a given quantity of eVTOLs considering the battery state of charge and the maximum flight range constraints to minimize the waiting time of passengers. Through changing the number of charging facilities, they also have analyzed the relationship between the number of charging facilities and the waiting time of passengers. Although careful considerations have been made in the aspect of charging limitations in their research, Z. Wu et al. [6] ignore the safety constraints during the operation of eVTOLs. Furthermore, due to the fact that in the real world, customers will not always wait at the vertiport until an eVTOL comes, it is unrealistic to simply use the waiting time of passengers as an objective function. Charging measures, which are split into those that are delayed and canceled, have also been considered by Z. Guo et al. [7] in their joint routing and charging strategy model. The model considering flexible charging measures was further solved through a branch and price algorithm taking the recovery cost and cancelling cost as objective functions. However, rather than maximizing services for passengers, improving service quality, and safety insurance, this study [6] focuses more on how to effectively reallocate a given number of eVTOL aircraft, flights, and charging tasks during UAM disruptions to minimize costs and restore normal operations. Moreover, the performance of eVTOLs, passenger transport demand, and electricity prices were considered to simulate the real-world environment in the joint routing and charging model developed by S. Paul et al. [8] utilizing a Graph Neural Network. Compared with the studies of Z. Wu et al. [6] and Z. Guo et al. [7], the focus and main contribution of S. Paul et al.’s [8] research is on having proposed a model that could simulate real-world conditions better and solved it with an algorithm with considerable computation efficiency. However, no specific optimization plan for fleet size was proposed in all the studies mentioned above.

Based on the previous studies, a multiple objective planning model considering detailed real-world factors, including population density and income distribution, has been developed by S. Roy et al. [9] to obtain an optimal fleet size and schedule that maximize the total served demand and the revenue. The main focus of S. Roy et al. [9] is on how to determine the required fleet size by taking into account factors such as population density, income distribution, ticket prices, and operational costs, while also expanding revenue. However, the detailed consideration of safety constraints was not a concern in S. Roy [9] et al.’s research. An on-demand scheduling plan with a heterogeneous fleet has been developed by S. H. Kim et al. [10], utilizing particle swarm optimization and a genetic algorithm that integrates with a greedy algorithm. The performance of these two algorithms has also been compared. The article approaches the scheduling issue from a profitability perspective, introducing the concept of a heterogeneous fleet, which allows the fleet to consist of two or more types of eVTOLs while determining pricing strategies and maximizing revenue. However, constraints such as battery capacity, flight range, operational time, etc., have not been considered.

Moreover, this joint fleet size optimization with a heterogeneous fleet and scheduling problem has also been studied by M. Lindner et al. [11]. Apart from considering the performance of eVTOLs and energy consumption, M. Lindner et al. have also set a precedent in setting a slight tolerance for deviations from the predetermined arrival time. A hybrid simulation goal programming model was developed by S. Rajendran et al. [12] to solve the scheduling problem. In addition, to evaluate the performance of this schedule, the willingness to fly rate, the percentage of demand fulfillment, customer waiting time, etc., were taken into consideration in their study. Although the research conducted by M. Lindner et al. [11] and S. Rajendran et al. [12] has carefully considered optimization with a heterogeneous fleet and the scheduling problem in the aspect of eVTOL performance and operation, it still lacks a detailed consideration of safety.

An adaptive control system (ACS) model for the arrival sequencing and scheduling of eVTOL aircraft in a multi-vertiport system has been proposed by S. Quan et al. [13] to address the conflict between limited operational resources and traffic demand in urban air mobility. The limited battery energy supply of eVTOL aircraft, different operational processes, and safety challenges in terminal areas have been considered in this research. The energy and safety constraints have been considered in detail in their dynamic time-saving path planning model. Furthermore, to ensure the safety of eVTOLs during takeoff and landing, an optimization model for the arrival sequencing and scheduling of eVTOL aircraft in urban air mobility (UAM) has also been developed by I. Kleinbekman et al. [14]. The study focuses on calculating the optimal required time of arrival (RTA) for eVTOL aircraft to ensure safe separation and minimal delays under limited battery energy and vertiport (such as rooftop landing points) capacity conditions. Moreover, detailed factors such as the eVTOL type and performance, safety intervals, energy constraints, and the capacity of vertiports have been considered in the scheduling model developed by Pradeep, P et al. [15] Different to the studies of S. Quan et al. [13] and I. Kleinbekman et al. [14], a heterogeneous eVTOL fleet and vertiport capacity have been introduced to make it reflect the real-world conditions more precisely. However, although the studies of S. Quan et al. [13], I. Kleinbekman et al. [14], and Pradeep, P [15] have considered the safety issues in detail, their research is limited to the scheduling problem of eVTOLs in terminal areas rather than passenger service and operation issues.

The aforementioned studies approach the scheduling issue mostly from a profitability perspective and charging strategy. However, it is unrealistic to overlook safety concerns such as the time differences caused by the minimum safe spacing between aircraft to ensure the safe takeoff and landing of multiple eVTOLs [16]. Ignoring constraints in these aspects may result in collisions between aircraft and obstacles or other aircraft. Therefore, it is of great importance to have safety limitations concerned while scheduling eVTOLs. Even the studies that have considered safety constraints in detail just limit their research to the operation of the terminal area, rather than integrating safety into the entire flight scheduling to serve passengers to the greatest extent. Moreover, as most passengers will not wait until the arrival of eVTOLs no matter the waiting time, it is also unrealistic to choose the waiting time of customers as an objective function. Furthermore, the aforementioned articles do not provide a detailed calculation of the energy consumption of eVTOLs during different stages of the flight in the process of considering the eVTOL scheduling problem, but merely assume that energy consumption is linearly related to the flight distance. The contributions of this paper are delineated as follows:

Time-varying Demand:

• In the process of solving scheduling problems, consider the maximum tolerable waiting time of passengers as an input parameter rather than an objective function that needs to be optimized.
• Refine the real-time calculation of energy consumption during different flight stages between vertiports and integrate the energy consumption condition into the scheduling model.

Safety Concerns:

• Integrating the safety time interval constraint during takeoff and landing into the optimization of this timetable.
• Integrating the constraint that eVTOLs should maintain a remaining battery level of at least 30% of the total battery capacity into this joint optimization model.

The subsequent sections are outlined as follows. Section 2 provides the methodology and an illustration of the proposed scheduling model. The assumptions and hypotheses of this problem are displayed in Section 2.2. The algorithm utilized and details of this scheduling model are demonstrated in Section 2.3. Section 3 presents case studies based on data from Beijing, Tianjin, and Xiong’an (Hebei). Lastly, Section 4 encapsulates the conclusions.

2. Methodology

2.1. Symbols

The variables utilized in this article are demonstrated in

Table 1. The definition of variables.

Symbol	Definition (Unit)	Type
$Amn$	The number of vertiports (count)	Input Parameter
$AP{R}_i$	The aspect ratio of eVTOL i	Input Parameter
$AS{P}_{ij}$	The number of passengers that can be served at vertiport j	Intermediate Parameter
$ac{c}_i$	The acceleration of eVTOL i (m/s2)	Intermediate Parameter
$B$	Amount (count)	Input Parameter
$C$	The operation cost of eVTOLs (CNY)	Input Parameter
$c_{D_0}$	The parasite drag coefficient	Input Parameter
$c_L$	Lift coefficient	Input Parameter
$Dis$	Distance (km)	Input Parameter
$Dev$	The deviation of the CNS system	Input Parameter
$d_i$	The maximum range of eVTOL i (km)	Input Parameter
$E$	Energy (kW·h)	Input Parameter
$Fi{t}_i$	The fitness value of individual i	Output Parameter
$FOM$	Figure of merit	Input Parameter
$Force$	The drag force of eVTOL i (N)	Input Parameter
$f_i$	The function of passenger demand in vertiport i over time	Function
$g$	Gravity acceleration (m/s2)	Constant
$h$	Altitude (m)	Input Parameter
$L_s$	The safety interval distance between eVTOLs (m)	Input Parameter
$Ind$	The initial distance between eVTOL i and j (km)	Input Parameter
$lat$	Latitude	Input Parameter
$lon$	Longitude	Input Parameter
$MH{A}_i$	The maximum cruising altitude of eVTOL i (m)	Input Parameter
$Mac$	The maintenance cost per kilometer (CNY/km)	Input Parameter
$m_i$	The mass of eVTOL i (kg)	Input Parameter
$N_i$	The number of eVTOLs that charge at vertiport i (count)	Intermediate Parameter
$n$	Before eVTOL i arrives, the number of eVTOLs that arrived at vertiport j during (t − μ, t) (count)	Intermediate Parameter
$nsl$	The number of time slots while cruising (count)	Input Parameter
$O$	Oswald’s efficiency factor	Constant
$Pce$	The electricity price (CNY/kW·h)	Input Parameter
$Pa{x}_{ij}$	The number of passengers vertiport j has sent before i-th flight arrives (count)	Intermediate Parameter
$Po{b}_i$	The probability that individual i is selected	Intermediate Parameter
$Px{d}_{jk}$	Passenger demand between vertiport j and k (count)	Intermediate Parameter
$Pow$	Power (kW)	Input and Intermediate Parameter
$Pois\left(\lambda \right)$	Poisson random number with parameter λ	Function
$Pu{r}_i^P$	The purchase price of eVTOL i (CNY)	Input Parameter
$pc{l}_i$	The passenger capacity of eVTOL i (count)	Input Parameter
$pro{b}^{collision}$	The collision probability between eVTOLs	Intermediate Parameter
$R_i$	The electricity price of ith-tier	Input Parameter
$r_{ij}$	A binary variable which equals 1 if eVTOL i needs to be charged at vertiport j, 0 otherwise	Decision Variable
$rd{a}_i$	The rotor disk area of eVTOL i (m2)	Input Parameter
$ST{I}_i$	The number of eVTOLs for which takeoff or landing time dissatisfies the safety time interval of eVTOL i (count)	Intermediate Parameter
$Se{v}_{ij}$	The number of passengers served by eVTOL i in vertiport j (count)	Output Parameter
$slt$	The length of each time slot (s)	Input Parameter
$Tim$	Time duration (h)	Input and Intermediate Parameter
$Th{r}_i$	The electricity threshold of ith-tier (kW·h)	Input Parameter
$t$	Time moment	Intermediate Parameter
$VSP$	The speed of the eVTOL (km/h)	Intermediate Parameter
$Vo{l}_i$	The number of eVTOLs vertiport i could accommodate (count)	Input Parameter
$WA{R}_i$	Wing area of eVTOL i (m2)	Input Parameter
$WS{P}_i$	Wing span of eVTOL i (m)	Input Parameter
$x$	The aerial or ground distance that the eVTOL has traveled (m)	Intermediate Parameter
$\alpha$	The climb angle (°)	Intermediate Parameter
$\delta$	The safe remaining electricity level (%)	Input Parameter
$\epsilon$	The minimum percentage of demand that eVTOLs should satisfy during peak hours (%)	Input Parameter
${\eta }_{prop}$	The propulsion system efficiency (%)	Input Parameter
${\theta }_{app}$	The approach angle (°)	Input Parameter
${\theta }_{dep}$	The departure angle (°)	Input Parameter
$\mu$	The maximum waiting time that passengers could endure (h)	Input Parameter
$\rho$	The density of dry air (kg/m3)	Constant

.

2.2. Description and Assumptions of the Problem

This study resolved the scheduling problem under time-dependent demand. Firstly, the number of different types of eVTOL and the timetable are identified randomly. When an eVTOL arrives at an airport while operating according to the initialized flight schedule, the decision to charge is based on the remaining battery level and the timetable. Furthermore, additional dwelling time will also be determined after charging. Finally, the passengers served and the operational cost are calculated and optimized utilizing a genetic algorithm. Moreover, different to the previous studies [6,12], in this research passengers whose waiting time exceeds a threshold μ will choose other models of transportation instead of eVTOLs.

The hypothesis proposed in this model are summarized as follows:

(1)

There is no need for an eVTOL to return to its original vertiport after completing the flight during the entire the day.

(2)

The function of passenger demand f_i during the workday is the same with holidays.

(3)

Passengers whose waiting time exceeds μ will leave. eVTOLs with their arrival time t can only serve passengers between (t − μ, t).

(4)

An eVTOL can only be charged at vertiports and its power consumption is linearly related to the distance it traveled.

(5)

The state of motion of eVTOLs during takeoff and landing is deemed to be uniform accelerated rectilinear motion which has a constant acceleration ($ac{c}_i$).

(6)

An eVTOL will accelerate to its cruising speed ($VS{P}_i^{cr}$) after arriving at the cruising height ($h^{cr}$) and will remain at this speed.

2.3. Model and Algorithm Development

2.3.1. Model Development

Objective function:

$$\mathrm{Served} \mathrm{passenger} \mathrm{amount}: \max {\displaystyle \sum _{k=1}^{B_t}{\displaystyle \sum _{i=1}^{B_e}{\displaystyle \sum _{j=1}^{Amn}Se{v}_{ij}(t)}}}$$

(1)

$${\mathrm{Cost}: \min C}^c(t)+C^p+C^M$$

(2)

Equations (1) and (2) are the objective functions of the maximum served passengers and minimum total cost, where $B_t$ is the number of eVTOL types; $B_e$ is the number of eVTOLs belonging to each type; and $Amn$ is the number of vertiports that the eVTOL has passed. The total cost is composed of purchasing the eVTOL, charging, and maintenance.

(a)

Energy Calculation in Objective Function

As the binary variable $r_{ij}$ symbolizes whether eVTOL i will be charged at vertiport j, $N_j$, which symbolizes the number of eVTOLs that charge at vertiport j, can be calculated by summing the $r_{ij}$ of each eVTOL $B_e$ for each type $B_t$ (Equation (3)). The charging cost is the product of the total charging time $Ti{m}^c$ of all eVTOLs that charge at each vertiport (related to $N_j$), charging power $Po{w}^c$, and the electricity price $Pce(t)$ (Equation (4)). To determine whether an eVTOL could be charged, $r_{ij}$ was introduced. Herein, an eVTOL must be charged if its remaining electricity $E_{ij}$ after landing at vertiport j falls below the safe reserved electricity level ($\delta \times E_i^0$, where $E_i^0$ is the battery capacity of eVTOL i) or cannot satisfy the energy required for the next flight (Equation (5)). $E_{ij}$ the remaining electricity after landing at the present vertiport j is the sum of the supplied energy (obtained through multiplying the charging power with the charging time of eVTOL i in vertiport k) $Po{w}^c\times Ti{m}_{ik}^C$ in vertiport k, the electricity left for eVTOL i after landing at the previous vertiport k ($E_{ik}$), and the energy required for the flight between vertiport k and j ($E_{kj}^r$) as Equation (6) illustrates.

$$N_j={\displaystyle \sum _{k=1}^{B_t}{\displaystyle \sum _{i=1}^{B_e}{r}_{ij}}}$$

(3)

$$C^c(t)={\displaystyle \sum _{j=1}^A{\displaystyle \sum _{i=1}^{N_j}Ti{m}_{ij}^c\times Pce(t)}}\times Po{w}^c$$

(4)

$$r_{ij}=\{\begin{array}{l}1\hspace{1em}\hspace{1em}{E}_{ij}\le \max \left\{\delta \times E_i^0,E_{ij}^r\right\}\\ 0\hspace{1em}\hspace{1em}\mathrm{else}\end{array}$$

(5)

$$E_{ij}=E_{ik}+Po{w}^c\times Ti{m}_{ik}^C\times r_{ik}-E_{kj}^r$$

(6)

$E_{kj}^r$ is composed of the electricity required for hovering ($E_{kj}^{hv}$), climbing ($E_{kj}^{cb}$), cruising ($E_{kj}^{cr}$), and descending ($E_{kj}^{ds}$) (Equation (7)). The motion state of an eVTOL while cruising is deemed as uniform motion; therefore, the discharging power while cruising $Po{w}_{cr}^{dis}$ is also a constant. The energy required during the cruising stage is demonstrated in Equation (8) where $Di{s}_{kj}^v$ is the distance between vertiport j and k and $VS{P}_i^{cr}$ is the cruising speed. However, as the motion state while climbing and descending is uniform accelerated rectilinear motion, the power in these stages ($Po{w}_{cb}^{dis}$, $Po{w}_{ds}^{dis}$) keeps changing. Thus, integration measures have to be taken to obtain the energy consumed during the climbing and descending stages (Equations (9) and (11)). It is noteworthy that the increase in flights in the air will cause conflict between eVTOLs which may further contribute to additional hovering time to wait while cruising $Ti{m}_w$ (Equation (10)). Therefore, the energy consumed while waiting has also been considered in Equation (10). The value of the electricity price in this study was set according to the tiered electricity price, which varies on the basis of the energy consumed at different time moments t (Equation (12)).

$$E_{kj}^r=E_{kj}^{cr}+E_{kj}^{cb}+E_{kj}^{hv}+E_{kj}^{ds}$$

(7)

$$E_{kj}^{cr}=Po{w}_{cr}^{dis}\times {\displaystyle \frac{Di{s}_{kj}^v}{VS{P}_i^{cr}}}$$

(8)

$$E_{kj}^{cb}={\displaystyle \underset{T^{cb}}{\int }Po{w}_{cb}^{dis}(t)dt}$$

(9)

$$E_{kj}^{hv}=Po{w}_{hv}^{dis}\times (Ti{m}^{hv}+Ti{m}_w)$$

(10)

$$E_{kj}^{ds}={\displaystyle \underset{T^{ds}}{\int }Po{w}_{ds}^{dis}(t)dt}$$

(11)

$$Pce(t)=\{\begin{array}{l}{R}_1\hspace{1em}\hspace{1em}{\displaystyle \sum _{j=1}^A{\displaystyle \sum _{i=1}^{N_j}Ti{m}_{ij}^c\times Po{w}^c<Th{r}_1}}\\ R_2\hspace{1em}\hspace{1em}Th{r}_1\le {\displaystyle \sum _{j=1}^A{\displaystyle \sum _{i=1}^{N_j}Ti{m}_{ij}^c\times Po{w}^c}}\le Th{r}_2\\ R_3\hspace{1em}\hspace{1em}{\displaystyle \sum _{j=1}^A{\displaystyle \sum _{i=1}^{N_j}Ti{m}_{ij}^c\times Po}}{w}^c>Th{r}_2\end{array}$$

(12)

(b)

Power Calculation in Objective Function

The discharging power of cruising ($Po{w}_{cr}^{dis}$), hovering ($Po{w}_{hv}^{dis}$), climbing ($Po{w}_{cb}^{dis}$), and descending ($Po{w}_{ds}^{dis}$) in Equations (8)–(11) can be acquired through Equations (13)–(21) [17].

$$Po{w}_{cr}^{dis}={\displaystyle \frac{Forc{e}_i^D\times VS{P}_i^{cr}}{{\eta }_{prop}}}$$

(13)

$$Forc{e}_i^D={\displaystyle \frac{1}{2}}{(VS{P}_i^{cr})}^2\times \rho \times WA{R}_i\times (c_{D_0}+{\displaystyle \frac{c_L^2}{\pi \times AP{R}_i\times O}})$$

(14)

$$Po{w}_{hv}^{dis}={\displaystyle \frac{{(Forc{e}_i^T)}^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}}{FOM\times \sqrt{2\times \rho \times rda}}}$$

(15)

$$Forc{e}_i^T=m_i\times ac{c}_i+m_{i}g$$

(16)

$$ac{c}_i={\displaystyle \frac{{(VS{P}_i^{cr})}^2}{2\times (h^{hv}+x_i^{cb})}}$$

(17)

$$Po{w}_{cb}^{dis}(t)=Po{w}_{hv}^{dis}\times \left({\displaystyle \frac{VS{P}_i^{cb,v}}{2VS{P}_i^{hv}}}+\sqrt{{({\displaystyle \frac{VS{P}_i^{cb,v}}{2VS{P}_i^{hv}}})}^2+1}\right)$$

(18)

$$Po{w}_{ds}^{dis}(t)=Po{w}_{hv}^{dis}\times ({\displaystyle \frac{VS{P}_i^{ds,v}}{2VS{P}_i^{hv}}}+\sqrt{{({\displaystyle \frac{VS{P}_i^{ds,v}}{2VS{P}_i^{hv}}})}^2-1})$$

(19)

$$VS{P}_i^{cb,v}=VS{P}_i^{ds,v}=ac{c}_i\times t\times \sin \alpha$$

(20)

$$VS{P}_i^{hv}=ac{c}_i\times Ti{m}^{hv}$$

(21)

It is worth mentioning that all these equations are derived based on hypothesis (5) and Figure 1 [18]. The $VS{P}_i^{cr}$ in Equation (13) is the cruising speed of eVTOL I and ${\eta }_{prop}$ is the efficiency of the propulsion system. The drag force ($Forc{e}_i^D$) of aircraft while cruising is composed of parasitic drag $c_{D_0}$ and induced drag, while the induced drag is related to the lift coefficient $c_L$, aspect ratio APR, and Oswald’s efficiency factor (e); thus, $Forc{e}_i^D$ can be calculated in Equation (14) [17]. The hovering power can be acquired from Equations (15) and (16) where m_i is the mass of the eVTOL, ρ is the density of air, and rda is the rotor disk area. Since the whole state of motion during takeoff is uniform accelerated rectilinear motion, the thrust force $Forc{e}_i^T$ is the sum of the external force ($m_i\times ac{c}_i$) and gravity ($m_i\times g$) as illustrated in Equation (16). The acceleration ($ac{c}_i$) utilized in Equation (16) to calculate the external force can be obtained from Equation (17) where $h^{hv}$ is the hovering height and $x_i^{cb}$ is the aerial distance traveled by the eVTOL while climbing. Based on O. Ugwueze et al.’s research [17], Equations (18) and (19) were developed to calculate the climbing and descending power, which are functions related to the time moment t. The $Ti{m}^{hv}$ in Equations (18) and (19) is the hovering duration. Furthermore, $VS{P}_i^{cb,v}$ is the vertical speed of eVTOL i while climbing. As the climbing speed in the air can be acquired through multiplying $ac{c}_i$ by t, $VS{P}_i^{cb,v}$ could be obtained by multiplying the climbing speed by the sine of the climb angle α as Equations (20) and (21) illustrate.

(c)

Time and Distance Calculation in Objective Function

The distance, such as $x_i^{cb}$, and the time variables, such as $Ti{m}^{hv}$, $Ti{m}^{cb}$, etc., are related in Equations (7)–(11) and calculated in Equations (22)–(31) based on Figure 1.

$$x_i^{cb}=x_i^{ds}=x_i^{ob}+{\displaystyle \frac{(h^{cr}-h^{cb})}{\sin {\theta }_{dep}}}$$

(22)

$$x_i^{ob}={\displaystyle \frac{(h^{cb}-h^{hv})}{\sin \alpha }}$$

(23)

$$\alpha =\{\begin{array}{l}\arctan ({\displaystyle \frac{h^{cb}-h^{hv}}{x_i^g}})\hspace{1em}\hspace{1em}\mathrm{obstacle} \mathrm{surpassing}\\ {\theta }_{dep}({\theta }_{app})\hspace{1em}\hspace{1em}\mathrm{climb} \mathrm{to} \mathrm{cruising} \mathrm{altitude}\end{array}$$

(24)

$$x_i^g=WS{P}_i$$

(25)

$$Ti{m}^{hv}={\displaystyle \frac{\sqrt{2\times ac{c}_i\times h^{hv}}}{ac{c}_i}}$$

(26)

$$Ti{m}^{cb}={\displaystyle \frac{VS{P}_i^{cr}-VS{P}_i^{hv}}{ac{c}_i}}$$

(27)

$$Ti{m}^{ds}={\displaystyle \frac{VS{P}^{hv}-VS{P}_i^{cr}}{ac{c}_i}}$$

(28)

$$Ti{m}_{ij}^{cr}={\displaystyle \frac{Di{s}_{ij}^v}{VS{P}_i^{cr}}}$$

(29)

$$Di{s}_{ij}^v=6371\times 2\times (\arctan (\sqrt{a_{ik}^v})+\arctan (\sqrt{a_{kj}^v}))$$

(30)

$$a_{ik}^v=\sin {\left({\displaystyle \frac{(la{t}_k^n-la{t}_i^v}{2}}\right)}^2+\cos \left(la{t}_i^v\right)\times \cos \left(la{t}_k^n\right)\times \sin {\left({\displaystyle \frac{lo{n}_k^n-lo{n}_i^v}{2}}\right)}^2$$

(31)

$$a_{kj}^v=\sin {\left({\displaystyle \frac{(la{t}_j^v-la{t}_k^n}{2}}\right)}^2+\cos \left(la{t}_k^n\right)\times \cos \left(la{t}_j^v\right)\times \sin {\left({\displaystyle \frac{lo{n}_j^v-lo{n}_k^n}{2}}\right)}^2$$

(32)

Descend could be considered as the reverse process of climb; therefore, the aerial distance of climbing $x_i^{cb}$ and descending $x_i^{ds}$ and the angles of departure ${\theta }_{dep}$ and approach ${\theta }_{app}$ are deemed to be the same in Equation (22). The climbing distance in the air $x_i^{cb}$ consists of the aerial distance of obstacle surpassing $x_i^{ob}$ and the distance required to climb to the cruising altitude $h^{cr}$. Due to the need for the eVTOL to hover before climbing, the climbing height during the obstacle surpassing stage $h^{cb}$ should be subtracted by the hovering height $h^{hv}$ (Equation (23)). Subsequently, to obtain the aerial distance of obstacle surpassing, the value of $h^{cb}-h^{hv}$ is further divided by the sine of the climb angle α acquired through Equation (24) where $x_i^g$ is the ground distance traveled while surpassing obstacles equal to eVTOL i’s wing span (WSP_i) and could be changed based on different conditions and regulations as Equation (25) demonstrates. The distance required to climb to the cruising altitude in Equation (22) can also be calculated using the same method with its climb angle α equal to the departure angle ${\theta }_{dep}$.

As the eVTOL accelerates from $VS{P}_i^{hv}=0$, the hovering duration $Ti{m}^{hv}$ can be obtained from Equation (26). Moreover, during the climbing period the eVTOL accelerates from $VS{P}_i^{hv}$ to $VS{P}_i^{cr}$; therefore, the time for climbing $Ti{m}^{cb}$ is calculated in Equation (27). Similarly, the time for descending $Ti{m}^{ds}$ can also be calculated in Equation (28). The cruising time of eVTOL i can be acquired through dividing the distance ($Di{s}_{ij}^v$) between vertiport i and j by its cruising speed $VS{P}_i^{cr}$ (Equation (29)) where $Di{s}_{ij}^v$ is the distance between vertiport i and j. For the avoidance of injuring citizens on the ground after the collision of eVTOLs, the flying route should avoid not only prohibited zones but also densely populated areas. Suppose that $lo{n}_k^n$ and $la{t}_k^n$ are the longitude and latitude coordinates of the detour location k. Thus, the route from vertiport i to j has been split into two parts, one from vertiport i to detour point k, the other from detour point k to vertiport j, with both parts calculated by using the Haversine formulation in Equations (30)–(32) where $a^v$ is the radian difference between the vertiport and the detouring point.

Finally, the purchasing cost is the sum of the purchase price of different types of eVTOL (Equation (33)). The maintenance cost $C^M$ is calculated based on the travel distance $Di{s}_i^T$ of each eVTOL as Equation (34) demonstrates.

$$C^P={\displaystyle \sum _{j=1}^B{\displaystyle \sum _{i=1}^{b_j}Pu{r}_i^P}}$$

(33)

$$C^M={\displaystyle \sum _{j=1}^B{\displaystyle \sum _{i=1}^{b_j}Mac\times Di{s}_i^T}}$$

(34)

(d)

Constraints of Demand

As passengers can only wait for μ, the total demand that still remains at vertiport j before eVTOL i arrives is the total demand, acquired through integrating the time demand function $f_i$, which conforms to Poisson distribution as Equation (35) demonstrates [10], during the time period $t_i-\mu$ to $t_i$, minus the demand that has been already served by other eVTOLs arriving before eVTOL i (Equation (36)). Furthermore, the number of passengers eVTOL i serves at vertiport j cannot surpass the total number of passengers waiting at this vertiport ($AS{P}_{ij}$) or the eVTOL’s capacity limit $pc{l}_i$ (Equation (37)). To alleviate the ground transportation pressure to the greatest extent, the total served demand should satisfy at least ε percent of the total demand [1] (Equation (38)).

$$f_j(t)~\mathit{Pois}(\lambda )=e^{\lambda (e^{\mathit{it}}-1)}$$

(35)

$$AS{P}_{ij}(t)={\displaystyle \underset{t_{ij}-\mu }{\overset{t_{ij}}{\int }}{f}_j\left(t\right) dt}-{\displaystyle \sum _{i=1}^{n}Se{v}_{ij}}(t)$$

(36)

$$Se{v}_{ij}(t)\le \min \{pc{l}_i,AS{P}_{ij}(t)\}$$

(37)

$${\displaystyle \sum _{k=1}^B{\displaystyle \sum _{i=1}^{b_k}{\displaystyle \sum _{j=1}^{A}Se{v}_{ij}(t)}}}\ge \epsilon \times {\displaystyle \sum _j{\displaystyle \sum _{k}Px{d}_{jk}}}$$

(38)

(e)

Constraints of Energy

As it is demonstrated in J. Chen et al.’s [19] research on safety concerns, δ percent of reserved battery level is required while operating; the remaining electricity level should lie between $\delta \times E_i^0$ and its maximum capacity $E_i^0$ (Equation (39)). Therefore, the charging duration of eVTOL i in vertiport j ($Ti{m}_{ij}^C$) should not only satisfy the power required to fly from vertiport i to j but should also meet the constraint of safety regarding reserved electricity (Equation (40)). Equation (37) states that the charging time $Ti{m}_{ij}^C$ equals 0 when the eVTOL does not need charging where M is an infinite number. Moreover, the distance between vertiport i and j $Di{s}_{ij}^v$ cannot exceed the eVTOL’s range limitation (Equation (38)).

$$\delta \times E_i^0\le E_{ij}\le E_i^0$$

(39)

$$\max \left\{\delta \times E_i^0,E_{ij}^r\right\}\le Po{w}^c\times Ti{m}_{ij}^C\times r_{ij}+E_{ij}\le E_i^0$$

(40)

$$0\le Ti{m}_{ij}^C\le M\times r_{ij}$$

(41)

$$0<Di{s}_{ij}^v\le d_i$$

(42)

Equations (43) and (44) constrain the departure and arrival time, which lie between the start $Ti{m}_1$ and final service times $Ti{m}_2$ where $Ti{m}_{ij}^C$ is the charging time duration, $Ti{m}_{ij}^d$ is the dwelling time duration, and $t_{ij}^d$ is the departure time of eVTOL i at vertiport j, which is equal to $t_{ij}^a$, the arrival time, plus the charging and dwelling time at vertiport j (Equation (45)). Herein, $t_{ij}^a$ is the sum of the arrival time of eVTOL i at its previous vertiport k ($t_{ik}^a$), the charging time $Ti{m}_{ik}^C$, dwelling time $Ti{m}_{ik}^d$ at vertiport k, and the hovering, climbing, cruising, and descending times between the two vertiports (Equation (46)). It is worth mentioning that since hovering is included in both the takeoff and landing periods, it is calculated twice. Moreover, for safety considerations the safety time interval $Ti{m}_s$ while taking off and landing and the waiting time $Ti{m}_w$ while cruising if there is a conflict between eVTOL i and another eVTOL are also considered in Equation (46). As it is demonstrated in hypothesis (1) that each vertiport could accommodate at most 4 eVTOLs taking off and landing at the same time, if the number of eVTOLs that dissatisfy the safety interval constraints of Equation (43) or Equation (44) is more than 4, all of these eVTOLs that dissatisfy constraints will be delayed for $\left[{\displaystyle \frac{ST{I}_i}{4}}\right]Ti{m}_s$.

$$Ti{m}_1\le t_{ij}^a\le Ti{m}_2$$

(43)

$$Ti{m}_1\le t_{ij}^d\le Ti{m}_2$$

(44)

$$t_{ij}^d=t_{ij}^a+Ti{m}_{ij}^C+Ti{m}_{ij}^d$$

(45)

$$t_{ij}^a=t_{ik}^a+Ti{m}_{ik}^C+Ti{m}_{ik}^d+Ti{m}_{kj}^{cr}+2Ti{m}^{hv}+Ti{m}^{cb}+Ti{m}^{ds}+\left[{\displaystyle \frac{ST{I}_i}{4}}\right]Ti{m}_s+Ti{m}_w$$

(46)

(f)

Constraints of Safety

To ensure the safety of eVTOL i and k, which take off or land at vertiport j continuously, the minimum safety time interval $Ti{m}_s$ has been introduced to constrain their departure $t_{ij}^d$,$t_{kj}^d$ or arrival times $t_{ij}^a$, $t_{kj}^a$ (Equation (47) and Equation (48)). The safety intervals during takeoff and landing are considered to be equal and are obtained based on Figure 1 [18] and Equation(49). According to the document published by the EASA, eVTOLs are deemed to have departed from a vertiport after reaching a certain altitude. Therefore, in this article, an eVTOL could take off or land after the last eVTOL has totally departed or approached, which indicates that the safety time interval $Ti{m}_s$ is composed of the time needed to hover and to climb to the predetermined departure or approach altitude (Equation (49)), while the time required to climb to the predetermined altitude is further divided into obstacle surpassing and the climbing stage from obstacle surpassing to the departure or approach altitude (Equation (49)). In order to calculate the time of this stage, the aerial distance $x^s$ of this period is needed. In this research, $x^s$ could be calculated utilizing the division between $h^s-h^{cb}$ and $\sin {\theta }_{dep}$ based on Figure 1 as Equation (50) demonstrates; the same goes for the approach.

$$t_{kj}^d-t_{ij}^d\ge Ti{m}_s$$

(47)

$$t_{kj}^a-t_{ij}^a\ge Ti{m}_s$$

(48)

$$Ti{m}_s=Ti{m}^{hv}+\frac{\sqrt{{(VS{P}_i^{hv})}^2+2\times ac{c}_i\times (x_i^{ob}+x^s)}}{ac{c}_i}$$

(49)

$$x^s={\displaystyle \frac{h^s-h^{cb}}{\sin {\theta }_{dep}}}$$

(50)

To identify the value of $Ti{m}_w$ Equations (51)–(58) have been established based on Lili W. [20] et al.’s research. Since $Ti{m}_w$ has a close relation with the recognition of conflicts which depends on the value of the safety interval distance $L_s$, $L_s$ should firstly be determined on the basis of Equations (51)–(55). Equation (51) calculates the collision probability between eVTOL i and j where $VS{P}_{ij}^{cr,y}$ is the relative velocity between eVTOL i and j; Dev is the deviation of the CNS system—Dev₁, Dev₂, and Dev₃ are the deviations of the CNS systems of RNPn1, RCPn2, and RSPn3, correspondingly (Equation (52)); and Ind is the initial distance between two eVTOLs. The time period of eVTOL i and j while cruising has been divided into nsl number of time slots each with slt seconds (Equation (53)). In order to ensure the safety of the eVTOL while cruising, the probability of a collision obtained from Equation (51) must be 0 as Equation (55) demonstrates, and thus the value of the safety interval distance $L_s$ could be determined. For eVTOLs whose distance between each other $Di{s}_{ij}^e$ lies below the safety interval $L_s$, one should wait for the other until their distance meets the safety interval constraint (Equation (56)). The distance between eVTOLs $Di{s}_{ij}^e$ is calculated with the same method in Equation (30) and Equation (31) where $lo{n}_i^e$ and $la{t}_i^e$ are the longitude and latitude of the eVTOL (Equations (53) and (54)).

$$pro{b}^{collision}(L_s)={\displaystyle \frac{{\displaystyle \underset{-L_s}{\overset{L_s}{\int }}{e}^{-\frac{{(y-Ind)}^3}{4\times nsl\times 0.2603\times Dev}}}dy}{3\sqrt{nsl\times \pi \times 0.2603\times Dev}}}$$

(51)

$$Dev=De{v}_1^2+De{v}_2^2\times {(VS{P}_{ij}^{cr,y})}^2+De{v}_3^2\times {(VS{P}_{ij}^{cr,y})}^2$$

(52)

$$nsl={\displaystyle \frac{Ti{m}_{ij}^{cr}}{slt}}$$

(53)

$$VS{P}_{ij}^{cr,y}=VS{P}_i^{cr}-VS{P}_j^{cr}$$

(54)

$$pro{b}^{collision}(L_s)=0$$

(55)

$$Ti{m}_w=\{\begin{array}{l}{\displaystyle \frac{L_s}{V_i^{cr}}}\hspace{1em}\hspace{1em}Di{s}_{ij}^e<L_s\\ 0\hspace{1em}\hspace{1em}\hspace{1em}\mathrm{Otherwise}\end{array}$$

(56)

$$Di{s}_{ij}^e=6371\times 2\times \arctan (\sqrt{a_{ij}^e})$$

(57)

$$a_{ij}^e=\sin {\left({\displaystyle \frac{(la{t}_j^e-la{t}_i^e}{2}}\right)}^2+\cos \left(la{t}_i^e\right)\times \cos \left(la{t}_j^e\right)\times \sin {\left({\displaystyle \frac{lo{n}_j^e-lo{n}_i^e}{2}}\right)}^2$$

(58)

Furthermore, to ensure the safety of ground citizens and avoid the prohibited areas, the detour point should be $Di{s}^m$ km away from these areas (Equation (59)). Suppose that ($lo{n}^p,la{t}^p$) is the lon–lat coordinates of populated or prohibited areas and ($lo{n}_k^n,la{t}_k^n$) is the lon–lat coordinates of the detouring point, then the distance between the detour points and prohibited areas $Di{s}_{jk}^d$ is calculated in Equations (59)–(61).

$$Di{s}_{jk}^d\ge Di{s}^m$$

(59)

$$Di{s}_{jk}^d=6371\times 2\times \arctan (\sqrt{a_{jk}^d})$$

(60)

$$a_{jk}^d=\sin {\left({\displaystyle \frac{(la{t}^p-la{t}_k^n}{2}}\right)}^2+\cos \left(la{t}_k^n\right)\times \cos \left(la{t}^p\right)\times \sin {\left({\displaystyle \frac{lo{n}^p-lo{n}_k^n}{2}}\right)}^2$$

(61)

Equation (62) limits the dwelling time of each eVTOL. Finally, the cruising altitude should not exceed the eVTOL’s maximum cruising altitude ($MH{A}_i$) as Equation (63) demonstrates.

$$0\le Ti{m}_{ij}^d\le 1$$

(62)

$$h^{cr}\le MH{A}_i$$

(63)

2.3.2. Algorithm Design

In this research, a Joint Optimization of Cost and Schedule Particle Swarm Optimization algorithm (JOCS-PSO) has been utilized to solve this joint routing and charging model with safety concerns. The detailed process of this algorithm is demonstrated in Algorithm 1. The JOCS-PSO algorithm utilized in this research consists of two parts where the Cost optimization part is used to initialize the configuration of different types of eVTOL and has the total cost optimized under the schedule acquired from the inner algorithm, while the Scheduling Part is used to maximize the served passenger amount through iteratively updating the schedule on the basis of this eVTOL configuration acquired in the Cost optimization part. The relationship between the Cost and Schedule Parts is demonstrated in Figure 2. To simplify calculations, the total cost is converted into a negative value to seek its maximum.

In the Cost Part, each individual containing B (the number of types of eVTOL) variables has been initialized and will be transferred to the inner layer. In the Scheduling Part employed in this study, each individual consists of three characters. The first character represents the sequence of vertiports that each eVTOL passes in time order, with 0 used as a separator between different eVTOLs. The second indicates the charging duration for each eVTOL at the corresponding vertiports defined in the first character. The third character denotes the additional dwelling time for each eVTOL at the corresponding vertiports as Figure 3 illustrates. The whole schedule of these eVTOLs is made up of these three characters.

The aforementioned characters in Figure 3 mean that the sequence of vertiports for the first eVTOL, in chronological order, is A, B, C, and D. The charging time in each vertiport is 0.18 h, 0,32 h, 0.05 h, and 0.19 h, correspondingly, and their dwelling times are 0.02 h, 0.11 h, 0.13 h, and 0.11 h. After the inner layer has reached its iteration times, the schedule with the maximum served demand will be transferred back to the outer layer to calculate the total cost.

After all the iteration times of cost and scheduling optimization were reached, the served passenger amount and the total cost under this eVTOL configuration and schedule was evaluated and optimized at the same time using Pareto dominance.

Algorithm 1. Joint optimization algorithm of cost and schedule

$\mathrm{Input}: \mathrm{distance} \mathrm{between} \mathrm{vertiports} Di{s}_{ij}^v$$, \mathrm{Safe} \mathrm{interval} Ti{m}_s$$, \mathrm{amount} \mathrm{of} \mathrm{eVTOL} \mathrm{types} B, \mathrm{charging} \mathrm{power} Po{w}^c$$, \mathrm{eVTOL} \mathrm{parameter}, \mathrm{Operation} \mathrm{time} Ti{m}_1$$, Ti{m}_2$

For i ∈ iteration times of outer layer do

Initialize the amount of different types of eVTOL

For i ∈ iteration times of inner layer do

Initialize the schedule (Figure 3) based on the eVTOL configuration

If onstraints Equations (35)–(63) are dissatisfied

Rectify the characters of this schedule

end

Sum the served demand of all the eVTOL based on Equation (1)

If  total served demand > max served demand

Best schedule = present schedule

max served demand = total served demand

end

Renew the schedule

end

If  max served demand > ε × Pxd

Calculate the cost based on Equations (3), (33) and (34) and the schedule

else

Total cost = inf

end

If  Total cost < min cost

min cost = Total cost

end

Renew the amount of different types of eVTOL

end

Output: Number of different types of eVTOL, timetable and charging demand

Initialization: The initialization method of the eVTOL configuration is performed through generating the number of each type of eVTOL randomly, while the initialization method of the schedule is performed through simulating the real-world operation of eVTOL under the constraints of Equations (35)–(63). The pseudo code is demonstrated in Algorithm 2.

Algorithm 2. Initialization of schedule

Input: distance between vertiports $Di{s}_{ij}^v$, Safe interval $Ti{m}_s$, amount of eVTOL types B, charging power $Po{w}^c$, eVTOL parameter, Operation time $Ti{m}_1$, $Ti{m}_2$

For i = 1: the total amount of eVTOL

$t_{ij}^d=t_{ij}^a=Ti{m}_1$ (Set the original depart time)

While $t_{ij}^a\le Ti{m}_2$

Select the initial vertiport among A, B, C, D, E, F randomly

Select the next vertiport among A, B, C, D, E, F randomly whose distance between the initial vertiport satisfies with the distance constraints $0<Di{s}_{ij}^v\le d_i$ and the energy constraint $\delta \times E_i^0\le E_{ij}\le E_i^0$

If $E_{ij}-E_{ij}^r<\delta \times E_i^0$

Charge eVTOL until it meets $\max \left\{\delta \times E_i^0,E_{ij}^r\right\}\le Po{w}^c\times Ti{m}_{ij}^C\times r_{ij}+E_{ij}\le E_i^0$

End if

Generate $Ti{m}_{ij}^d$ within $0\le Ti{m}_{ij}^d\le 1$

$t_{ij}^d=Ti{m}_1+Ti{m}_{ij}^C+Ti{m}_{ij}^d$

Find the amount of eVTOL that dissatisfies with $t_{kj}^d-t_{ij}^d\ge Ti{m}_s$

While the amount > 1

$t_{ij}^d=Ti{m}_1+Ti{m}_{ij}^C+Ti{m}_{ij}^d+(amount-1)\times Ti{m}_s$

Find the amount of eVTOL that dissatisfies with $t_{kj}^d-t_{ij}^d\ge Ti{m}_s$ again

End While

$t_{ik}^a=t_{ij}^d+Ti{m}_{jk}^{cr}+2Ti{m}^{hv}+Ti{m}^{cb}+Ti{m}^{ds}$

Write the passing vertiports, $Ti{m}_{ij}^d$, $Ti{m}_{ij}^C$ into the schedule as Figure 3

End while

End for

Renewing: To expand the sample size and the search space of this self-developed algorithm, several flight schedule schemes have been initialized with each as an individual. For the aim of accelerating convergence, searching speed $spd$, individual $idv$, and social learning $scl$ factors have been introduced to renew the eVTOL configuration and the schedule. The method for updating flight schedules employed by the cost and scheduling joint optimization algorithm developed in this study is demonstrated in Equations (64) and (65) where cfs^k is the schedule in iteration time k and ω is the inertia weight. r₁ and r₂ are two matrixes that are randomly generated; pBest and gBest are the best solution of each individual and the whole population. The renewal plan of eVTOL configuration is also as follows, with two numbers r₁ and r₂ generated randomly within 0 to 1.

$$cf{s}^{k+1}=cf{s}^k+sp{d}^{k+1}$$

(64)

$$sp{d}^{k+1}=\varpi \times sp{d}^k+idv\times r_1\times (pBest-cf{s}^k)+scl\times r_2\times (gBest-cf{s}^k)$$

(65)

The detailed values of candidate flight schedule amount, iteration time, inertia weight, individual $idv$, and social learning $scl$ factors are displayed in

Table 2. Algorithm parameters of schedule and cost optimization algorithm.

Algorithm	Parameter	Value
Cost Optimization Part	Iteration time	50
	Candidate flight schedule amount	10
	Inertia weight	0.8
	Individual learning factor	1.5
	Social learning factor	1.5
Scheduling Optimization Part	Iteration time	50
	Candidate eVTOL configuration amount	10
	Inertia weight	0.5
	Individual learning factor	1
	Social learning factor	1

.

3. Case Study

Beijing–Tianjin–Xiong’an (Hebei) was taken as a case study in this research which has a huge amount of demand due to the congestion of ground transportation as far as we are concerned [21]. Two types of eVTOL (B = 2) with different capacities (Xpeng X2 Guangzhou, China and Geely Aerofugla AE200, manufacturer, Chengdu, China) that have successfully achieved airworthiness were selected for service in this study.

3.1. Settings of Basic Parameters

3.1.1. Settings of Route and Airspace

Vertiports were selected based on the demand in each different city, with three vertiports in Beijing, two in Tianjin, and one in Xiong’an (Hebei) as shown in Figure 4 [21]. For the safety concerns of citizens and the protection of government and military facilities, appropriate detouring measures should be undertaken while cruising as illustrated in Figure 5a where the distance between the route and these areas should be over 10 km according to the regulations of the CAAC [22]. The red area in Figure 5a represents the prohibited area; the yellow area represents the densely populated region. The whole flight process is demonstrated in Figure 5b. The distance $D_{ij}$ between vertiports considering detouring measures is illustrated in

Table 3. The distance between vertiports.

	A	B	C	D	E	F
A	0	30.20	142.74	132.05	134.44	120.98
B	30.20	0	112.89	102.93	107.59	96.79
C	142.74	112.89	0	17.39	40.12	92.69
D	132.05	102.93	17.39	0	24.64	99.04
E	134.44	107.59	40.12	24.64	0	121.26
F	120.98	96.79	92.69	99.04	121.26	0

. The takeoff and landing procedure parameters in Figure 1 are shown in

Table 4. The takeoff and landing procedure parameters.

Variables	Value
$h^{cb}$	30.5 m
$h^{cr}$	300–1000 m
$h^{ds}$	30.5 m
$h^{hv}$	3 m
$h^s$	152 m
${\theta }_{app}$	7.125°
${\theta }_{dep}$	7.125°

[18]. The variable $x_i^g$ in Figure 1 equals the wing span $W{S}_i$ of the eVTOL. The departure and approach point utilized to calculate the safety time interval $Ti{m}_s$ in this research is 152 m above the vertiport, which is the end point of the obstacle limitation surface (OLS) and the takeoff climb/approach surface according to the EASA [17]. Furthermore, based on the report by the EASA [17] and the wing span of the eVTOL used in this research, the maximum area required for the eVTOL to take off and land is 900 m², while, according to the research conducted by Y. Wang [18], the vertiports in Beijing–Tianjin–Xiong’an are at least 4200 m². Thus, the vertiport can accommodate at least four pads with a sufficient distance interval between each other. Owing to the sufficient distance interval, the vortices generated by eVTOLs during simultaneous takeoff and landing on these four platforms will not affect each other. Therefore, the vertiport can accommodate four eVTOLs taking off and landing at the same time.

Moreover, in order to reduce congestion and avoid collision or conflicts between different types and directions of eVTOL, the cruising altitude $h^{cr}$ of the eVTOL in

Table 5. The altitude assignment to different eVTOLs.

eVTOL Type	Altitude Layer (m)	Flying Direction
Geely Aerofugla AE200 X01	750–1000	earlier alphabetical codes to later (A–F)
	500–750	later alphabetical codes to earlier (F–A)
Xpeng X2	400–500	earlier alphabetical codes to later (A–F)
	300–400	later alphabetical codes to earlier (F–A)

, acquired on the basis of the CAAC [23] and the performance of the eVTOL [24,25], is divided into four altitude layers according to the cruising altitude of these two eVTOLs. Two of the layers with a higher altitude (500–1000 m) are assigned to Geely Aerofugla AE200 X01 for cruising as it has a higher speed than Xpeng X2. One of these two layers was allocated to flights flying from vertiports with earlier alphabetical codes to later alphabetical codes, while the other layer is the opposite. The assignment of altitude to Xpeng X2 can be obtained similarly. The allocation of altitude is illustrated in Table 5 in detail. Furthermore, each altitude layer is deemed to have the ability to accommodate 20 eVTOLs at the same time; namely, the airspace is enough for all the eVTOLs to keep the safety interval distance at the same time while cruising. Thus, the waiting time in the air $T_w$ is 0 ($T_w=0$) in this case study. Moreover, it is worth mentioning that the energy consumption between vertiports in Equations (7)–(11) have also been calculated on the basis of Table 5.

3.1.2. Settings of eVTOL-Related Parameters

Two types of eVTOL (B = 2) with different capacities that are airworthiness were selected for service in this study. The basic parameters of these two kinds of eVTOL are displayed in

Table 6. The parameters of different eVTOLs.

Parameter	Xpeng X2	Geely Aerofugla AE200 X01
$\mathrm{Battery} \mathrm{capacity}{E}^0$ (kw·h)	120	250
$\mathrm{Purchase} \mathrm{price} C^P$ (million CNY)	0.9	2
$\mathrm{Lift} \mathrm{coefficient} c_L$	1.5	1.5
$\mathrm{Wing} \mathrm{aspect} \mathrm{ratio} A{R}_i$	7.0	7.0
$\mathrm{Wing} \mathrm{span} W{S}_i$ (m)	4.79	14.5
$\mathrm{Wing} \mathrm{area} W{A}_i$ (m2)	22.89	60.10
Rotor disk area $rd{a}_i$ (m2)	61.80	435.50
$\mathrm{Maximum} \mathrm{takeoff} \mathrm{mass} m_i$ (kg)	760	2500
$\mathrm{Figure} \mathrm{of} \mathrm{merit} FM$	0.75	0.75
$\mathrm{Maximum} \mathrm{cruising} \mathrm{altitude} M{H}_i$ (m)	500	1000
Speed $V$ (km/h)	130	264
Oswald’s efficiency factor $O$	0.85	0.85
Passenger capacity p (count)	2	5
Range $d$ (km)	75	200
$\mathrm{Propulsion} \mathrm{system} \mathrm{efficiency} {\eta }_{prop}$	0.85	0.85

[24,25].

According to Thipphavong D P et al.’s [26] research, the charging power of eVTOLs is between 200–600 kW; thus, in this article the charging power P is set as 200 kW (P = 200 kW). Furthermore, the remaining electricity level in the battery must be higher than 30% (δ = 30%) according to J. Chen et al.’s [19] research for safety concerns. The electricity price $P^e\left(t\right)$ was determined based on the Beijing Non-Residential Electricity Charging Standards [27] (Equation (66)).

$$P^e(t)=\{\begin{array}{l}0.48 \hspace{1em}\hspace{1em}{\displaystyle \sum _t{\displaystyle \sum _{i=1}^{A_t}{\displaystyle \sum _{j=1}^{N_i}{T}_{ij}^c\times Po}}\le 240}\\ 0.53 \hspace{1em}\hspace{1em}241\le {\displaystyle \sum _t{\displaystyle \sum _{i=1}^{A_t}{\displaystyle \sum _{j=1}^{N_i}{T}_{ij}^c\times Po}}\le 400}\\ 0.78\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}{\displaystyle \sum _t{\displaystyle \sum _{i=1}^{A_t}{\displaystyle \sum _{j=1}^{N_i}{T}_{ij}^c\times Po}}\ge 401}\end{array}$$

(66)

3.1.3. Settings of Operational Control

While considering daily eVTOL scheduling problems, it is necessary to ascertain the operating time and the relationship between passenger flow and time variation. As demonstrated by Dagi G. et al. [28], the flight rule for UAM is mainly the Visual Flight Rule (VFR), which results in the fact that eVTOLs cannot operate during the evening or night owing to their low visibility. Therefore, in this research, the operational time of eVTOLs is from 6:30 a.m. to 17:30 p.m., which indicates $T_1$ = 6.5 and $T_2$ = 17.5. Part of the passenger demand between vertiports at different times during this period (6:30–17:30) was determined based on Yihui W et al.’s [29] study and is shown in Figure 6. Furthermore, it is worth mentioning that as eVTOLs can only operate from 6:30 a.m. to 17:30 p.m., the passenger demand was also calculated in this time period. Moreover, in order to provide passengers with a better travel experience and alleviate ground transportation demand to the greatest extent, the maximum waiting time for passengers is 0.15 h (μ = 0.15 h) and the total number of passengers served by eVTOL services must not be less than 45% of the total demand (ε = 45%). Lastly, to ensure the safety of eVTOLs during takeoff and landing, the safe interval T_s is set as 5 min according to the document published by the CAAC [16]. To simplify the calculation, the maintenance cost $C^M$ is deemed to be 20% [28] of the purchasing cost $C^P$.

3.2. Economic Efficiency Analysis

The programming and solving process in this study was conducted according to the above-mentioned passengers’ demands and parameter setting based on MATLAB 2019a by using a computer with 16 GB memory and AMD Ryzen 7 5700X 3.40 GHz CPU.

The convergence curve of PSO in optimizing the cost and served demand utilizing Pareto dominance is displayed in Figure 7. As the principle of Pareto dominance can only identify solutions that have at least one objective function value better than other solutions, the final Pareto frontiers may have more than one solution as Figure 7 illustrates. It can be observed from Figure 7 that the total cost increases with the total served demand. This can be attributed to the fact that as the served demand increases, more eVTOLs will be needed. Moreover, to evaluate the cost effectiveness of eVTOLs, in this article, some of these Pareto frontiers were selected as a control group. The total cost and served demand of different configurations of eVTOLs are illustrated in

Table 7. Cost and served demand of different eVTOL configurations.

Configuration (Count)	First Operation Day Cost (CNY Hundred Million)	Total Cost During Lifecycle (CNY Hundred Million)	First Operation Day Served Demand (Count)	Total Served Demand During Lifecycle (Count)	Service Cost per Passenger (CNY)
166 Xpeng X2 + 171 Geely Aerofugla AE200	5.95	344.35	22,682	124,183,950	277.29
264 Xpeng X2 + 268 Geely Aerofugla AE200	9.56	555.39	33,151	181,501,725	305.71

.

The total cost during an eVTOL’s entire lifecycle C^A in Table 7 is calculated on the basis of Equation (67). As the lifetime of an eVTOL is estimated to be 15 years [30], the purchase cost $C^P$ is calculated once. However, due to the daily maintenance and charging of eVTOLs, the maintenance $C^M$ and charging costs $C^c$ should be multiplied by 15 years. Moreover, as the passenger demand is deemed to be unchanged on workdays and holidays according to hypothesis (2) the annual served demand can be acquired through multiplying the daily served demand by 15 years.

$$C^A=C^P+15\times 365\times (C^M+C^c)$$

(67)

It can be concluded from Table 7 that as the number of eVTOLs increases, the cost and served demand increase sharply, an increase of 66.3% in total cost and 46.2% in daily served demand. Moreover, the extra cost of serving each additional person can also be obtained on the basis of Table 7, utilizing the total cost during an eVTOL’s lifecycle divided by the served demand during its lifecycle. The extra cost of serving each additional person is around CNY 277.29 and 305.71 under these two eVTOL configurations, correspondingly. Through incurring a little additional expense, one can save 1 to 1.5 h, which might be the best choice for people who are in a hurry especially in the morning rush hours. Moreover, to save cost, the eVTOL configuration with 166 Xpeng X2 and 177 Geely Aerofugla AE200 will be further analyzed below.

3.3. Results

Owing to the reasons mentioned above, the optimal number of eVTOLs is 337 with 171 Geely Aerofugla AE200 X01 and 166 Xpeng X2 with a total cost of CNY 5.95 hundred million in the first operation day. Furthermore, the total served passengers are 22,682 under this scheduling plan, which represents 46% of the total demand of 49,308, which proves that the ground transportation pressure can be greatly relieved especially during morning rush hours under this eVTOL configuration. The scheduling timetable of part of the eVTOLs is illustrated in Figure 8.

As has been proved by some researchers [6,7], to serve more demand, passengers usually have to wait a relatively long time. Therefore, the relationship between tolerance waiting time and served demand is not discussed here. Under the condition that the passenger waiting time does not exceed 0.15 h, it can be concluded from Figure 8 that the operational time limitation has not been violated by any eVTOL. Meanwhile, for a single eVTOL, it can perform at least 15 flights (passing 15 vertiports) within a day. Furthermore, owing to the distance between vertiports and the range of Xpeng X2, Xpeng X2 can only operate between vertiports within the same city (Figure 8). However, as the speed of Xpeng X2 is relatively high, the number of vertiports it visits within a single day is higher than Geely Aerofugla AE200.

Furthermore, the additional dwell time apart from charging time at any vertiport for any eVTOL does not exceed 0.2 h, causing the utilization rate to exceed 77%, which indicates that most the operational time of the eVTOL has been fully used. Moreover, the charging time of each eVTOL which needs to be charged accounts for over 80% of the total time resident in each vertiport. This could be attributed to the fact that after being fully charged, prolonged extra dwelling time is unnecessary and will only lead to a decrease in the amount of service demand and the loss of passengers.

The charging demand to support the operation of eVTOLs has also been evaluated in this research (

Table 8. Charging demand of each vertiport.

Vertiport	Total Energy Required (MW·h)	Charging Cost (CNY Thousand)
A	136.32	66.80
B	945.12	519.82
C	2704.32	1757.81
D	2406.24	1515.93
E	1028.16	596.33
F	2686.56	1692.53

). It can be observed from Table 8 that the charging cost at C, D, and F vertiports is much higher than any other vertiports. This can be attributed to their high demand in C and D, located within the Fifth Ring of Beijing, and the transportation hubs in these areas such as the airports in F [20].

3.4. Safety Margin Analysis

The safety interval, as an input parameter to reflect the safety margin, is important in considering the scheduling problem. A higher safety interval will contribute to a higher safety margin especially during takeoff and landing. However, it will also result in a low number of service passengers and low service quality owing to the long flight delay. On the other hand, a lower safety interval will sacrifice safety for a higher number of service passengers and higher service quality. Therefore, it is necessary to balance service and safety in this eVTOL scheduling problem.

The number of different types of eVTOL needed to serve at least 45% of the total demand and the total cost under different safety intervals (1 min, 3 min, 6 min) has been calculated to determine the acceptable safety margin range. It is worth mentioning that according to the document produced by the EASA [18], the takeoff decision point (TDP) can be placed 30.5 m above the vertiport, which is the end point of the obstacle surpassing stage. An eVTOL could choose to stop takeoff at any time before it reaches the TDP and land back at the vertiport. Therefore, the safety time interval selected must not be less than the time it takes for the eVTOL to reach the TDP, which is 1 min. The results are given in

Table 9. Number of eVTOLs and total cost needed under different safety time intervals.

Safety Time Interval (Min)	eVTOL Amount (Count)	Cost During Lifecycle (CNY Hundred Million)	Served Demand During Lifecycle (Count)	Service Cost per Passenger (CNY)
1	165 Xpeng X2 + 204 Geely	368.57	138,610,575	265.9
3	146 Xpeng X2 + 178 Geely	367.65	136,360,350	269.62
6	166 Xpeng X2 + 171 Geely	344.35	124,183,950	277.29

.

Table 9 shows the number of eVTOLs required to increase the safety time interval. However, the served demand decreases. This can be attributed to the fact that the longer the safety interval, the longer the eVTOL has been idle and the greater the passenger loss; thus, to satisfy the 45% constraints, more eVTOLs will be needed. Moreover, despite the increase in the number of eVTOLs, the total cost is still reduced. This might be attributed to the fact that owing to the increased safety intervals, a large number of eVTOLs are idle for the majority of the day and thus do not incur charging and maintenance costs. Furthermore, it can also be observed that the service cost per passenger increases with the safety time interval. This is due to the reduction in the total number of eVTOLs, which lowers the acquisition costs that account for a significant portion of the total costs, while increasing the number of passengers served.

4. Conclusions

This study developed a joint optimization model of scheduling and charging to maximize the revenue and the number of served passengers on the basis of six vertiports in Beijing–Tianjin–Xiong’an (Hebei). Two kinds of eVTOL, Xpeng X2 and Geely Aerofugla AE200 X01, were utilized to serve the real-time demand. Moreover, safety factors such as the time interval between eVTOLs during takeoff and landing have also been considered. Furthermore, to provide passengers with a better travel experience, a maximum waiting time has also been introduced. The proposed model was resolved using a JOCS-PSO algorithm. The conclusions below can be acquired based on this research.

(1)

As the safety time interval increases from 1 min to 6 min, the number of eVTOLs needed and the service cost per passenger increase by 4% and 3%, while the total cost and served demand during the lifespan of the eVTOL decrease by 7% and 10%. This can be attributed to the longer idling time of eVTOLs, which may lead to the greater loss of passengers. In conclusion, the increase in the safety time interval will cause a lower utilization rate of eVTOLs and the decrease in the total cost and served demand and thus increase the service cost per passenger. Therefore, as the safety time increases, a higher ticket price will be needed to prevent deficiency.

(2)

Under the condition of 1 min safety time interval, the optimal number of these two types of eVTOL is 165 of Xpeng X2 and 204 of Geely Aerofugla AE200, correspondingly, with the total operation cost of CNY 368.57 hundred million during the entire lifespan of these eVTOLs. The total served demand of these eVTOLs during their entire life is 138,610,575 with a CNY 265.9 per capita cost, slightly higher than the economic cost of driving from Beijing to Tianjin. However, when compared with driving, eVTOLs can save 1–1.5 h, which is necessary for passengers in a hurry especially in rush hours. Therefore, there is a certain priority for passengers in rush hours to select eVTOLs as their transportation mode. Moreover, the total energy required within a day in Beijing, Tianjin, and Xiong’an is 398.72 MW·h, 60.52 MW·h, 52.58 MW·h, correspondingly, and could be satisfied owing to the strong power supply capability.

(3)

The number of flights an eVTOL could execute under this scheduling plan from 6:30 a.m. to 17:30 p.m. is at least 15 and the number of vertiports a short-range eVTOL passes through during the day is more than a long-range eVTOL. Moreover, the utilization rate of a single eVTOL is over 77%; no safety constraints have been dissatisfied. Furthermore, the charging time of each eVTOL that needs to be charged accounts for over 80% percent of the total time resident in each vertiport. This indicates that each eVTOL does not stay for too long at the vertiport except for charging. In conclusion, most the operational time of the eVTOL has been used to serve passengers safely and thus proves the effectiveness of this eVTOL schedule and the proposed algorithm.