Optimal Tilt-Wing eVTOL Takeoff Trajectory Prediction Using Regression Generative Adversarial Networks

Abstract

Electric vertical takeoff and landing (eVTOL) aircraft have attracted tremendous attention nowadays due to their flexible maneuverability, precise control, cost efficiency, and low noise. The optimal takeoff trajectory design is a key component of cost-effective and passenger-friendly eVTOL systems. However, conventional design optimization is typically computationally prohibitive due to the adoption of high-fidelity simulation models in an iterative manner. Machine learning (ML) allows rapid decision making; however, new ML surrogate modeling architectures and strategies are still desired to address large-scale problems. Therefore, we showcase a novel regression generative adversarial network (regGAN) surrogate for fast interactive optimal takeoff trajectory predictions of eVTOL aircraft. The regGAN leverages generative adversarial network architectures for regression tasks with a combined loss function of a mean squared error (MSE) loss and an adversarial binary cross-entropy (BC) loss. Moreover, we introduce a surrogate-based inverse mapping concept into eVTOL optimal trajectory designs for the first time. In particular, an inverse-mapping surrogate takes design requirements (including design constraints and flight condition parameters) as input and directly predicts optimal trajectory designs, with no need to run design optimizations once trained. We demonstrated the regGAN on optimal takeoff trajectory designs for the Airbus $A^3$ Vahana. The results revealed that regGAN outperformed reference surrogate strategies, including multi-output Gaussian processes and conditional generative adversarial network surrogates, by matching simulation-based ground truth with 99.6% relative testing accuracy using 1000 training samples. A parametric study showed that a regGAN surrogate with an MSE weight of one and a BC weight of 0.01 consistently achieved over 99.5% accuracy (denoting negligible predictive errors) using 400 training samples, while other regGAN models require at least 800 samples.

1. Introduction

The attention on electric vertical takeoff and landing (eVTOL) aircraft has grown significantly because eVTOL is suitable for missions that require flexible maneuverability and precise control, such as urban air mobility (UAM) [1,2]. UAM represents a safe and efficient air transportation system where everything from small package delivery drones using unmanned aerial vehicles (UAVs) to passenger-carrying air taxis through eVTOL aircraft can operate above populated areas [3,4]. UAVs have broad applications, including the sixth-generation communication using UAVs and an aerial base station for supporting the Internet of Things deployment in remote and disaster areas [5]. Hua et al. [6] proposed 3D non-stationary modeling for UAV-to-ground communications and managed to demonstrate its validity and practicality against measured results. eVTOL aircraft’s unique characteristics of precise delivery, lower cost, and reduced noise have motivated significant developments [7,8]. Conceptions under development can be categorized into a few major aircraft types, including lift+cruise (such as Aurora Flight Sciences eVTOL [9]) and tilt-wing (such as Airbus ${\mathrm{A}}^3$ Vahana [10]).

Among the aforementioned eVTOL types, the tilt-wing configuration enables the aircraft to combine the flexibility of a helicopter for vertical takeoff and landing with the efficiency of airplanes during cruising. Thus, the transition optimization on tilt-wing eVTOL aircraft attracts special interest since transition has major effects on the success of flight tasks [11]. Pradeep and Wei [12,13] put an emphasis on the formulation of multiphase optimal control problems with energy consumption index for tilt-wing and multirotor eVTOL vehicles. The proposed multiphase optimal control problem formulation and the numerical solution allowed an eVTOL air taxi to fulfill the specified required arrival time while attaining the most energy-efficient trajectory for arrival. This capability played a vital role in enabling safe and efficient future operations of eVTOL aircraft, facilitating passenger transportation and cargo delivery. Chauhan and Martins [14] constructed an Airbus ${\mathrm{A}}^3$ Vahana [10] model and optimized its takeoff-to-cruise trajectory with the objective of minimizing energy consumption. They concluded that the optimal takeoff trajectory involved stalling the wings or flying near the stall angle of attack. Moreover, in the absence of acceleration constraints, the optimized trajectories involved a rapid transition to forward flight, followed by climbing at a relatively constant speed, and then accelerating to the desired cruising speed. When considering passenger comfort and incorporating acceleration constraints, the transition, climb, and acceleration phases exhibited more gradual and less distinct behavior, as expected. However, completed work on transition optimization as well as conventional engineering design optimizations rely on iterative simulation model evaluations, which prohibit real-time decision making.

Thus, surrogate models have emerged as an effective alternative for fast interactive decision making in lieu of time-consuming simulation models [15,16]. Surrogate-based design optimization is a methodology that leverages surrogate models to approximate the behavior of design candidates. To elaborate, surrogate models are mathematical models that are trained using a limited number of data points obtained from computationally expensive simulations or physical experiments. These surrogate models are then used in optimization algorithms to efficiently make predictions, explore design space, and find optimal solutions. Surrogate-based design optimization is advantageous over conventional simulation-based optimization for computational efficiency [17,18], the effective exploration of high-dimensional design spaces [19,20], multi-objective optimization capabilities [21,22], sensitivity analysis [23], and uncertainty analysis [24,25,26].

The Gaussian process (GP) [27] is one of the most commonly used surrogate models due to its capability and flexibility for modeling and predicting unknown functions based on observed data [28]. The concept of a multi-output GP (MOGP) has emerged to further handle problems with multiple outputs [29,30]. By extending the concept of a GP, an MOGP is a surrogate modeling technique that allows for the joint modeling of multiple correlated outputs and the simultaneous capture of their dependencies. In particular, by jointly modeling the outputs by estimating the covariance matrix (also known as the correlation function), an MOGP leverages the shared information among the outputs, which reduces the overall modeling complexity. For more mathematical and surrogate modeling details, please refer to Section 2.3. Meanwhile, deep learning surrogates are going through revolutionary developments and pushing forward cutting-edge research in the design optimization community. Thelen et al. [31] computed design variable derivatives by analytically linking objective definition, mesh, and geometry, and they demonstrated the derivatives through the multi-fidelity Broyden–Fletcher–Goldfarb–Shannon algorithm for high-dimensional aerodynamic and aeroelastic design optimizations. Tao and Sun [32] developed a multi-fidelity surrogate-based optimization framework based on deep belief networks. The results for airfoil and wing designs under uncertainty indicated that the multi-fidelity surrogate model performed well by significantly improving optimization efficiency. Renganathan et al. [33] performed aerodynamic design optimization by fusing deep neural networks (DNN) and GP as one surrogate to incorporate the predictive power of DNN and the confidence interval of GP. Thus, the DNN-GP surrogate enabled relatively high-dimensional Bayesian optimization on aerodynamic design and outperformed adjoint-based optimization in terms of efficiency.

Even with these developments, surrogate modeling still cannot completely address the “curse of dimensionality”, which necessitates dimensionality reduction techniques for further improvements. O’Leary-Roseberry et al. [34] constructed adaptive residual networks on reduced-dimensional space exploited by principal component analysis to directly predict optimal designs based on design requirements. They successfully demonstrated outstanding performance over full-space feed-forward networks on aerodynamic wing design cases. Meanwhile, a prosperous family of machine learning models, generative adversarial networks (GAN), was developed as generative models [35] with the feature of automatic design space reduction. The competition between a generator and a discriminator allows them to generate new data that follow the same patterns as the training data. Therefore, when fed with realistic airfoil or wing design shapes, a GAN model enables implicit design space dimensionality reduction by filtering out unrealistic shapes and generating only realistic designs [36]. Du et al. [37] developed the B-spline-based generative adversarial networks (BSplineGAN) for intelligent airfoil parameterization with the UIUC airfoil database as training data. BSplineGAN automatically reduced the design space while maintaining sufficient shape variability, which was verified by fitting optimizations to arbitrary UIUC airfoils. They constructed DNN surrogates on the reduced space exploited by BSplineGAN and verified predictive performance on aerodynamic design cases in a fast, interactive manner.

Derived from the original GAN, GAN variants also handle regression tasks. Du and Martins [38] introduced a novel multi-fidelity surrogate modeling architecture, super-resolution GAN (SRGAN), for predicting airfoil pressure distributions based on low-fidelity counterparts. Specifically, the SRGAN generator generated super-resolution pressure distributions, while the discriminator aimed to distinguish between pressure distributions of the generated super-resolution shapes and the high-resolution data set. Thus, training the generator minimized the difference between the super-resolution shapes and corresponding high-fidelity (i.e., high-resolution) data and maximized the similarity of super-resolution with the high-fidelity data. The results showed that the SRGAN outperformed low-fidelity simulations and direct DNN by accurately capturing the locations and magnitudes of strong high-fidelity shocks. A conditional generative adversarial networks (cGAN) incorporates additional conditioning information into the generative process [39]. The goal of a cGAN is to generate samples that not only resemble the real data but also adhere to specified conditions. The conditioning information serves as a guide for the generator to generate samples that align with the desired conditions. Aggarwal et al. [40] conducted experiments demonstrating the effectiveness of a cGAN model for regression problems. The cGAN was successfully demonstrated on a real-world ailerons data set for an F-16 airplane. The cGAN managed to predict the control input on the aircraft’s ailerons, which was described by 40 continuous inputs. Ye et al. [41] proposed a novel GAN-based regression model (regGAN), adopting a combined loss function on a mean squared error (MSE) and a binary cross-entropy (BC) loss, which showed outstanding predictive performance on frying oil deterioration when provided with time and temperature as input parameters.

Prior research has demonstrated the potential of machine learning surrogates, including GAN variants. The SRGAN and cGAN exhibited predictive potential, while the regGAN realized predictions from physical input space to output space. We demonstrate the regGAN performance in predicting optimal eVTOL takeoff trajectories directly based on design requirements and compare results with the MOGP and cGAN surrogates. We summarize the contribution of this paper as follows. First, we develop and introduce the regGAN surrogate into the takeoff trajectory design area for the first time. Second, we realize the surrogate-based optimal takeoff trajectory design for eVTOL aircraft to fill the lack of research study in this field. Third, we introduce the inverse mapping concept (from design requirements, including design constraints and flight condition parameters, directly into optimal designs) to eVTOL trajectory design for the first time.

We organize the rest of the paper as follows. Section 2 introduces the optimization framework and simulation models used in this work, followed by MOGP, cGAN, and regGAN setups for surrogate modeling. We demonstrate the regGAN as well as other surrogates on eVTOL optimal takeoff trajectory predictions in Section 3. We end this paper with conclusions in Section 4.

2. Methodology

In this section, we introduce the open-source optimization toolbox Dymos (https://github.com/OpenMDAO/RevHack2020/tree/master/problems/evtol_trajectory/evtol_dymos, accessed on 1 August 2022) [42] within OpenMDAO (https://github.com/OpenMDAO/OpenMDAO, accessed on 1 August 2022) [43] and simulation models. Additionally, we detail an MOGP toolkit (MOGPTK (https://github.com/GAMES-UChile/mogptk, accessed on 1 August 2022)) [44] followed by DNN and GAN series models (namely, GAN, cGAN, and regGAN) within TensorFlow [45] for surrogate modeling.

2.1. Dymos Framework

OpenMDAO [43], an open-source framework for multidisciplinary design, analysis, and optimization (MDAO), empowers engineers and scientists to efficiently analyze and optimize complex engineering systems. It is primarily designed for gradient-based optimization, which is based on efficient coupled derivative computation and problem sparsity exploitation [43]. Its component-based architecture enables the integration of models from diverse disciplines, facilitating seamless collaboration and reusability. OpenMDAO excels in interdisciplinary coupling, managing data flow, solving systems of equations, and offering efficient automatic differentiation for optimization and sensitivity analysis. Additionally, OpenMDAO supports parallel computing, leveraging multiple cores and distributed computing clusters for efficient analyses.

Dymos [42] is an advanced software toolkit built upon the OpenMDAO framework by the National Aeronautics and Space Administration for optimizing complex aerospace systems. Dymos expands OpenMDAO capabilities by offering specialized tools and algorithms for dynamic optimization. With a primary focus on dynamic systems, Dymos transcends solving problems that involve the evolution of systems over time. This is valuable for trajectory optimizations, where various constraints and objectives need to be considered. By employing numerical methods, such as direct transcription and indirect shooting methods, Dymos can discretize and optimize the trajectory based on specified nodes or iteratively refine initial guesses, respectively. The software also provides an intuitive interface and visualization tools for easy problem formulation and manipulation while supporting integration with external software and models. Dymos serves as a powerful and versatile toolkit for dynamic optimization and multidisciplinary analysis, enabling enhanced design and optimization of complex aerospace systems.

In this work, we utilize one of the implicit collocation techniques demonstrated in Dymos: higher-order Gauss–Lobatto collocation (LGL). Implicit collocation techniques are especially suitable for gradient-based optimization including analytic derivatives [42]. The equation of motion propagates explicitly through time-marching and requires a fixed number of analysis points in time with enough density to capture the system dynamics accurately. Implicit collocation schemes are able to attain this with fewer analysis points in a grid fixed throughout an optimization. In Dymos, implicit collocation methods further disintegrate each phase into one or more segments. For each segment, to ensure the value in time can be modeled as a polynomial, each state or control variable is assumed to be $C^1$ continuous. Additionally, within each segment, state and control values are allocated at specific points or nodes, indicated by the state discretization nodes and control discretization nodes, respectively. Inside each phase, there are several variables described in the following, such as design variables (initial time $t_0$, duration $t_p$), state variables at the state discretization nodes (${\overline{x}}_d$), dynamic control variables at the control discretization nodes (${\overline{u}}_d$), and the design parameters ($\overline{d}$). Herman and Conway proposed the high-order LGL method, which is the primary method used in OTIS [46] and is a general extension of the Hermite-Simpson collocation to higher orders [47]. In a phase, the segments are discretized at the LGL nodes so that, for each segment, the total number of nodes is an odd value. The state discretization nodes are the nodes with an even index, where the index starts at zero for the first node. Dynamic controls are provided at designated control discretization nodes. In LGL collocation, all nodes within the phase are involved in the control discretization subset. The evaluation of an LGL phase in Dymos progresses with the following steps: Initially, the values of the design variables ($t_0$, $t_p$, ${\overline{x}}_d$, ${\overline{u}}_d$, and $\overline{d}$) are given by the optimizer and user. Secondly, the OpenMDAO system offering the state dynamics is assessed at the state discretization nodes, contributing the state time derivatives at the state discretization nodes:

$${\dot{\overline{x}}}_d=f_{\mathrm{ode}}(t_d,{\overline{x}}_d,{\overline{u}}_d,\overline{d}),$$

(1)

where $f_{\mathrm{ode}}$ is the governing ordinary differential equation (ODE) for simulation models.

Then, using Hermite interpolation [48], the states and state rates are interpolated to the collocation nodes:

$${\overline{x}}_c=\left[A_i\right]x_d+\frac{t_{\mathrm{seg}}}{2}\left[B_i\right]{\dot{\overline{x}}}_d,$$

(2)

$${\overline{x}}_c^{\prime }=\frac{2}{\mathrm{seg}}\left[A_d\right]x_d+\left[B_d\right]{\dot{\overline{x}}}_d,$$

(3)

where $t_{\mathrm{seg}}$ is the duration of the segment where each nodes originates. At the collocation nodes, the state system provided by the system at this point is evaluated and is given by

$${\dot{\overline{x}}}_c=f_{\mathrm{ode}}(t_c,{\widehat{x}}_c,{\widehat{u}}_c,d).$$

(4)

Consequently, there are two independent state rate values existing: one is the actual value derived from the ODE, and the other is the approximate value calculated from the slope of the interpolating polynomials. For each state in the phase, the “collocation defects, $\Delta$”, normalized by segment time space, consists of the difference between two state rates. Moreover, these values are constrained at zero for a nonlinear programming problem and can be expressed as

$$\overline{\Delta }={\overline{x}}_c^{\prime }-\frac{{\overline{t}}_{\mathrm{seg}}}{2}{\dot{\overline{x}}}_{\mathrm{col}}.$$

(5)

2.2. Simulation Models

We consider the verified simulation models covering aerodynamics, propulsion, propeller–wing interaction, and dynamics developed by the Dymos team [42].

2.2.1. Aerodynamics

The aerodynamics model assumes that the forward and rear wings are equivalent and have the same reference area for simplicity, and there is no flow interaction between the two wings. The rotations of the two wings are assumed to be the same, so that the angles of attack are equivalent as well as the lift and drag generated by the two wings. During the transition from vertical to horizontal flight, separated-flow conditions are considered. A model developed by Tangler and Ostowari [49] for wing aerodynamics is implemented to predict the lift and drag beyond the linear-lift region. The poststall lift coefficient is

$$C_L=A_1\sin 2\alpha +A_2\frac{{\cos }^2\alpha }{\sin \alpha },$$

(6)

where

$$A_1=\frac{C_1}{2},$$

(7)

$$A_2=(C_{L_s}-C_1\sin {\alpha }_s\cos {\alpha }_s)\frac{\sin {\alpha }_s}{{\cos }^2{\alpha }_s},$$

(8)

and

$$C_1=1.1+0.018AR,$$

(9)

where $\alpha$ is the wing angle of attack (in radians), ${\alpha }_s$ is the angle of attack at stall (in radians), $C_{L_s}$ is the lift coefficient at stall, and $AR$ is the wing aspect ratio.

The drag coefficient at a wing angle of attack between 27.5 and 90 degrees is given by

$$C_D=B_1\sin \alpha +B_2\cos \alpha ,$$

(10)

where

$$B_1=C_{D_{\max }},$$

(11)

$$B_2=\frac{C_{D_s}-C_{D_{\max }}\sin {\alpha }_s}{\cos {\alpha }_s},$$

(12)

and

$$C_{D_{\max }}=\frac{1.0+0.065AR}{0.9+t/c},$$

(13)

where $C_{D_s}$ is the drag coefficient at stall, and $t/c$ is the airfoil thickness-to-chord ratio. For the poststall drag coefficient at a wing angle of attack below 27.5 degrees, the equation is

$$C_D=0.008+1.107{\alpha }^2+1.792{\alpha }^4,$$

(14)

where $\alpha$ is in radians.

The well-known finite wing corrections from lifting-line theory for unswept wings in incompressible flow are utilized to modify the lift and drag of the airfoil prior to stalling,

$${\alpha }_{\mathrm{wing}}=\frac{{\alpha }_{\mathrm{airfoil}}}{1+\left({\alpha }_{\mathrm{airfoil}}/(\pi ·AR·e)\right)},$$

(15)

where ${\alpha }_{\mathrm{wing}}$ is the finite-wing lift–curve slope, ${\alpha }_{\mathrm{airfoil}}$ is the airfoil lift–curve slope, and e is the span efficiency factor. The prestall lift curve is assumed to be linear, and the wing stall angle of attack is 15 deg.

In order to obtain the total drag of the wing before stall, induced drag has been added using the well-known formula based on the lifting-line theory for prestall parasite drag coefficients,

$$C_{D_i}=\frac{C_L^2}{\pi ·AR·e},$$

(16)

where $C_{D_i}$ is the induced drag coefficient, $C_L$ is the wing’s lift coefficient, and $AR=8$ for each wing of the configuration.

Additional drag on an assumed drag area for the fixed landing gear and the fuselage is assumed to be independent of the freestream angle of attack. The induced drag for the configuration can be expressed as

$$D_{\mathrm{induced}}=2\left(1.4\frac{L_{\mathrm{wing}}^2}{\pi q{b}^2·0.95}\right)=2\left(\frac{L_{\mathrm{wing}}^2}{\pi q{b}^2·0.68}\right),$$

(17)

where $L_{\mathrm{wing}}$ is lift per wing, q is freestream dynamic pressure, b is wing span.

2.2.2. Propulsion

Momentum theory [50] has been utilized to compute thrust from the propellers as a function of power.

$$P_{\mathrm{disk}}=T{V}_{\infty \perp }+\kappa T\left(-\frac{V_{\infty \perp }}{2}+\sqrt{\frac{V_{\infty \perp }^2}{4}+\frac{T}{2\rho A_{\mathrm{disk}}}}\right),$$

(18)

where $P_{\mathrm{disk}}$ is the power supplied to the propeller disk excluding profile power, T is the thrust, $V_{\infty \perp }$ is the freestream velocity component normal to the propeller disk, $\rho$ is the air density, $A_{\mathrm{disk}}$ is the disk area of the propeller, and $\kappa$ is the correction factor utilized to incorporate induced power losses associated with non-uniform flow, tip effects, and other simplifications made in momentum theory ($\kappa =1$ for ideal power). Power is used as a design variable in the optimization problems, and the Newton–Raphson method [51] has been chosen to solve the nonlinear equation for thrust with power as an input.

Applying blade-element theory to a rotor operating in nonaxial forward flight [52,53] estimates the profile power coefficient $C_{P_p}$ and profile power $P_p$

$$C_{P_p}=\frac{\sigma C_{d_{0p}}}{8}(1+4.6{\mu }^2),$$

(19)

and

$$C_{P_p}=\frac{P_p}{\rho A_{\mathrm{disk}}{R}^3{\Omega }^3},$$

(20)

where R is the radius of the propeller, $\Omega$ is the angular speed, $\sigma$ is the solidity, $C_{d_{0p}}$ is a representative constant profile drag coefficient, and $\mu$ is defined as

$$\mu =\frac{V_{\infty \Vert }}{\Omega R},$$

(21)

where $V_{\infty \Vert }$ is the freestream velocity component parallel to the disk. For our cases, we assume that $\Omega =181$ rad/s for $R=0.75$ m, $\sigma =0.13$, and $C_{d_{0p}}=0.012$.

With electrical power as an input, we use a factor $k_{\mathrm{elec}}$ to account for electrical and mechanical losses related to batteries, electrical systems, and motors. $P_{\mathrm{disk}}$ can be expressed as

$$P_{\mathrm{disk}}=k_{\mathrm{elec}}{P}_{\mathrm{electrical}}-P_p,$$

(22)

where $P_{\mathrm{electrical}}$ is the power from the batteries. We set the limit on the maximum electrical power available to 311 kW. The loss factor $k_{\mathrm{elec}}$ ranges from 0.7 to 0.9 and is also one of our design requirements considered as input parameters for surrogate modeling.

Furthermore, when the freestream flow is not completely normal to the propeller disks, the normal force N is calculated as

$$N=\frac{4.25{\sigma }_{e}sin(\beta +8^{\circ })f{q}_{\perp }{A}_{\mathrm{disk}}}{1+2{\sigma }_e}tan{\alpha }_{in},$$

(23)

where $q_{\perp }$ is the dynamic pressure based on the freestream velocity component normal to the propeller disk, $\beta$ is the blade pitch angle at $0.75R$ and is assumed to change linearly from 10 deg at a flight speed of 0 m/s to 35 deg at a cruise speed of 67 m/s, $A_{\mathrm{disk}}$ is the propeller area, and ${\alpha }_{in}$ is the incidence angle. The remaining terms are calculated in the following manner. The effective solidity ${\sigma }_e$ is

$${\sigma }_e=\frac{2B{c}_{\mathrm{b}}}{3\pi R},$$

(24)

where B is the number of blades per propeller, $c_{\mathrm{b}}$ is the average chord length of the blades, and R is the propeller radius. The thruster factor f is

$$f=1+\frac{\sqrt{1+T_c}-1}{2}+\frac{T_c}{4(2+T_c)},$$

(25)

where $T_c$ is a thrust coefficient, defined as

$$T_c=\frac{T}{q_{\perp }{A}_{\mathrm{disk}}},$$

(26)

and T is the thrust.

2.2.3. Propeller–Wing Interaction

The momentum theory [50] is implemented for modeling the interaction between a wing and the flow induced by propellers. The induced speed at the disk $v_i$ is given by

$$v_i=\left(-\frac{V_{\infty \perp }}{2}+\sqrt{\frac{V_{\infty \perp }^2}{4}+\frac{T}{2\rho A_{\mathrm{disk}}}}\right).$$

(27)

To explain the propeller–wing interaction in this work, the chordwise component of the freestream velocity for the entire wing by a range of factors $k_{\mathrm{in}}$ is varied between 0.3 and 1, multiplied by the induced speed at the disk $v_i$. The effective $k_{\mathrm{in}}$ for the wing is anticipated to be close to 1 but slightly less than it, especially when the aircraft is at low speed. Moreover, cases will become infeasible if $k_{\mathrm{in}}$ is too low; therefore, we set the range $k_{\mathrm{in}}\in [0.3,1.0]$ in this paper to make sure all cases are feasible under the design requirements.

2.2.4. Dynamics

A two degrees of freedom (DOF) representation and the forward Euler method [54] are used for simulating the trajectory of the aircraft. The horizontal and vertical components of the aircraft velocity are solved as functions of time, considering the control variables, which include the wing-tilt angle and electrical power. The horizontal component of velocity at each time step is calculated as follows:

$$v_{x_{i+1}}=v_{x_i}+\frac{Tsin\theta -D_{\mathrm{fuse}}sin(\theta +{\alpha }_{\infty })-D_{\mathrm{wings}}sin(\theta +{\alpha }_{\mathrm{EFS}})-L_{\mathrm{wings}}cos(\theta +{\alpha }_{\mathrm{EFS}})-Ncos\theta }{m}\Delta t,$$

(28)

where i is index of a time step, $\Delta t$ is a time step length, $\theta$ is a wing angle relative to the vertical, ${\alpha }_{\infty }$ is a freestream angle of attack, ${\alpha }_{\mathrm{EFS}}$ is the effective freestream angle of attack experienced by the wings due to the propeller influence, m is the mass of the aircraft, T is the total thrust, $D_{\mathrm{fuse}}$ is the drag of the fuselage, $D_{\mathrm{wings}}$ is the total drag of the two wings, $L_{\mathrm{wings}}$ is the total lift of the two wings, and N is the total normal force generated by the propellers. Figure 1 displays the forces and the angles on the aircraft. Likewise, the vertical component of velocity at each time step is calculated as follows:

$$v_{y_{i+1}}=v_{y_i}+\frac{Tcos\theta -D_{\mathrm{fuse}}cos(\theta +{\alpha }_{\infty })-D_{\mathrm{wings}}cos(\theta +{\alpha }_{\mathrm{EFS}})+L_{\mathrm{wings}}sin(\theta +{\alpha }_{\mathrm{EFS}})+Nsin\theta -mg}{m}\Delta t,$$

(29)

where g is the gravitational acceleration. To determine the time step length, $\Delta t$, the total takeoff time, which is a design variable in optimization problem formulation, is divided by 500. A convergence study concluded the optimized takeoff flight time would be typically below 50 s, resulting in adequately small time steps of less than 0.1 s [14].

2.3. Multi-Output Gaussian Processes

A GP is a Bayesian nonparametric model that is designed by parameterizing a covariance kernel, meaning that constructing expressive kernels allows for an effective representation of complex systems [44]. The GP assumes the function $f\left(\mathbf{x}\right)$ is distributed as [55]

$$f\left(\mathbf{x}\right)\sim \mathcal{GP}\left(m\left(\mathbf{x}\right),k\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)\right),$$

(30)

where a Gaussian process is a distribution over functions and is characterized by a mean function $m\left(\mathbf{x}\right)$ and a covariance function $k\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)$ between input parameters $\mathbf{x}$ and ${\mathbf{x}}^{\prime }$. The mean function represents the expected function value at input $\mathbf{x}$

$$m\left(\mathbf{x}\right)=\mathbb{E}\left[f\right(\mathbf{x}\left)\right].$$

(31)

In other words, the mean function represents the average function value at input $\mathbf{x}$ across all functions in the distribution. To simplify posterior computations and rely only on the covariance function for inference, the prior mean function is frequently assumed to be a constant. In this work, we apply the constant mean function ($m\left(\mathbf{x}\right)=b_1$) and the linear mean function ($m\left(\mathbf{x}\right)=a\mathbf{x}+b_2$), where a, $b_1$, and $b_2$ are arbitrary constant coefficients. The covariance function $k\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)$ captures the relationship between the function values at different input points $\mathbf{x}$ and ${\mathbf{x}}^{\prime }$:

$$k\left(\mathbf{x},{\mathbf{x}}^{\prime }\right)=\mathbb{E}\left[(f\left(\mathbf{x}\right)-m\left(\mathbf{x}\right))\left(f\left({\mathbf{x}}^{\prime }\right)-m\left({\mathbf{x}}^{\prime }\right)\right)\right].$$

(32)

The function k is commonly called the kernel of the Gaussian process. The kernel defines the shape and characteristics of the GP model. We implement the Matérn kernel and the square exponential kernel function in our study, and they are described as follows.

• Matérn kernel
  The Matérn covariance between two points with the distance $\tau =|\mathbf{x}-{\mathbf{x}}^{\prime }|$ is

  $$k_{\nu }\left(\tau \right)={\sigma }^2\frac{2^{1-\nu }}{\Gamma \left(\nu \right)}\left(\sqrt{2\nu }\frac{\tau }{\rho }\right)K_{\nu }\left(\sqrt{2\nu }\frac{\tau }{\rho }\right),$$

  (33)

  where $\Gamma$ is the gamma function, $K_{\nu }$ is the modified Bessel function of the second kind, and $\rho$ and $\nu$ are non-negative covariance parameters. The Matérn kernel is stationary since the covariance only depends on distances between points.
• Squared exponential kernel

  $$k_{SE}=s^{2}exp\left(-\frac{1}{2}\left(\frac{|\mathbf{x}-{\mathbf{x}}^{\prime }|}{\mathbf{\lambda }}\right)\right),$$

  (34)

  where s is the scale factor. In the case where we use a single length scale parameter, the squared exponential (SE) kernel is an example of a radial basis function.

The GP has been extended to multiple series (or channels), which are referred to as MOGP. The key difference lies in their implementation and capabilities. The GP is designed for single-output regression, predicting a single target variable given inputs, while an MOGP extends the GP to handle multiple outputs simultaneously. An MOGP is useful when dealing with correlated outputs, enabling the modeling of dependencies between different target variables. The complexity of a GP relates to the number of data points, while an MOGP introduces additional complexity due to modeling correlations between multiple outputs, but the MOGP remains efficient and more capable compared with training individual GPs for each output. Thus, a key feature of the MOGP is to harness applicable information across outputs to provide more accurate predictions than separately modeling correlated outputs [28].

In this study, the multi-output spectral mixture [56] kernel is used as the multivariate kernel for the MOGP model, and can be expressed as

$$k_{ij}\left(\tau \right)=\sum _{q=1}^Q{\alpha }_{ij}^{\left(q\right)}exp\left(-\frac{1}{2}{\left(\tau +{\theta }_{ij}^{\left(q\right)}\right)}^{\top }{\Sigma }_{ij}^{\left(q\right)}\left(\tau +{\theta }_{ij}^{\left(q\right)}\right)\right)cos\left({\left(\tau +{\theta }_{ij}^{\left(q\right)}\right)}^{\top }{\mu }_{ij}^{\left(q\right)}+{\varphi }_{ij}^{\left(q\right)}\right),$$

(35)

where ${\alpha }_{ij}^{\left(q\right)}=w_{ij}^{\left(q\right)}{\left(2\pi \right)}^{\frac{n}{2}}{\left|{\Sigma }_{ij}^{\left(q\right)}\right|}^{1/2}$ and the superindex ${(·)}^{\left(q\right)}$ denotes the parameter of the $q^{\mathrm{th}}$ component of the spectral mixture. In addition, the cross-spectral density between channels i and j is modeled as a complex-valued SE function with the following parameters:

• Covariance: ${\Sigma }_{ij}=2{\Sigma }_i{\left({\Sigma }_i+{\Sigma }_j\right)}^{-1}{\Sigma }_j$;
• Mean: ${\mu }_{ij}={\left({\Sigma }_i+{\Sigma }_j\right)}^{-1}\left({\Sigma }_i{\mu }_j+{\Sigma }_j{\mu }_i\right)$;
• Magnitude: $w_{ij}=w_i{w}_{j}exp\left(-\frac{1}{4}{\left({\mu }_i-{\mu }_j\right)}^{\top }{\left({\Sigma }_i+{\Sigma }_j\right)}^{-1}\left({\mu }_i-{\mu }_j\right)\right)$;
• Delay: ${\theta }_{ij}={\theta }_i-{\theta }_j$;
• Phase: ${\varphi }_{ij}={\varphi }_i-{\varphi }_j$;

where the so-constructed magnitudes $w_{ij}$ ensure positive definiteness.

We use the MOGPTK [44] to establish MOGP models. MOGPTK is an open-source Python package that provides a natural way to train and use the models. MOGPTK is built upon GPflow [57], an extensive GP framework with a wide variety of implemented kernels, likelihoods, and training strategies. The main components of MOGPTK include MOGP modeling, data handling, parameter initialization, and parameter interpretation [44].

2.4. Deep Neural Networks

The DNN makes the core of deep learning. In this section, we introduce the DNN model basic setup and the Adam optimizer used for model training.

2.4.1. DNN Model Setup

A DNN [58] consists of multiple layers of interconnected nodes, known as neurons, to process input data and make predictions. A DNN starts with an input layer (reading input data), followed by one or more hidden layers, each of which consists of multiple neurons. The final hidden layer connects to the output layer, which produces the predictions. Moreover, neurons in the same layers share a weight matrix ($\mathbf{W}$) and a bias vector ($\mathbf{b}$) to be tuned. The equation within each neuron can be mathematically expressed as

$$\mathbf{o}=\sigma (\mathbf{W}·\mathbf{h}+\mathbf{b}),$$

(36)

where $\mathbf{h}$ is a data vector from previous layers, $\sigma (·)$ represents an activation function used in the current layer, and $\mathbf{o}$ is an output vector of the current layer. During training, $\mathbf{W}$ and $\mathbf{b}$ are adjusted to optimize the model’s performance. We utilize two commonly used activation functions in this work to introduce nonlinearity into the DNN. The rectified linear unit (relu) [59] activation function follows the equation ${\sigma }_{\mathtt{relu}}\left(\mathbf{h}\right)=max(\mathbf{0},\mathbf{h})$, while the sigmoid [60] activation function has an expression of ${\sigma }_{\mathtt{sigmoid}}\left(\mathbf{h}\right)=1/\left(1+e^{-\mathbf{h}}\right)$.

In this work, we harness two types of loss functions for training DNN surrogates to minimize the difference between predicted responses and corresponding actual values. The MSE [61] loss function is formulated as

$${\mathcal{L}}_{\mathrm{MSE}}=\frac{1}{N}\sum _{i=1}^N{\left({\left|{\mathbf{v}}_i-{\widehat{\mathbf{v}}}_i\right|}_2\right)}^2,$$

(37)

where N is the number of training data points in the data set, ${|·|}_2$ is $L_2$ norm, ${\mathbf{v}}_i$ and ${\widehat{\mathbf{v}}}_i$ are actual values and predicted values for the ith training sample. Binary crossentropy (BC) [62] is formulated as

$${\mathcal{L}}_{\mathrm{BC}}=-\frac{1}{N}\sum _{i=1}^N\left({\mathbf{v}}_i·log{\widehat{\mathbf{v}}}_i+\left(\mathbf{1}-{\mathbf{v}}_i\right)·log\left(\mathbf{1}-{\widehat{\mathbf{v}}}_i\right)\right),$$

(38)

where ${\mathbf{v}}_i$ is the ith actual binary label ($\mathbf{0}$ or $\mathbf{1}$), and ${\widehat{\mathbf{v}}}_i$ is the predicted probability (between $\mathbf{0}$ and $\mathbf{1}$). The loss function penalizes the model more severely when it predicts the opposite class with high confidence. In this work, the DNN unknown parameters are tuned by minimizing loss functions through the gradient-based Adam optimizer (see the following section) enabled by backpropagation within Tensorflow [45].

2.4.2. Adam Optimizer

The Adam optimizer combines the principles of gradient descent with momentum [63] and root mean square propagation (RMSP) [64] algorithms [65]. The gradient descent algorithm has been enhanced by the moment algorithm [63], with the updates as

$$w_{t+1}=w_t-\alpha m_t,$$

(39)

where

$$m_t=\beta m_{t-1}+(1-\beta )\left[\frac{\partial L}{\partial w_t}\right],$$

(40)

where the subscripts $t-1$, t, and $t+1$ are the indices of previous, current, and next optimization steps, respectively. Initially, $m_0=0$, w is a weight to be determined, $\alpha$ is learning rate, $\frac{\partial L}{\partial w_t}$ is the derivative of loss function with respect to weight at current optimization step, and $\beta$ is a constant moving average parameter.

RMSP [64] is an adaptive learning algorithm with a mathematical formulation as

$$w_{t+1}=w_t-\frac{\alpha }{{(v_t+ϵ)}^{1/2}}{\left[\frac{\partial L}{\partial w_t}\right]}^2,$$

(41)

where

$$v_t=\beta v_{t-1}+(1-\beta )\left[\frac{\partial L}{\partial w_t}\right].$$

(42)

A small positive constant ($ϵ={10}^{-7}$) is utilized. The remaining variables remain the same as in Equation (40).

The Adam optimizer integrates the moment and RMSP algorithms as follows.

$$w_{t+1}=w_t-{\widehat{m}}_t\left(\frac{\alpha }{\sqrt{{\widehat{v}}_t}+ϵ}\right),$$

(43)

where

$${\widehat{m}}_t=\frac{m_t}{1-{\beta }_1^t},{\widehat{v}}_t=\frac{v_t}{1-{\beta }_2^t},$$

(44)

where $m_t$ and $v_t$ follow the updating process as described above.

2.5. Generative Adversarial Network Models

A GAN is a novel DNN architecture that is leading in state-of-the-art generative modeling strategies. This section introduces GAN setups along with their variations, namely, cGAN and regGAN.

2.5.1. Generative Adversarial Networks

A GAN is in nature a type of generative model that consists of two neural networks: a generator and a discriminator (Figure 2) [35]. The generator initially produces random outputs based on random input variables, but as it receives feedback from the discriminator, it learns to generate more realistic samples through training. The discriminator, on the other hand, is trained on a data set of existing samples and the generated samples created by the generator. The discriminator is to distinguish between existing and generated data and provide feedback to the generator by assigning probabilities to the existing and generated samples, indicating how likely they are to be real. During the training process, the generator and the discriminator attempt to compete with each other in an adversarial manner. This competition drives the improvement of both models through the model training process. Thus, as the training proceeds, the generator is more capable of producing similar samples as existing data that the discriminator cannot differentiate.

Mathematically, training a GAN model can be expressed as a minimax problem [35] on the loss function as follows.

$$\underset{G}{min}\underset{D}{max}{\mathcal{L}}_{\mathrm{GAN}}(D,G)={\mathbb{E}}_{\mathbf{x}\sim P_{\mathrm{data}}}\left[logD\left(\mathbf{x}\right)\right]+{\mathbb{E}}_{\mathbf{z}\sim P_{\mathbf{z}}}\left[log(1-D(G\left(\mathbf{z}\right)\left)\right)\right],$$

(45)

where $\mathbf{x}$ is sampled from the existing data distribution $P_{\mathrm{data}}$, $\mathbf{z}$ is sampled from the noise variable distribution $P_{\mathbf{z}}$, and G and D are the generator and discriminator. Thus, a trained GAN model generates reasonable designs with ample shape variability within prior noise variable distributions.

2.5.2. Conditional Generative Adversarial Networks

A cGAN is an extension of the original GAN that incorporates additional conditioning information into the generative process [39]. For the generator and the discriminator in cGANs, both networks receive not only random noises as input but also additional conditioning information $\mathbf{y}$. This information offers extra guidance to the generator to produce samples that align with specified conditions. Conditioning information $\mathbf{y}$ can be in several forms, such as text descriptions and class labels, depending on the utilization. The conditioning is achieved by providing $\mathbf{y}$ as an additional input layer to both the discriminator and generator. Within the generator, the prior input noise $P_{\mathbf{z}}\left(\mathbf{z}\right)$ and $\mathbf{y}$ are combined in a joint hidden representation, benefiting from the adversarial training framework’s flexibility in composing this hidden representation. In addition, the discriminator takes as inputs $\mathbf{x}$ and $\mathbf{y}$ pairs as well as shapes generated by the generator to complete the adversarial competition with the generator. Similarly, as a GAN loss function (Equation (45)), the training on loss functions could also be expressed as a minimax problem,

$$\underset{G}{min}\underset{D}{max}{\mathcal{L}}_{\mathrm{cGAN}}(D,G)={\mathbb{E}}_{\mathbf{x}\sim P_{\mathrm{data}}\left(\mathbf{x}\right)}\left[logD\left(\mathbf{x}\right|\mathbf{y})\right]+{\mathbb{E}}_{\mathbf{z}\sim P_{\mathbf{z}}\left(\mathbf{z}\right)}\left[log(1-D(G\left(\mathbf{z}\right|\mathbf{y}\left)\right))\right].$$

(46)

In addition, cGANs can be adapted for regression tasks by encoding the regression labels (i.e., model observations) as additional conditioning information (Figure 3). The cGAN generator takes random noise variables as the original GAN, together with arbitrary model observations as an additional input group. The goal is to generate data samples that correspond to the model observations. The discriminator, which tries to distinguish between real and generated data samples, will now take both real data samples $\mathbf{x}$ and generated data samples $\mathbf{g}$ as input. During training, the cGAN generator learns to map the random noise variables and model observations to corresponding data samples, which helps in regression tasks by generating samples consistent with the given target values. This approach offers a way to tackle problems with complex relationships between inputs and targets, and incorporate the adversarial feature of GANs into regression processes. In this work, $\mathbf{x}$ represents the eVTOL optimal takeoff trajectory control points to be predicted, $\mathbf{y}$ includes design requirements, and we set $\mathbf{z}$ with only one element and assign a Uniform (0, 1) distribution to $\mathbf{z}$. During regression tasks, an average of predictions over a fixed set of $\mathbf{y}$ and 100 randomly sampled $\mathbf{z}$ is used as the predicted model response [39], which in this work, is the optimal takeoff trajectory ($\mathbf{x}$). The network architecture settings for the generator and discriminator within the cGAN follow the same architecture as the regGAN, which will be discussed in Section 2.5.3.

2.5.3. Regression Generative Adversarial Networks

To further investigate GAN models for regression tasks, we developed the regGAN for surrogate modeling. The regGAN (Figure 4) shares a similar structure as the cGAN (Figure 3), but the regGAN has only design requirements ($\mathbf{y}$) as the generator input. Moreover, training regGAN incorporates a combined loss function of MSE (Equation (37)) and BC (Equation (38)), such that the regGAN couples a predictive feature of surrogate models and a generative feature of GAN models. The combined loss function is applied in model training.

$$\begin{array}{cc}\hfill \underset{G}{min}\underset{D}{max}{\mathcal{L}}_{\mathrm{regGAN}}(D,G)& =min\left({\mathcal{L}}_{\mathrm{context}}+{\mathcal{L}}_{\mathrm{adversarial}}\right)\hfill \end{array}$$

(47)

$$\begin{array}{cc}\hfill & =w_{\mathrm{BC}}·\left(\underset{G}{min}\underset{D}{max}{\mathcal{L}}_{\mathrm{BC}}(D,G)\right)+w_{\mathrm{MSE}}·\left(\underset{G}{min}{\mathcal{L}}_{\mathrm{MSE}}\left(G\right)\right)\hfill \end{array}$$

(48)

$$\begin{array}{cc}\hfill & =w_{\mathrm{BC}}·\underset{G}{min}\underset{D}{max}\left({\mathbb{E}}_{\mathbf{x}\sim P_{\mathrm{data}}}\left[logD\left(\mathbf{x}\right)\right]+{\mathbb{E}}_{\mathbf{y}\sim P_{\mathbf{y}}}\left[log(1-D(G\left(\mathbf{y}\right)\left)\right)\right]\right)\hfill \end{array}$$

(49)

$$\begin{array}{cc}\hfill & +w_{\mathrm{MSE}}·\underset{G}{min}\frac{1}{N}\sum _{i=1}^N{\left({\left|{\mathbf{x}}_i-G{\left(\mathbf{y}\right)}_i\right|}_2\right)}^2.\hfill \end{array}$$

(50)

where $w_{\mathrm{BC}}$ and $w_{\mathrm{MSE}}$ are constant weights on BC and MSE loss functions; we conduct a parametric study in Section 3 for the best predictive performance.

The MSE loss (also known as context loss) guides the regGAN predictions to match model observations in a similar way as traditional surrogates, while the BC loss (also known as adversarial loss due to the competitive training on GAN models) on the discriminator engenders similar patterns between predicted results and observations. In this work, we leverage the regGAN for eVTOL takeoff trajectory inverse mapping, i.e., directly from design requirements to optimal takeoff trajectories. Thus, the regGAN generator reads the design requirements ($\mathbf{y}$) and predicts the corresponding optimal takeoff trajectories ($\mathbf{x}$). Meanwhile, the discriminator attempts to distinguish between predicted/generated trajectories and actual optimal takeoff trajectories. Thus, the nature of the eVTOL optimal takeoff trajectory prediction makes it necessary to develop and introduce the regGAN surrogate. First, the regGAN makes use of the predictive capability of the DNN via the generator. Second, the optimal takeoff trajectory profiles typically follow realistic patterns. For instance, optimal power and wing angle profiles follow ascending trends in general, meaning that the battery provides more and more power to balance the weight and gradually moves forward by turning the wing from vertical to horizontal positions. In this way, the BC loss based on the GAN architecture facilitates the training by automatically filtering out unrealistic trajectory profiles.

The regGAN generator architecture setting is as follows (

Table 1. Neural architectures and training setups of the generator and the discriminator within regGAN.

	Generator	Discriminator
Input layer	5 neurons, no activation	21 neurons, no activation
Hidden layer 1	100 neurons, relu activation	100 neurons, relu activation
Hidden layer 2	100 neurons, relu activation	100 neurons, relu activation
Output layer	21 neurons/1 neuron, sigmoid activation	1 neurons, sigmoid activation
Training algorithm	Adam optimizer	Adam optimizer
Training parameters	${\beta }_1=0.9$, ${\beta }_2=0.999$	${\beta }_1=0.9$, ${\beta }_2=0.999$
Learning rate	0.001	0.001
Batch size	20	20
Epochs	1000	1000

). The generator consists of one input layer, two hidden layers, and one output layer. Five neurons in the input layer correspond to five different design requirements. Both hidden layers have 100 neurons, followed by relu activation functions. The output layer has twenty-one neurons corresponding to the twenty-one time-sequence optimal design variables or one neuron for time prediction, followed by sigmoid activation functions. The discriminator includes one input layer, two hidden layers, and one output layer. The input layer has 21 neurons corresponding to the time-sequence optimal design variables. Both hidden layers have 100 neurons with relu activation functions. The output layer has only one neuron followed by sigmoid activation functions in order to distinguish between generated and true optimal design variables. The observations in training data sets by the generator are scaled within the range of [0, 1] via the MinMaxScaler within Scikit-learn. We finish training the regGAN using the Adam optimizer (Section 2.4.2) within Tensorflow with the momentum coefficients of ${\beta }_1=0.9$ and ${\beta }_2=0.999$. The learning rates for the generator and the discriminator networks are both 0.001. The batch size is 20 and the maximum number of training iterations is 1000. We use the same architectures and setups to train three separate regGAN surrogates to predict the power profile, wing angle profile, and total takeoff time. The same architectures and setups also apply to cGAN surrogates for fair comparisons.

2.6. Verification Metric

In this work, we implement the mean $L_1$ relative error (${\overline{ϵ}}_{L_1}$) to evaluate the surrogate predictive performance. The $L_1$ relative error, also known as the mean absolute percentage error, is a metric used to measure the relative $L_1$ norm difference between predictions compared with actual values.

$${\overline{ϵ}}_{L_1}=\frac{1}{N_{\mathrm{test}}}\sum _{i=1}^{N_{\mathrm{test}}}\frac{\left|{\mathbf{y}}_{\mathrm{pred},i}-{\mathbf{y}}_{\mathrm{true},i}\right|}{\left|{\mathbf{y}}_{\mathrm{true},i}\right|}\times 100\%,$$

(51)

where $N_{\mathrm{test}}$ is the number of testing samples, $y_{\mathrm{true},i}$ and $y_{\mathrm{pred},i}$ are the true observations and predicted values of the $i^{\mathrm{th}}$ set of design requirements. The relative accuracy is calculated as

$$AC{C}_{L_1}=1-{\overline{ϵ}}_{L_1}.$$

(52)

The $L_1$ relative accuracy measures the average relative match between predicted values and true observations, expressed in percentage. It gives an indication of to what extend surrogate predictions agree with true observations relative to true observations. The percentage exhibits deeper insights through the errors relative to model response magnitudes.

3. Results and Discussion

In this section, we formulate the optimization problem and vary the operating parameter bounds to consider various takeoff scenarios. We showcase and compare the optimal takeoff trajectory predictions by MOGP, cGAN, and regGAN surrogates with simulation-based optimal counterparts. The MOGP represents the traditional surrogates (such as polynomial chaos expansions) and shows promising capability in handling multiple outputs. The cGAN makes use of the GAN architecture for regression tasks that share similar principles as the regGAN. Note that the training and testing data sets are generated by simulation-based optimal designs, and the testing accuracy was calculated based on the testing data set. In addition, we do not consider buildings, although the application scenario is set up for UAM, since we mainly focus on the UAM scope of precise transportation/delivery through optimal energy-efficient takeoff designs.

Table 2. Trajectory optimization problem formulation.

	Function or Variable	Description	Quantity
minimize	E	Electrical energy consumed
w.r.t.	$\mathbf{P}$	Electric power using 21 cubic curve control points	21
	$\mathit{\theta }$	Wing angle to vertical using 21 cubic curve control points	21
	$t_{\mathrm{flight}}$	Takeoff time	1
		Total design variables	43
subject to	$y_{\mathrm{final}}\ge 305$ m	Final vertical displacement constraint	1
	$x_{\mathrm{final}}\le 1400$ m	Final horizontal displacement constraint	1
	$v_x=67$ m/s	Final horizontal speed constraint	1
	$y\ge 0$ m	Vertical displacement constraint	1
	$a\le a_{max}$	Acceleration constraint	1
	$\alpha \le {\alpha }_{lim}$	Positive stall-angle constraint	1
	$\alpha \ge -{\alpha }_{lim}$	Negative stall-angle constraint	1
		Total constraints	7
Conditions	$k_{\mathrm{in}}$	Propeller-induced velocity factor
	${\alpha }_{lim}$	Angle of attack constraint value
	$a_{max}$	Maximum acceleration magnitude
	$k_{\mathrm{elec}}$	Electrical and mechanical loses factor
	$S_{\mathrm{ref}}$	Wing size factor

formulates the trajectory optimization problem. The objective is to minimize the electrical power consumed to reach a minimum vertical displacement of 305 m and a minimum horizontal speed of 67 m/s. The design requirements (i.e., design constraints and flight condition parameters) include an angle of attack constraint ${\alpha }_{lim}\in [10,15]$ deg, maximum acceleration magnitude $a_{max}\in [0.2,0.4]g$ (g is gravitational acceleration), propeller-induced velocity factor $k_{\mathrm{in}}\in [0.3,1.0]$, electrical and mechanical loss factor $k_{\mathrm{elec}}\in [0.7,0.9]$, and wing size factor $S_{\mathrm{ref}}\in [0.9,1.0]$. Note that the $a_{\max }$ constraint is a concept to consider for passenger comfort for future real-world transportation, although eVTOL aircraft currently have not been widely applied for such tasks yet. The design variables are the time-sequence electrical power ($\mathbf{P}$) and wing angle to vertical ($\mathit{\theta }$), both of which have 21 cubic curve control points and a total takeoff time ($t_{\mathrm{flight}}$). As mentioned in Section 2, we use the open-source Dymos package within OpenMDAO. Figure 5 shows the optimal takeoff trajectory profiles, which verify the Dymos package in terms of the time history of design variables and takeoff conditions.

3.1. MOGP Surrogate Modeling

We use 1000 random Latin hypercube sampling (LHS) points as training data for MOGP surrogate modeling and 300 LHS testing samples to verify predictive performance.

Table 3. MOGP model conditions and mean ± standard deviation of testing accuracy on testing data set. SE stands for SE kernel function.

Model	Kernel Function	Mean Function	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, ${\mathit{t}}_{\mathbf{flight}}$ (%)	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, P (%)	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, $\mathit{\theta }$ (%)
M1	SE	Constant	99.6 ± 0.313	92.5 ± 2.93	94.7 ± 2.30
M2	SE	Linear	99.5 ± 0.372	92.5 ± 2.89	94.7 ± 2.32
M3	Matérn	Constant	96.7 ± 2.65	88.7 ± 4.85	94.6 ± 2.20
M4	Matérn	Linear	96.7 ± 2.61	91.4 ± 2.73	94.6 ± 2.20

shows the mean testing accuracy for each design variable group using different kernel and mean functions. The results show that the SE kernel function has better overall predictive performance than the Matérn kernel. The SE-based MOGP predicts the $t_{\mathrm{flight}}$, $\mathbf{P}$, and $\mathit{\theta }$ at mean testing accuracies of 99.5%, 92.5%, and 94.7%, respectively. The prediction difference between the constant and linear mean functions is negligible.

We chose an arbitrary case to further reveal the predictive performance of each MOGP surrogate.

Table 4. Design requirements for visualized case.

Design Requirements	Values
Propeller-induced velocity factor, $k_{\mathrm{in}}$	87.75 (%)
Angle of attack constraint, ${\alpha }_{\lim }$	±12.25833333 (deg)
Acceleration, $a_{\max }$	0.38766667 (g)
Electrical and mechanical loses factor, $k_{\mathrm{elec}}$	0.761
Wing size, $S_{\mathrm{ref}}$	0.9285

shows the design requirements for the visualization case.

Table 5. Testing accuracy of visualized case for MOGP models 1–4.

Model	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, ${\mathit{t}}_{\mathbf{flight}}$ (%)	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, P (%)	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, $\mathit{\theta }$ (%)
M1	99.9	94.4	91.0
M2	99.8	94.7	91.0
M3	95.9	91.5	90.7
M4	95.7	94.3	90.5

shows the testing accuracy on design variables for MOGP models 1–4. We compare the predicted optimal trajectory profiles by MOGP models 1–4 in Figure 6, where all models achieve similar predictive accuracies. The testing accuracies, shown in Table 5, are all over 90%; however, visualization exhibits obvious discrepancies between surrogate-based and simulation-based optimal trajectories. The results indicate that MOGP surrogates intend to use higher power early but not the maximum power. The early power usage results in higher thrust and greater acceleration in the early takeoff phase by MOGP surrogates than by the simulation-based optimal design. The acceleration by the surrogate-based optimal design in the early phase even violates the maximum acceleration constraint. So, we conclude that the MOGP cannot realize a sufficient accuracy level using 1000 training samples.

3.2. cGAN Surrogate Modeling

We implement one single-noise cGAN model and compare the predictive performance of this model against the Dymos simulation-based optimal trajectory as well as the MOGP predictions shown in Section 3.1. The model architecture and hyperparameter settings are as mentioned in Section 2.5.3. Here, we train the model with MSE and BC loss functions separately. To enable the regression feature with the cGAN, we take the mean value of 100 random noise inputs for each set of design requirements.

Table 6. cGAN models exhibit better predictive performance over MOGP in terms of mean ± standard deviation in testing accuracy of design variables.

Model	Loss Function	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, ${\mathit{t}}_{\mathbf{flight}}$ (%)	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, P (%)	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, $\mathit{\theta }$ (%)
cGANBC	BC	98.4 ± 0.858	99.0 ± 0.653	98.6 ± 0.424
cGANMSE	MSE	98.3 ± 0.843	98.8 ± 0.745	98.7 ± 0.426

displays the mean and standard deviation of the testing accuracy for cGAN models. The mean of both cGAN models shows that the cGAN has an overall $98\%$ accuracy. The low standard deviation of the cGAN model represents the robustness of the cGAN model. Both BC-based and MSE-based cGAN surrogates outperform the MOGP in terms of accuracy.

We compare and visualize an arbitrary set of predicted results among MSE-based cGAN, BC-based cGAN, and MOGP surrogates using the same design requirements as Table 4.

Table 7. Testing accuracy of visualized case for cGAN model BC and MSE.

Model	Loss Function	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, ${\mathit{t}}_{\mathbf{flight}}$ (%)	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, P (%)	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, $\mathit{\theta }$ (%)
cGANBC	BC	98.8	99.8	98.7
cGANMSE	MSE	98.9	99.7	98.7

presents the testing accuracy of the cGAN models for the visualization case. Figure 7 shows the optimal takeoff trajectory profiles for the cGAN models, with BC and MSE having similar predictive accuracy. The results indicate that the cGAN outperforms the MOGP’s predictive accuracy in this visualization case. Specifically, the cGAN captures the general trend of ground truth well due to the predictive power of the DNN surrogate as well as the addition of the discriminator to drive prediction matching training data set shapes. However, there are still some unexpected wiggles when we look at the acceleration profile. In addition, the cGAN approximates ground truth based on an average of a number of predictions (100 predictions in this work) over the same set of input parameters (design requirements in this work), where the prediction may vary slightly due to Monte Carlo properties. Hence, we develop and introduce the regGAN surrogate for further predictive improvement as follows.

3.3. regGAN Surrogate Modeling

The predictive performance of the regGAN surrogates is compared against Dymos simulation-based optimal trajectory predictions as well as MOGP and cGAN surrogates. The same 1000 training samples and 300 testing samples are used for regGAN model training and verification, respectively. The architecture and the hyperparameters are introduced in Section 2.5.3. We first focus on the regGAN surrogate with a single MSE loss function by setting a zero weight on BC in the combined loss (CL1) (

Table 8. Mean ± standard deviation of testing accuracy for regGAN model CL1.

Model	${\mathit{w}}_{\mathbf{MSE}}$	${\mathit{w}}_{\mathbf{BC}}$	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, ${\mathit{t}}_{\mathbf{flight}}$ (%)	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, P (%)	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, $\mathit{\theta }$ (%)
CL1	1	0	99.5 ± 0.246	99.7 ± 0.196	99.7 ± 0.117

). The table shows that the regGAN CL1 model reaches 99.5% accuracy with robust predicted results, which can be recognized by the low standard deviation of testing accuracy. To further compare results, we implement regGAN on the same visualization case in Table 4.

Table 9. Testing accuracy of visualization case for regGAN model BC and MSE.

Model	${\mathit{w}}_{\mathbf{MSE}}$	${\mathit{w}}_{\mathbf{BC}}$	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, ${\mathit{t}}_{\mathbf{flight}}$ (%)	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, P (%)	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, $\mathit{\theta }$ (%)
CL1	1	0	99.4	99.8	99.8

shows that the regGAN has a mean testing accuracy of over 99.5% in the visualization case. Figure 8 shows that the optimal takeoff trajectory profiles predicted by regGAN match Dymos results closer than the cGAN model BC and MSE and MOGP model 2, as expected.

To investigate the regGAN surrogate’s performance, we utilized a combination of two loss functions during the model training.

Table 10. regGAN model conditions and mean ± standard deviation of testing accuracy for design variables.

Model	${\mathit{w}}_{\mathbf{MSE}}$	${\mathit{w}}_{\mathbf{BC}}$	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, ${\mathit{t}}_{\mathbf{flight}}$ (%)	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, P (%)	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, $\mathit{\theta }$ (%)
CL2	1	0.01	99.7 ± 2.388	99.7 ± 0.118	99.7 ± 0.123
CL3	1	0.001	99.7 ± 0.193	99.7 ± 0.160	99.7 ± 0.124
CL4	1	0.0001	99.6 ± 0.216	99.7 ± 0.145	99.7 ± 0.128

shows the loss weights together with the mean and standard deviation of testing accuracy for each regGAN model. The results indicate that all regGAN CL models have over 99.6% mean testing accuracy and high predictive robustness, as revealed by the low standard deviations.

Table 11. Testing accuracy for visualization case of regGAN model CL2–4.

Model	${\mathit{w}}_{\mathbf{MSE}}$	${\mathit{w}}_{\mathbf{BC}}$	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, ${\mathit{t}}_{\mathbf{flight}}$ (%)	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, P (%)	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, $\mathit{\theta }$ (%)
CL2	1	0.01	99.6	99.9	99.8
CL3	1	0.001	99.9	99.8	99.8
CL4	1	0.0001	99.9	99.9	99.5

shows that the testing accuracy by regGAN models CL2–4 on the visualization case agrees well with the corresponding mean testing accuracy. Figure 9 shows that all predicted trajectories by regGAN models CL2–4 match the reference profiles, meaning that an accuracy of 99.6% is reliable in lieu of simulation models. regGAN models CL2–4 slightly outperform regGAN model CL1, but the matches towards reference profiles are at a comparable level.

We also conducted a parametric study for different combined loss weights to further investigate regGAN performance. We use a different number of training samples, ranging from 50, 100, 200, 400, 600, 800, to 1000.

Table 12. Parametric study for regGAN models exhibits deeper insights on mean ± standard deviation of testing accuracy with respect to the number of training samples for each design variable group. Note that we only vary $w_{\mathrm{BC}}$ while keeping $w_{\mathrm{MSE}}$ as 1.

Model	${\mathit{w}}_{\mathbf{BC}}$	Samples	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, ${\mathit{t}}_{\mathbf{flight}}$ (%)	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, P (%)	${\overline{\mathit{ϵ}}}_{{\mathit{L}}_1}$, $\mathit{\theta }$ (%)
CL1	0	50	95.2 ± 0.279	97.4 ± 1.58	98.3 ± 0.758
		100	98.6 ± 0.906	98.5 ± 1.05	98.9 ± 0.414
		200	99.3 ± 0.368	99.2 ± 0.461	99.4 ± 0.212
		400	99.5 ± 0.314	99.5 ± 0.389	99.6 ± 0.167
		600	99.2 ± 0.477	99.4 ± 0.294	99.6 ± 0.148
		800	99.7 ± 0.281	99.7 ± 0.146	99.7 ± 0.119
		1000	99.5 ± 0.246	99.7 ± 0.196	99.7 ± 0.117
CL2	0.01	50	97.4 ± 1.51	97.7 ± 0.852	97.6 ± 0.955
		100	97.4 ± 0.732	98.6 ± 0.829	98.9 ± 0.434
		200	99.3 ± 0.443	97.8 ± 0.980	98.9 ± 0.427
		400	99.6 ± 0.251	99.6 ± 0.257	99.6 ± 0.145
		600	99.6 ± 0.204	99.7 ± 0.215	99.6 ± 0.153
		800	99.5 ± 0.314	99.7 ± 0.158	99.6 ± 0.147
		1000	99.7 ± 0.239	99.7 ± 0.118	99.7 ± 0.123
CL3	0.001	50	95.3 ± 2.77	97.7 ± 1.29	98.0 ± 0.759
		100	97.2 ± 1.72	96.4 ± 2.08	96.6 ± 1.31
		200	99.3 ± 0.380	99.0 ± 0.557	99.1 ± 0.308
		400	99.6 ± 0.343	99.4 ± 0.378	99.4 ± 0.230
		600	99.7 ± 0.228	98.9 ± 0.767	99.4 ± 0.239
		800	99.5 ± 0.228	99.7 ± 0.194	99.7 ± 0.127
		1000	99.7 ± 0.193	99.7 ± 0.160	99.7 ± 0.124
CL4	0.0001	50	97.1 ± 1.64	98.5 ± 1.11	96.4 ± 1.78
		100	98.5 ± 0.757	98.9 ± 0.648	97.7 ± 0.847
		200	99.2 ± 0.384	98.8 ± 0.706	98.9 ± 0.494
		400	99.4 ± 0.346	99.4 ± 0.358	99.3 ± 0.249
		600	99.2 ± 0.611	99.3 ± 0.509	99.4 ± 0.195
		800	99.6 ± 0.287	99.5 ± 0.363	99.5 ± 0.182
		1000	99.6 ± 0.216	99.7 ± 0.145	99.7 ± 0.128

indicates that using 100 and 200 training samples is able to obtain the above with an overall 97% accuracy. Note that when provided with 50 training samples, CL1 has the lowest mean testing accuracy and highest standard deviation time $t_{\mathrm{flight}}$ prediction, while CL2, CL3, and CL4 have better overall predictive performance. Moreover, regGAN model CL2 ($w_{\mathrm{MSE}}=1$ and $w_{\mathrm{BC}}=0.01$) achieves over 99.5% accuracy starting with 400 training samples, while all the other regGAN surrogates require at least 800 training samples. In addition, regGAN model CL2 consistently shows lower predictive standard deviations (such as 0.251%, 0.257%, and 0.145% on $t_{\mathrm{flight}}$, $\mathbf{P}$, and $\mathit{\theta }$) on the testing data set, further confirming a better and more robust predictive performance. Based on the visualization case, a mean testing accuracy of over 99.5% can be considered sufficiently accurate with negligible differences.

In sum, the MOGP surrogates could not match ground truths or capture the general trend of optimal takeoff trajectories well using 1000 training samples, which may be due to the Gaussian assumption of GP series models. The cGAN surrogates achieve better performance over MOGP using the 1000 training samples, with a closer match towards the general trend of ground truth labels, mainly because of the guidance of the discriminator for generating similar data patterns as the training data set. The regGAN surrogates outperform MOGP and cGAN since regGAN is trained with a combined loss function of MSE and BC adversarial losses. The BC adversarial loss leads regGAN to handle the general trend of observations, while the MSE loss directly drives the match between predictions and observations.

4. Conclusions

In this paper, we investigated surrogate-based optimal takeoff trajectory predictions for electric vertical takeoff and landing (eVTOL) drones within the scope of urban air mobility. We developed the regression generative adversarial network (regGAN), which outperformed the Gaussian process (MOGP) and the conditional generative adversarial network (cGAN) by achieving over 99.5% accuracy. We summarize the main contribution of this work as follows.

First, we implemented the surrogate-based inverse mapping concept into eVTOL optimal trajectory design for the first time. Specifically, surrogate models took design requirements as input and predicted optimal trajectories. We realized fast interactive eVTOL takeoff trajectory design without running any optimizations since the trained surrogates directly predicted optimal trajectories. However, reducing training costs is essential since each training sample requires a simulation-based trajectory design optimization.

Second, we introduced the MOGP, a representative traditional surrogate model, into the eVTOL takeoff trajectory design for rapid predictions. The results showed that the MOGP with a square exponential kernel function could accurately capture the inverse mapping using 1000 training samples. We then implemented a cGAN for eVTOL inverse mapping since cGAN also makes use of the GAN architecture for regression tasks. The results revealed that the cGAN achieved over 98% generalization accuracy in predicting optimal designs using the same 1000 training samples as the MOGP, which means the cGAN outperformed the MOGP on predictive performance (around 92% accuracy). In addition, the visualization case verified that the cGAN could match the general trend of optimal designs well with actual observations from Dymos but missed detailed features.

Third, we introduced the regGAN into takeoff trajectory design for the first time and achieved over 99.6% accuracy in predicting optimal design variables with the same 1000 training samples as the MOGP and cGAN. By varying the weights of different loss functions, the regGAN could achieve over 99.6% accuracy. Moreover, results indicated that the best regGAN surrogate architecture consistently achieved over 99.5% accuracy if provided with 400 or more random training samples. This confirmed the outstanding predictive performance and potential generality of the regGAN.

In future work, we are planning to explore and develop other novel deep learning architectures in regGAN. In addition, the simulation models used in this work are not high-fidelity models but effective for describing the physics; we will increase the fidelity of simulation models in future work, which may lead to a higher training cost for surrogate modeling. Moreover, we will consider takeoff time as another constraint to make sure the total takeoff will not take unreasonable time.