Online Deterministic 3D Trajectory Generation for Electric Vertical Take-Off and Landing Aircraft

Abstract

The use of non-piloted eVTOL aircraft in non-segregated airspace requires reliable and deterministic automatic flight guidance systems for the aircraft to remain predictable to all the users of the airspace and maintain a high level of safety. In this paper we present a 3D trajectory generation module based on clothoid transition segments in the horizontal plane and high order polynomial transition segments in the vertical plane. The expressions of the coefficients of the polynomial are derived offline are used to generate the trajectory online, making the system capable of running in real time without requiring enormous computational power. For the horizontal plane, we focus on the flyby transition, and therefore present a thorough analysis of the flyby geometry and the limitations linked to this geometry and the construct of three-segment trajectory generation around a fixed turn rate. We present feasible solutions for these limitations, and show simulation results for the combined horizontal and vertical plane concepts, allowing the system to generate complex 3D trajectories.

1. Introduction

Electric Vertical Take-Off and Landing (eVTOL) aircraft have received much attention recently for their versatility and wide range of applications. This includes urban air mobility, medical evacuations, firefighting, and so on. In this framework, there is an increasing need for reliable automatic flight guidance systems for both manned and unmanned eVTOL aircraft. The reliability requirement for trajectory generation and control implies the ability to withstand the strict certification procedure of safety-critical systems. Furthermore, the systems are intended for use in non-segregated airspace. Therefore, they are expected to be predictable to the air traffic controllers, and to other users of the airspace, while maintaining the versatility expected of a human pilot. The development of such systems remains an active research topic in the international scientific community.

Current literature mainly focuses on non-deterministic and sampling-based trajectory generation methods, such as RRT or RRT* [1,2,3,4]. Zhao et al. [5] explore various computational intelligence-based methods for path planning including Fuzzy logic, neuron computing, genetic computing, and probabilistic computing. These approaches are promising for small UAVs operated in secluded and safe environments. However, as they are all probability-based methods, they are non-deterministic and are, therefore, unsuitable for eVTOL aircraft operated in non-segregated airspace.

Qu et al. [6] consider a Reinforcement learning-based grey wolf optimizer (RLGWO) for optimal path planning and use a cubic B-spline to ensure a continuous trajectory. This hybrid approach comes with its unique set of advantages, but runs into the same limitations as the other non-deterministic methods.

Al Satai et al. [7] propose using any-angle path-finding algorithms like Basic Theta*, Lazy Theta*, and Phi* paired with Bézier curves to smooth the generated paths to be converted to a flyable trajectory. This is a novel approach to efficiently construct optimal trajectories, but may be limited by its computational complexity for real-time performance. The approach presented in [8] looks at real-time path planning based on bi-level programming (BLP).The approach aims to integrate various path planning requirements, including convergence to a target, obstacle avoidance, path length optimization, flight path smoothing, and adaptability, to changes in UAV kinematic and sensory properties. However, it is also limited by the computational complexity.

Bousson et al. [9] address the 4D trajectory generation problem. In addition to the geographic location, they also consider the arrival time at the waypoints. They construct an optimization problem in which Lagrange and Chebyshev polynomials approximate the solution space. Two simulations verified that the approach is able to generate smooth trajectories. However, high-order polynomials are required to capture the system’s nonlinear behavior, and they are used in the optimization problem. This makes the algorithm infeasible for online use.

In [10], deterministic path planning for a forklift truck is investigated. Given a specified number of waypoints, linear line segments are combined with polar splines. The splines are constrained on the radius, slope, and curvature to achieve 2-order continuity of the generated trajectory. However, the approach only covers geometrical considerations. Characteristics of the dynamical system are neglected during the trajectory design process.

Schneider et al. [11,12] use line and circular arc segments to construct the trajectory. For the transition between segments, a clothoid function is used to ensure $C^2$-order continuity of the trajectory. The clothoid function is designed based on the aircraft characteristics, namely the roll dynamics. While the approach is proven to work well for fixed-wing aircraft, flight path construction is more restrictive when considering eVTOL configurations. The clothoid-based approach is widely used for its simplicity and the system presented in [11] has been successfully flight tested on a fixed-wing aircraft [13].

In this paper, we conduct an analysis of flyby maneuvers and the implications of the geometrical construction of these maneuvers and the resulting limitations, and propose an extension to the system developed in [11] to make it suitable for an eVTOL aircraft. One aspect analyzed is the turn distance. This is important because it controls how far apart two successive waypoints need to be. This is particularly important for eVTOL aircraft as they are typically intended to be flown in urban environments and other environments with short distances. The other aspect is the inherent limitation of the allowable maximum angle between two successive legs in the flight plan. We show that this limitation applies to trajectory generation systems designed by combining straight lines, transition segments, and circular arcs. Then, we present an analytical derivation of the problem, exposing the variables and the potential solutions. We pursue the solution involving a dynamic change of the design turn rate. This allows the system to accommodate a wider range of flight plans.

In addition to the clothoid-based generation for lateral trajectory, we use a higher order polynomial for the vertical trajectory, where, instead of the optimization framework presented in [9], the expressions for coefficients of the polynomial are calculated analytically beforehand, taking into account predefined continuity requirements, and no optimization is performed. This allows us to maintain a low computational complexity.

As shown in the literature review above, most of the current literature considers trajectory optimization, which not only requires considerable computational power, but also is difficult to certify. The purpose of this work is to present a system capable of generating deterministic 3D trajectories online with a low computational cost. Therefore, rather than considering the optimality of the generated path, we only consider the determinism. This is essential for use in non-segregated airspace where the behavior of the aircraft needs to be predictable to other users of the airspace. Additionally, this makes the system relatively easy to certify.

This paper is organized as follows: Section 2 provides a complete mathematical description of each element of the trajectory composition (straigt line, circular arc segment, and clothoid segment). Section 3 deals with the construction of the trajectory by merging the mathematical preliminaries with the aircraft dynamics. Section 4 details the flyby implications for eVTOLs. Finally, Section 5 illustrates our key findings and analysis.

2. Trajectory Construction: Overview

It is widely known that the shortest path between two points is a straight line. Ref. [14] shows that, if there are constraints that need to be satisfied, such as a given orientation at both the starting and end points, then the shortest path is a combination of straight lines and circular arcs also known as a Dubin’s path. This is the approach followed in [15], which is the basis of the work carried out in this paper.

Since the trajectory is made up of straight lines and circular arcs, there is a curvature discontinuity arising at the joints of the two segments. As shown in [16], a straight line has a curvature $\kappa =0$, while a circle has a curvature inversely proportional to its radius $\kappa =1/r$. This is also shown later in Section 2.2 and Section 2.3 of this paper.

A clothoid or Euler spiral is a curve with the property that its curvature linearly increases with the arc length [17,18]. This makes the clothoid a particularly good candidate for transition segments that connect the straight line and circle as shown in Figure 1 where the straight line is shown in blue, the clothoid in green and the circular arc in orange and stars are delimiting different sections.

In the following section we will first recapitulate how the curvature of a plane curve is calculated. In Section 2.2, Section 2.3 and Section 2.4, the curve parametrization will be given, and the curvature will be calculated, for each of the elements in Figure 1 (line, clothoid, circle).

2.1. On Plane Curves

This section introduces important mathematical preliminaries for plane curves. These mathematical relations will be used throughout this paper.

In general, a curve in two-dimensional Euclidean space can be parametrized as a piecewise differentiable function [19]:

$$\overrightarrow{c}=\left[\begin{array}{c}{c}_x\\ c_y\end{array}\right]=\overrightarrow{c}\left(\tau \right)\in {\mathbb{R}}^2$$

(1)

where $\tau$ is the running parameter of the parametrization.

The arc length s of the curve between two instants ${\tau }_1$ and ${\tau }_2$ can be calculated as [20]:

$$s\left({\tau }_1,{\tau }_2\right)={\int }_{{\tau }_1}^{{\tau }_2}\left|\right|\overrightarrow{\dot{c}}\left(\tau \right)\left|\right|d\tau$$

(2)

where $\left|\right|\overrightarrow{\dot{c}}\left(\tau \right)\left|\right|=\sqrt{{\dot{c}}_x{\left(\tau \right)}^2+{\dot{c}}_y{\left(\tau \right)}^2}$ is the 2-norm of the first derivative of the curve with respect to the running parameter.

According to [21], the pointwise curvature of a curve $\overrightarrow{c}\left(\tau \right)\in {\mathbb{R}}^2$ is the norm of the derivative of the unit tangent vector $\overrightarrow{T}$ at that point with respect to the arc length:

$$\kappa \left(\tau \right)=\left|\left|\frac{d\overrightarrow{T}}{ds}\right|\right|=\left|\left|\overrightarrow{T^{\prime }}\left(s\right)\right|\right|$$

(3)

where s is the arc length. The derivative of the unit tangent vector with respect to the arc length can be found as follows:

We know the tangent vector $\overrightarrow{t}$ is simply the first derivative of the position vector with respect to time, or in our case, the running parameter $\tau$:

$$\overrightarrow{t}=\frac{d\overrightarrow{c}}{d\tau }=\overrightarrow{\dot{c}}\left(\tau \right)$$

(4)

And the unit tangent vector with respect to $\tau$ is:

$$\overrightarrow{T}\left(\tau \right)=\frac{\overrightarrow{t}}{\left|\left|\overrightarrow{t}\right|\right|}=\frac{\overrightarrow{\dot{c}}\left(\tau \right)}{\left|\left|\overrightarrow{\dot{c}}\left(\tau \right)\right|\right|}$$

(5)

If we take the derivative of Equation (2) with respect to $\tau$, we obtain:

$$\frac{ds}{d\tau }=\left|\left|\overrightarrow{\dot{c}}\left(\tau \right)\right|\right|$$

(6)

Replacing the result of Equation (6) into Equation (5), we obtain:

$$\overrightarrow{\dot{c}}\left(\tau \right)=\overrightarrow{T}\left(\tau \right)·\frac{ds}{d\tau }$$

(7)

Using the chain rule, we can rewrite Equation (3) in terms of $d\tau$ as follows:

$$\kappa \left(\tau \right)=\left|\left|\frac{d\overrightarrow{T}}{ds}\right|\right|=\left|\left|\frac{d\overrightarrow{T}}{d\tau }·\frac{d\tau }{ds}\right|\right|$$

(8)

With the information we have from Equation (6), we obtain:

$$\kappa \left(\tau \right)=\left|\left|\frac{d\overrightarrow{T}}{d\tau }·\frac{d\tau }{ds}\right|\right|=\left|\left|\frac{d\overrightarrow{T}}{d\tau }·\frac{1}{\left|\left|\overrightarrow{\dot{c}}\left(\tau \right)\right|\right|}\right|\right|=\frac{\left|\left|\overrightarrow{\dot{T}}\left(\tau \right)\right|\right|}{\left|\left|\overrightarrow{\dot{c}}\left(\tau \right)\right|\right|}$$

(9)

Another representation of the curvature that will prove useful is derived in [16], and shown in Equation (10). This gives us the angle $\phi$ representing the change in direction covered over an arc length s.

$$\kappa \left(\tau \right)={\left.\frac{d\phi }{ds}\right|}_{\tau }$$

(10)

Now that we have derived the mathematical expressions for curvature and arc length, we can apply them to the segments that define our trajectory, starting with straight line segments.

2.2. On Straight Line Segments

From the fundamental principles of geometry in [22], we can deduce the parametrization of a straight line followed by an aircraft flying with a course angle ${\chi }_T^R$ as shown in Figure 2:

$$\overrightarrow{c}\left(\tau \right)=\left[\begin{array}{c}{c}_{x_{line}}\\ c_{y_{line}}\end{array}\right]=\left[\begin{array}{c}{x}_i+\tau ·cos\left({\chi }_T^R\right)\\ y_i+\tau ·sin\left({\chi }_T^R\right)\end{array}\right]$$

(11)

where $\tau$ is the distance traveled on the line, measured from the coordinates of starting position of the line given by $x_i$ and $y_i$.

Curvature and Arc Length of the Line Segment

We can apply Equation (9) to find the curvature. We therefore need the first derivatives of the unit tangent vector and the parametrization of the line itself.

The derivative of the line curve is:

$$\overrightarrow{\dot{c}}\left(\tau \right)=\left[\begin{array}{c}cos\left({\chi }_T^R\right)\\ sin\left({\chi }_T^R\right)\end{array}\right]$$

(12)

The tangent vector is calculated using Equation (5):

$$\overrightarrow{T}\left(\tau \right)=\left[\begin{array}{c}\frac{{\dot{c}}_x}{\sqrt{{\dot{c}}_x^2+{\dot{c}}_y^2}}=cos\left({\chi }_T^R\right)\\ \frac{{\dot{c}}_y}{\sqrt{{\dot{c}}_x^2+{\dot{c}}_y^2}}=sin\left({\chi }_T^R\right)\end{array}\right]$$

(13)

We can now immediately see that, as expected, the curvature calculated with Equation (9) is zero as $\overrightarrow{\dot{T}}\left(\tau \right)={[0,0]}^T$. Additionally, as we mentioned before that $\tau$ represent the distance traveled along the line, it also makes it the arc length, such that:

$$s_{line}(0,\tau )=\tau$$

(14)

2.3. On Circle Segment

A circle centered at $(x_i,y_i)$ as shown in Figure 3 can be parametrized as [23]:

$$\overrightarrow{c}\left(\tau \right)=\left[\begin{array}{c}{c}_{x_c}\\ c_{y_c}\end{array}\right]=\left[\begin{array}{c}{x}_i+r_c·cos\left(\tau \right)\\ y_i+r_c·sin\left(\tau \right)\end{array}\right]$$

(15)

where $x_i$ and $y_i$ are the coordinates of circle center and $r_c$ is the radius of the circle.

2.3.1. Curvature of the Circle Segment

Applying Equation (9) we find the curvature.

The derivative of the circle parametrization with respect to $\tau$ is:

$$\overrightarrow{\dot{c}}\left(\tau \right)=\left[\begin{array}{c}-r_c·sin\left(\tau \right)\\ r_c·cos\left(\tau \right)\end{array}\right]$$

(16)

The tangent vector is calculated using Equation (5):

$$\overrightarrow{T}\left(\tau \right)=\left[\begin{array}{c}\frac{{\dot{c}}_x}{\sqrt{{{\dot{c}}_x}^2+{{\dot{c}}_y}^2}}=-sin\left(\tau \right)\\ \frac{{\dot{c}}_y}{\sqrt{{{\dot{c}}_x}^2+{{\dot{c}}_y}^2}}=cos\left(\tau \right)\end{array}\right]$$

(17)

And its derivative with respect to $\tau$ is:

$$\overrightarrow{\dot{T}}\left(\tau \right)=\left[\begin{array}{c}-cos\left(\tau \right)\\ -sin\left(\tau \right)\end{array}\right]$$

(18)

Replacing these solutions in Equation (9) we find the curvature:

$${\kappa }_c=\frac{\left|\left|\left[\begin{array}{c}-cos\left(\tau \right)\\ -sin\left(\tau \right)\end{array}\right]\right|\right|}{\left|\left|\left[\begin{array}{c}-r_c·sin\left(\tau \right)\\ r_c·cos\left(\tau \right)\end{array}\right]\right|\right|}=\frac{1}{r_c}$$

(19)

2.3.2. Arc Length of the Circle Segment

The arc length of the circle segment parametrized as shown in Equation (15) is calculated using Equation (2):

$$s_c\left({\tau }_1,{\tau }_2\right)={\int }_{{\tau }_1}^{{\tau }_2}\left|\left|\left[\begin{array}{c}-r_c·sin\left(\tau \right)\\ r_c·cos\left(\tau \right)\end{array}\right]\right|\right|d\tau =r_c·({\tau }_2-{\tau }_1)$$

(20)

Assuming ${\tau }_1=0$ and ${\tau }_2=\tau$, we obtain:

$$s_c(0,\tau )=r_c·\tau$$

(21)

2.4. On Clothoid Segment

As shown in [15,24], an exact parametrization of the clothoid curve, as shown in Figure 4, is given by the irreducible Fresnel integrals:

$$\overrightarrow{c}\left(\tau \right)=\left[\begin{array}{c}{c}_{x_{cl}}\\ c_{y_{cl}}\end{array}\right]=\left[\begin{array}{c}{x}_i+A·{\int }_0^{\tau }cos\left({\tau }^2\right)d\tau \\ y_i+A·{\int }_0^{\tau }sin\left({\tau }^2\right)d\tau \end{array}\right]$$

(22)

where A is the shaping parameter, which is further explained in Section 2.4.4, and $x_i$ and $y_i$ are the starting position of the clothoid curve.

The goal is to have a clothoid section starting with a curvature $\kappa =0$ and ending with the curvature of the circular arc used for the maneuver $\kappa =\frac{1}{r_c}$.

2.4.1. Curvature of the Clothoid Segment

To calculate the curvature at a given instant, we start by taking the first derivative of the clothoid curve with respect to $\tau$:

$$\overrightarrow{\dot{c}}\left(\tau \right)=\left[\begin{array}{c}A·cos\left({\tau }^2\right)\\ A·sin\left({\tau }^2\right)\end{array}\right]$$

(23)

Then we find the unit tangent vector from Equation (5):

$$\overrightarrow{T}\left(\tau \right)=\left[\begin{array}{c}\frac{{\dot{c}}_x}{\sqrt{{\dot{c}}_x^2+{\dot{c}}_y^2}}=cos\left({\tau }^2\right)\\ \frac{{\dot{c}}_y}{\sqrt{{\dot{c}}_x^2+{\dot{c}}_y^2}}=sin\left({\tau }^2\right)\end{array}\right]$$

(24)

Taking the derivative of the unit tangent vector with respect to $\tau$ we obtain:

$$\overrightarrow{\dot{T}}\left(\tau \right)=\left[\begin{array}{c}-2\tau sin\left({\tau }^2\right)\\ 2\tau cos\left({\tau }^2\right)\end{array}\right]$$

(25)

We substitute $\overrightarrow{\dot{c}}\left(\tau \right)$ and $\overrightarrow{\dot{T}}\left(\tau \right)$ into Equation (9) and obtain the parametrized curvature of the clothoid:

$${\kappa }_{cl}\left(\tau \right)=\frac{\left|\left|\left[\begin{array}{c}-2\tau sin\left({\tau }^2\right)\\ 2\tau cos\left({\tau }^2\right)\end{array}\right]\right|\right|}{\left|\left|\left[\begin{array}{c}A·cos\left({\tau }^2\right)\\ A·sin\left({\tau }^2\right)\end{array}\right]\right|\right|}=\frac{2·\tau }{A}$$

(26)

2.4.2. Arc Length of the Clothoid Segment

We can now calculate the arc length:

$$s_{cl}\left({\tau }_1,{\tau }_2\right)={\int }_{{\tau }_1}^{{\tau }_2}\left|\left|\left[\begin{array}{c}A·cos\left({\tau }^2\right)\\ A·sin\left({\tau }^2\right)\end{array}\right]\right|\right|d\tau =A·({\tau }_2-{\tau }_1)$$

(27)

Assuming that ${\tau }_1=0$ and ${\tau }_2=\tau$, we obtain:

$$s_{cl}(0,\tau )=A·\tau$$

(28)

Equation (28) shows the simplicity of the clothoid curve, which lies in the fact that it can be defined completely by the shaping parameter A.

2.4.3. Change in Direction

Now that we have the arc length and curvature of the clothoid segment, we can use the second representation of curvature shown by Equation (10), to calculate the angle representing the change in direction throughout the clothoid segment:

$$\frac{2\tau }{A}={\left.\frac{d\phi }{d\left(A\tau \right)}\right|}_{\tau }$$

(29)

This leads to:

$${\left.d\phi \right|}_{\tau }=\frac{2\tau }{A}·A·d\tau$$

(30)

Integrating Equation (30) with respect to $\tau$, we obtain a very simple representation of change in angle, which is really a change in tge course angle incurred while flying the clothoid segment:

$${\phi }_{cl}\left(\tau \right)={\tau }^2$$

(31)

2.4.4. Clothoid Approximation

As shown earlier in this section, the clothoid can be exactly represented using Fresnel integrals. The Fresnel integrals are transcendental functions that cannot be solved analytically [25]. Given that the system is intended to be used for real-time computations, ref. [15] introduce a power series approximation of the clothoid:

$$\left[\begin{array}{c}\Delta x_{cl}\left(\tau \right)\\ \Delta y_{cl}\left(\tau \right)\end{array}\right]=\left[\begin{array}{c}A·\sum _{m=0}^5\frac{{(-1)}^m}{\left(2m\right)!·(4m+1)}·{\tau }^{4m+1}\\ A·\sum _{m=0}^5\frac{{(-1)}^m}{(2m+1)!·(4m+3)}·{\tau }^{4m+3}\end{array}\right],$$

(32)

As shown in [15], the fifth-order power series is accurate enough for trajectory generation applications. A comparison of the exact representation and the approximation is shown in Figure 5. As we can see, the power series approximation is accurate even for extreme cases of a change in the angle of ${180}^{\circ }$. Note that the clothoid shaping parameter A seen in Equation (32) simply controls the rate of growth of the clothoid curve.

In this section we presented the generalities of trajectory generation, and introduced the clothoid segment that connects straight lines and circular arcs. Since the clothoid is represented by Fresnel integrals that are rather difficult to solve online, we use a fifth-order power series approximation. The simplicity of the clothoid segment lies in the fact that it has only a single variable, which is the shaping parameter A. This variable will be determined in the next section, and connected to the aircraft dynamics.

3. Trajectory Construction: Flyby

The flyby transition is mainly based on the developments in [15]. A geometrical overview is shown in Figure 6, where the clothoid segment is shown in green. The flyby transition has three segments: the turn-in, circle, and turn-out. In the horizontal plane, the turn-in and turn-out segments are clothoid sections and are assumed to be symmetrical. As introduced earlier in this paper, they are used to provide continuous curvature in the generated trajectory. In this section we want to derive the geometry of the maneuver as depicted in Figure 6 and use it, together with the relations derived in Section 2 for the three different segments (line, clothoid, circle), to derive an expression for the turn distance (see Figure 7) that only depends on the following parameters:

• Fixed parameters resulting from the flight plan:
  course angle change between two subsequent legs $\Delta {\chi }_T$, kinematic speed $V_T$;
• Fixed parameters resulting from the inner loop controller design, the resulting inner loop dynamics and aircraft performance:
  roll time constant $T_p$, roll rate command $p_{cmd}$;
• Trajectory design parameter:
  desired turn rate ${\dot{\chi }}_{T_d}$.

Figure 6. Flyby geometry—an overview.

Figure 6. Flyby geometry—an overview.

Note that the desired turn rate above is strictly positive as the direction of the turn is defined by the positions of waypoints in the flight plan. See Figure 6 as an example of a right turn, which is independent of the turn rate. This direction is multiplied directly to the turn rate command generated by the trajectory generation module as 1 for a right turn and $-1$ for a left turn.

3.1. Horizontal Plane

From an aircraft perspective, there are two main parts to the maneuver:

1.

Straight line: the aircraft is flying in a straight line with a kinematic bank angle $\mu =0$.

2.

Circle: The aircraft is performing a turn with a kinematic bank angle as given by Equation (33)

$$\mu =arctan\left(\frac{V_T·{\dot{\chi }}_{T_d}}{g_0}\right)$$

(33)

where $V_T$ is the kinematic speed, ${\dot{\chi }}_{T_d}$ is the desired turn rate, and $g_0$ is the gravitational acceleration.

The desired turn rate can be used to describe the radius of the maneuver circle:

$$r_c=\frac{V_T}{{\dot{\chi }}_{T_d}}$$

(34)

As we see, in one instant, we have a zero kinematic bank angle, and in the next instant, we require a non-zero bank angle in order to have that non-zero turn rate from Equation (33). This discontinuity is represented geometrically by the step in curvature of the trajectory. Therefore, the clothoid segment should help the aircraft build up the necessary bank angle, in order to achieve the desired turn rate. [15] describes a method to find the time needed for the aircraft to build up a bank angle $\mu$, given a constant roll rate p, and considering the aircraft turn performance through the roll time constant $T_p$. This is shown in Equation (35).

$$t_{cl}=2·T_p+\frac{{\mu }_{cmd}}{p_{cmd}}$$

(35)

We know the clothoid arc length from Equation (28). And from the basic physics of motion, given a distance and a speed, we can calculate the time it takes to travel that distance. This gives us a relationship between the transition time $t_{cl}$ and the clothoid shaping parameter A.

$$t_{cl}=\frac{s(0,{\tau }_{cl})}{V_T}=\frac{A·{\tau }_{cl}}{V_T}$$

(36)

where the arc length $s_{cl}(0,{\tau }_{cl})=A·{\tau }_{cl}$ calculated in Equation (28) was inserted.

Solving Equation (36) for the running parameter ${\tau }_{cl}$, allows us to express it as a function of the shaping parameter A and the physical transition time $t_{cl}$ during which we need to fly a clothoid segment.

$${\tau }_{cl}=\frac{V_T·t_{cl}}{A}$$

(37)

As we stated before, our goal is to have a clothoid section that starts with a curvature $\kappa =0$ and ends with a curvature $\kappa =\frac{1}{r_c}$, which is the curvature of the maneuver circle calculated by Equation (34). This means that we want the curvature of the clothoid, which we can parametrize with Equation (26), to equal $\frac{1}{r_c}$:

$${\kappa }_c=\frac{1}{r_c}\equiv {\kappa }_{cl}=\frac{2·{\tau }_{cl}}{A}$$

(38)

Replacing ${\tau }_{cl}$ in Equation (38) with the expression for ${\tau }_{cl}$ given in Equation (36) we obtain:

$$\frac{1}{r_c}=\frac{2·V_T·t_{cl}}{A^2}$$

(39)

We can now solve Equation (39) for A:

$$A=\sqrt{2·V_T·r_c·t_{cl}}$$

(40)

Substituting Equation (35), we obtain the final expression for the shaping parameter A:

$$A=\sqrt{2·V_T·r_c·\left(2·T_p+\frac{{\mu }_{cmd}}{p_{cmd}}\right)}$$

(41)

where ${\mu }_{cmd}$ and $r_c$ are calculated using Equation (33) and Equation (34), respectively.

Equation (41) gives us a parametrization of the clothoid shaping parameter A as a function of the aircraft dynamics and performance. As shown in Figure 5, this parameter controls the rate of growth of the clothoid. This means that, although the start and end points of the clothoid are defined by the curvatures of the straight line and maneuver circle, the space occupied by the clothoid, and therefore the distance flown by the aircraft, is dependent on the aircraft performance. And this can be traced back to Equation (36) which defines the time needed for the aircraft to build up a certain bank angle.

3.1.1. Turn Distance

Now that we have established some general features of the flyby, we need to calculate the turn distance $d_{turn}$ between the $TO$ waypoint and the first clothoid point $clp1$ as shown in Figure 7. This will be used to define the location ${\overrightarrow{r}}^{clp1}$ where the turn maneuver needs to start. The aim is to derive an expression that only depends on the parameters $V_T,\Delta {\chi }_T,T_p,p_{cmd}$, and the design parameter ${\dot{\chi }}_{T_d}$.

Figure 7. Turn distance.

Figure 7. Turn distance.

In Figure 7, we can already see that, geometrically, the turn distance is made up of the distance between the $TO$ and the auxiliary point $ap1$, and the distance $d_2$:

$$d_{turn}=d_{ap}+d_2$$

(42)

We see a triangle formed between the circle center $cc$, the $TO$ and the $ap1$ points. The left side of that triangle is exactly $d_{ap}$. The triangle has the angle $\overline{\beta }$ (see Equation (47) at the $ap1$ vertex and angle $\Delta {\chi }_T^{*}/2$ (due to the symmetry of the maneuver) at the $cc$ vertex. We can apply the law of sines [26] to find the $d_{ap}$ side:

$$d_{ap}=d_{cc}·\frac{sin\left(\frac{\Delta {\chi }_T^{*}}{2}\right)}{sin\left(\overline{\beta }\right)}$$

(43)

where $d_{cc}$ is the distance between the $TO$ and the circle’s center. This distance can also be calculated using the law of sines, considering that the angle at the $TO$ vertex is ${\alpha }_T/2$, and knowing that the side opposite to $d_{cc}$ is $r_c+d_1$:

$$d_{cc}=\left(r_c+d_1\right)·\frac{sin\left(\overline{\beta }\right)}{sin\left(\frac{{\alpha }_T}{2}\right)}$$

(44)

where ${\alpha }_T=\pi -\Delta {\chi }_T$ as can be seen in Figure 7 and where $d_1$ is the distance between the clothoid point $clp2$ and the auxiliary point $ap1$. This distance can be found by considering the right triangle formed between $clp2$, $ap1$, and the point marking the end of $d_3$. The side opposite to the $ap1$ vertex represents the lateral distance traveled throughout the clothoid segment $\Delta y_{cl,f}=\Delta y_{cl}({\tau }_{cl},A)$ (see Section 2.4.4). Keeping in mind that the angle opposite to $\beta$ is also equal to $\beta$ [27], we can find the $d_1$ distance considering the right triangle:

$$d_1=\frac{\Delta y_{cl,f}}{sin\left(\beta \right)}$$

(45)

In the expressions above, ${\alpha }_T$ is the angle between the legs of the flight plan and the angles $\beta$ and $\overline{\beta }$ are defined geometrically as:

$$\beta =\frac{\pi }{2}-{\phi }_{cl}$$

(46)

which is obtained from the angle sum of the right triangle defined by $ap1$, $cc$ and the perpendicular intersection to the $FROM-TO$ leg.

$$\overline{\beta }=\pi -\beta$$

(47)

which becomes clear when looking at the line connecting the $cc$ and $ap1$ points.

The distance $d_2$ needed in Equation (42) can be further broken down into a sum of $d_3$ and $d_4$. $d_3$ is a side of the right triangle explored earlier to find $d_1$. Since we already have $d_1$ and we know $\Delta y_{cl,f}$, we can apply the Pythagorean theorem [26] to find $d_3$:

$$d_3=\sqrt{d_1^2-\Delta y_{cl,f}^2}$$

(48)

The remaining distance, $d_4$ is geometrically the distance covered by the clothoid segment along the x-axis. It is therefore:

$$d_4=\Delta x_{cl}({\tau }_{cl},A)=\Delta x_{cl,f}$$

(49)

We can now substitute Equations (34) and (43)–(45) into Equation (42), as well as $d_2=d_3+d_4$ with $d_3$ and $d_4$ given by Equations (48) and (49). Then the resulting expression for $d_{turn}$ is:

$$\begin{array}{cc}\hfill d_{turn}& =\left(\frac{V_T}{{\dot{\chi }}_{T_d}}+\frac{\Delta y_{cl,f}}{sin\left(\beta \right)}\right)·\frac{sin\left(\frac{\Delta {\chi }_T^{*}}{2}\right)}{sin\left(\frac{{\alpha }_T}{2}\right)}+\Delta x_{cl,f}+\sqrt{{\left(\frac{\Delta y_{cl,f}}{sin\left(\beta \right)}\right)}^2-\Delta y_{cl,f}^2}\hfill \\ \hfill & =\left(\frac{V_T}{{\dot{\chi }}_{T_d}}+\frac{\Delta y_{cl,f}}{sin\left(\beta \right)}\right)·\frac{sin\left(\frac{\Delta {\chi }_T^{*}}{2}\right)}{sin\left(\frac{{\alpha }_T}{2}\right)}+\Delta x_{cl,f}+\frac{\Delta y_{cl,f}}{tan\left(\beta \right)}\hfill \end{array}$$

(50)

Figure 8 shows that the total change in course angle $\Delta {\chi }_T$ is the sum of the change in course incurred while flying the turn-in, circle, and turn-out segments:

$$\Delta {\chi }_T=\Delta {\chi }_T^{*}+2{\phi }_{cl}$$

(51)

If we now replace $\Delta {\chi }_T^{*}$ with $\Delta {\chi }_T-2{\phi }_{cl}$, $\beta$ with Equation (46), ${\alpha }_T=\pi -\Delta {\chi }_T$ and ${\phi }_{cl}={\tau }_{cl}^2$ (see Equation (31)) in Equation (50), we obtain an expression for $d_{turn}$ which only depends on $V_T,\Delta {\chi }_T,{\dot{\chi }}_{T_d},\Delta x_{cl,f},\Delta y_{cl,f}$, and ${\tau }_{cl}$:

$$d_{turn}=\left(\frac{V_T}{{\dot{\chi }}_{T_d}}+\frac{\Delta y_{cl,f}}{sin(\frac{\pi }{2}-{\tau }_{cl}^2)}\right)·\frac{sin\left(\frac{\Delta {\chi }_T-2{\tau }_{cl}^2}{2}\right)}{sin\left(\frac{\pi -\Delta {\chi }_T}{2}\right)}+\Delta x_{cl,f}+\frac{\Delta y_{cl,f}}{tan(\frac{\pi }{2}-{\tau }_{cl}^2)}$$

(52)

The clothoid distances $\Delta x_{cl,f}=\Delta x_{cl}({\tau }_{cl},A),\Delta y_{cl,f}=\Delta y_{cl}({\tau }_{cl},A)$ are given by Equation (32), where A, according to (41), only depends on $V_T,{\dot{\chi }}_{T_d},T_p,p_{cmd}$. Finally, we reformulate ${\tau }_{cl}$ by inserting Equations (34) and (39) into (37):

$${\tau }_{cl}=\frac{V_T{t}_{cl}}{\sqrt{2V_T\frac{V_T}{{\dot{\chi }}_{T_d}}{t}_{cl}}}=\sqrt{\frac{t_{cl}{\dot{\chi }}_{T_d}}{2}}$$

(53)

With $t_{cl}$ given by Equation (35) and ${\mu }_{cmd}$ by (33), we obtain an expression for ${\tau }_{cl}$ which only depends on $p_{cmd},T_p,{\dot{\chi }}_{T_d}$ and $V_T$:

$${\tau }_{cl}=\sqrt{\frac{{\dot{\chi }}_{T_d}}{2}\left(2T_p+\frac{arctan\left(\frac{V_T{\dot{\chi }}_{T_d}}{g_0}\right)}{p_{cmd}}\right)}$$

(54)

Hence, Equation (52), gives an expression for $d_{turn}$, together with (32), (41), and (54), that only depends on the design parameter ${\dot{\chi }}_{T_d}$, and the given fixed parameters $V_T,\Delta {\chi }_T,T_p,p_{cmd}$.

3.1.2. Implementation Considerations

To implement the equations presented in the previous section, we need to start by defining the speed, which is used in the said equations. This is the speed with which the maneuvers are planned. Looking at Equation (34), it is clear that the speed is proportional to the maneuver radius $r_c$. Therefore, the higher the speed, the bigger the maneuver circle will be. This means that, if the turns are planned with a speed lower than the speed of the aircraft, then the turn will not be achievable. On the other hand, if the turn is planned with a higher speed than the aircraft is carrying, the maneuver remains feasible. Therefore it is safer to always plan the maneuvers with a speed that is slightly higher than the predicted speed of the aircraft. This will allow us to implicitly account for wind, which would change the kinematic speed of the aircraft. We can calculate the planning speed as shown in Equation (55) where $V_{T,cmd}$ is the commanded kinematic speed, $V_T$ is the kinematic speed of the aircraft, and $V_{buffer}$ is a the safety buffer speed. The buffer is defined from the knowledge of the aircraft and wind speeds in the operational environment.

$$V_p=max\left(V_{T,cmd},V_T\right)+V_{buffer}$$

(55)

Using the planning speed, the calculation process can be summarized as follows:

1.

Calculate the maneuver circle $r_c$ using Equation (34), where ${\dot{\chi }}_{T_d}$ is defined during design.

2.

Calculate the corresponding bank angle ${\mu }_{cmd}$ using Equation (33).

3.

Calculate the clothoid transition time $t_{cl}$ using Equation (35), where $p_{cmd}$ is the maximum achievable roll rate.

4.

Calculate the clothoid shaping parameter A using Equation (40).

5.

Calculate the clothoid running parameter (dimensionless time) ${\tau }_{cl}$ using Equation (37).

6.

Calculate the clothoid displacement $\left[\Delta x_{cl}\left(\tau \right),\Delta y_{cl}\left(\tau \right)\right]$ using Equation (32).

7.

Calculate the turn distance $d_{turn}$ using Equation (52).

The turn distance is then used to calculate the ${\overrightarrow{r}}^{clp1}$ position (see Figure 7), which defines the point at which the turn mode is activated.

3.2. Vertical Plane

The flyby maneuver in the vertical plane retains the same structure as the horizontal plane. The difference is that the transition segment is made of a higher order polynomial instead of a clothoid section, in order to achieve a higher degree of continuity in the altitude change. These transition segments are used to connect segments of constant altitude with segments of linearly changing altitude as seen in Figure 9 where the linear sections are shown in blue and the polynomial transition section shown in orange. Note that the vertical plane does not consider the aircraft climb performance. It only considers the continuity constraints. It is assumed that the aircraft is able to complete the transition segment within the same duration as the horizontal turn.

3.2.1. Altitude Polynomial Definition

In order to define the polynomial connecting a straight line and circle, we first start by defining the conditions that we want the polynomial to fulfill. The conditions are based on the desired degree of continuity.

We want to have a continuity up to the fourth derivative of the altitude, or $C^4$ continuity. Therefore we have five conditions for the start of the polynomial, and five conditions for the end. Therefore, we need to have a ninth order polynomial, as it will have 10 coefficients, implying 10 degrees-of-freedom to apply constraints:

$$h_T^F\left(s_{hor}\right)=\sum _{m=0}^9{a}_m·s_{hor}^m$$

(56)

where $s_{hor}$ is the current distance traveled in the horizontal plane.

The next thing we need is the climb angle ${\gamma }_T^F$ along the curve. We start by writing down the expression for the climb angle at any point along a line between two altitudes as:

$${\gamma }_T^F{\left(s_{hor}\right)}_{line}=arctan\left(\frac{\Delta h_T^F}{s_{hor}}\right)$$

(57)

where $\Delta h_T^F/s_{hor}$ represents the slope of the line at the $s_{hor}$ point as seen in Figure 9, assuming that the starting point is defined by $s_{hor}=0$. Therefore, for the curve, we need to find the slope and take the arctangent of the later.

The slope of a curve is the first derivative of that curve. Therefore, we take the derivative of Equation (56) with respect to $s_{hor}$:

$${\dot{h}}_T^F\left(s_{hor}\right)=\sum _{m=1}^{9}m·a_m·s_{hor}^{m-1}$$

(58)

We then take the arctan of the slope to find the climb angle ${\gamma }_T^F$:

$${\gamma }_T^F\left(s_{hor}\right)=arctan\left(\sum _{m=1}^{9}m·a_m·s_{hor}^{m-1}\right)$$

(59)

For flight trajectory applications, ref. [15] gives an example table of conditions that need to be satisfied:

In

Table 1. Conditions for altitude polynomial [15].

Condition Parameter	Condition for $\mathit{s}=0$	Condition for $\mathit{s}={\mathit{s}}_{\mathit{tot},\mathit{hor}}$
$h_T^F$	$h_T^{P_1}$	$h_T^{P_2}$
${\left(h_T^F\right)}^{\prime }$	$tan{\gamma }_T^{P_1}$	$tan{\gamma }_T^{P_2}$
${\left(h_T^F\right)}^{\prime \prime }$	0	0
${\left(h_T^F\right)}^{\prime \prime \prime }$	0	0
${\left(h_T^F\right)}^{\prime \prime \prime \prime }$	0	0

, s is the running parameter of the polynomial and $s=0$ represents the starting point of the polynomial ($P_1$). In contrast, $s_{tot,hor}$ is the total distance expected to be traveled and $s=s_{tot,hor}$ represents the end point of the polynomial ($P_2$).

3.2.2. Polynomial Coefficients

Now that we have defined the general shape of the polynomial and the conditions that it needs to satisfy, we need to find the coefficients of the polynomial. We can do this by building a system of 10 equations, such that each equation satisfies one condition in Table 1 as shown in Equation (60).

$$\left\{\begin{array}{cc}{h}_T^F\left(s_{hor}\right)\hfill & ={\left.h_T^{P_1}\right|}_{s_{hor}=0}\hfill \\ h_T^F\left(s_{hor}\right)\hfill & ={\left.h_T^{P_2}\right|}_{s_{hor}=s_{tot,hor}}\hfill \\ h_T^F{\left(s_{hor}\right)}^{\prime }\hfill & ={\left.tan{\gamma }_T^{P_1}\right|}_{s_{hor}=0}\hfill \\ h_T^F{\left(s_{hor}\right)}^{\prime }\hfill & ={\left.tan{\gamma }_T^{P_2}\right|}_{s_{hor}=s_{tot,hor}}\hfill \\ h_T^F{\left(s_{hor}\right)}^{\prime \prime }\hfill & ={\left.0\right|}_{s_{hor}=0}\hfill \\ h_T^F{\left(s_{hor}\right)}^{\prime \prime }\hfill & ={\left.0\right|}_{s_{hor}=s_{tot,hor}}\hfill \\ h_T^F{\left(s_{hor}\right)}^{\prime \prime \prime }\hfill & ={\left.0\right|}_{s_{hor}=0}\hfill \\ h_T^F{\left(s_{hor}\right)}^{\prime \prime \prime }\hfill & ={\left.0\right|}_{s_{hor}=s_{tot,hor}}\hfill \\ h_T^F{\left(s_{hor}\right)}^{\prime \prime \prime \prime }\hfill & ={\left.0\right|}_{s_{hor}=0}\hfill \\ h_T^F{\left(s_{hor}\right)}^{\prime \prime \prime \prime }\hfill & ={\left.0\right|}_{s_{hor}=s_{tot,hor}}\hfill \end{array}\right.$$

(60)

The system of equations (60) is not overly difficult to solve because we are guaranteed that at least half of them will be simplified by the fact that they are evaluated at $s_{hor}=0$. After solving it, we obtain expressions for each of the coefficients, as shown in

Table 2. Coefficients of the altitude polynomial of Equation (56).

Coefficient	Expression
$a_0$	$h_{\mathrm{p}1}$
$a_1$	$h_{\mathrm{p}2}$
$a_2$	0
$a_3$	0
$a_4$	0
$a_5$	$-\frac{14\left(9h_{\mathrm{p}1}-9h_{\mathrm{p}2}+5h_{\mathrm{b}1}{s}_{\mathrm{tot}}+4h_{\mathrm{b}2}{s}_{\mathrm{tot}}\right)}{{s_{\mathrm{tot}}}^5}$
$a_6$	$\frac{28\left(15h_{\mathrm{p}1}-15h_{\mathrm{p}2}+8h_{\mathrm{b}1}{s}_{\mathrm{tot}}+7h_{\mathrm{b}2}{s}_{\mathrm{tot}}\right)}{{s_{\mathrm{tot}}}^6}$
$a_7$	$-\frac{20\left(27h_{\mathrm{p}1}-27h_{\mathrm{p}2}+14h_{\mathrm{b}1}{s}_{\mathrm{tot}}+13h_{\mathrm{b}2}{s}_{\mathrm{tot}}\right)}{{s_{\mathrm{tot}}}^7}$
$a_8$	$\frac{5\left(63h_{\mathrm{p}1}-63h_{\mathrm{p}2}+32h_{\mathrm{b}1}{s}_{\mathrm{tot}}+31h_{\mathrm{b}2}{s}_{\mathrm{tot}}\right)}{{s_{\mathrm{tot}}}^8}$
$a_9$	$-\frac{35\left(2h_{\mathrm{p}1}-2h_{\mathrm{p}2}+h_{\mathrm{b}1}{s}_{\mathrm{tot}}+h_{\mathrm{b}2}{s}_{\mathrm{tot}}\right)}{{s_{\mathrm{tot}}}^9}$

. Note that, for readability, the following simplifications are assumed:

• $h_T^{P_n}\to h_{p_n}$;
• $tan{\gamma }_T^{P_1}\to h_{b_n}$;
• $s_{tot,hor}\to s_{tot}$.

4. Trajectory Construction: Flyby Implications for eVTOLs

4.1. Minimum Distance between Waypoints

The developments in [15] are carried out assuming a fixed-wing aircraft. When considering eVTOLs applications, we need to also consider the small flight distances anticipated. Therefore, we need to calculate the minimum distance between waypoints that the trajectory generation system can handle.

When we look at the geometrical construction of the flyby maneuver in Figure 6, we can empirically establish that we need at least the $2·d_{turn}$ to account for the turn-in section and the turn-out section of the potential flyby we made at the previous waypoint. To explicitly show this, we can consider the flight shown in Figure 10, where we have two flyby transitions at the $TO$ and at the $NEXT$ waypoints.

We see that the distance between the $TO$ and $NEXT$ waypoints needs to necessarily be greater than twice that of the turn distance in case that the heading change and velocity are the same at both waypoints. In practice, the turn distance could be different for the two maneuvers, as it also depends on the angle between the legs ${\alpha }_T$ and the velocity $V_T$ at which the legs are flown.

In Section 3 and Section 3.1.1, we derived an expression for the turn distance for the fly-by which depends only on the given fixed parameters $V_T,\Delta {\chi }_{cl},T_p,p_{cmd}$ and the desired turn rate ${\dot{\chi }}_{T_d}$ (see Equation (52)). The roll time constant $T_p$ and roll rate command $p_{cmd}$ are given by the aircraft performance and inner loop dynamics. Hence, we can summarize: The minimum distance between the waypoints is influenced by

• Aircraft speed $V_T$: The slower the aircraft, the less distance we need to turn.
• Desired turn rate ${\dot{\chi }}_{T_d}$: The higher the turn rate, the less distance we need to turn.

The speed of the aircraft is given in the flight plan for each leg. That leaves the desired turn rate as the only variable that can be changed during design to minimize the minimum distance between waypoints, which is important for eVTOLs applications. We can now advance the following theorem:

Theorem 1.

For trajectory generation systems designed around a fixed desired turn rate, the minimum distance between waypoints can be reduced by increasing the fixed desired turn rate, so long as it stays within the performance range of the aircraft.

4.2. Maximum Angle ${\alpha }_T$ between Legs

Although this may seem trivial because an angle of ${\alpha }_T={180}^{\circ }$ between two legs will result in a straight line flight, we will see in this section that ${\alpha }_T$ is actually limited and not allowed to exceed a certain value. Furthermore, we will see that increasing the turn rate as given by Theorem 1 in order to reduce the minimum turn distance will decrease the maximum angle ${\alpha }_T$ allowed between the legs.

When we take a closer look at the geometry of flyby transitions, we see that the total change in course angle $\Delta {\chi }_T$ is the sum of the change in course incurred while flying the turn-in, circle, and turn-out segments:

$$\Delta {\chi }_T=\Delta {\chi }_T^{*}+2{\phi }_{cl}$$

(61)

where ${\phi }_{cl}$ is simply ${\tau }_{cl}^2$ according to Equation (31) and with ${\tau }_{cl}$ given in Equation (54), ${\phi }_{cl}$ becomes

$${\phi }_{cl}=\frac{{\dot{\chi }}_{T_d}}{2}\left(2T_p+\frac{arctan\left(\frac{V_T{\dot{\chi }}_{T_d}}{g_0}\right)}{p_{cmd}}\right)$$

(62)

As we see, the change in course angle incurred while flying the clothoid segment does not depend on the waypoints or the angle between them. Rather it depends only on the aircraft performance given by $T_p$ and $p_{cmd}$, speed $V_T$, and the desired turn rate ${\dot{\chi }}_{T_d}$, which is the only parameter that can be set during design.

This implies that $\Delta {\chi }_T$ is limited to $\Delta {\chi }_T\ge 2{\phi }_{cl}$ in order for Equation (61) to hold true, and only the desired turn rate can be modified in order to affect this limitation. We know that $\Delta {\chi }_T$ is directly related to the flight plan by $\Delta {\chi }_T=\pi -{\alpha }_T$, where ${\alpha }_T$ is the angle between two legs. This means that the limitation applies directly to the angle between legs in the flight plan that the system is able to handle. We can also see this graphically in Figure 8. The total course angle change $\Delta {\chi }_T$ needs to allow for two transition segments that will have a course angle change of ${\phi }_{cl}$ each, which is independent of the flight plan itself. We can write down the expression for minimum course angle change as:

$$min\left(\Delta {\chi }_T\right)=2·{\phi }_{cl}={\dot{\chi }}_{T_d}\left(2T_p+\frac{arctan\left(\frac{V_T{\dot{\chi }}_{T_d}}{g_0}\right)}{p_{cmd}}\right)$$

(63)

And the maximum angle between two legs in the flight plan is:

$$max\left({\alpha }_T\right)=\pi -min\left(\Delta {\chi }_T\right)$$

(64)

We can now advance the following theorem:

Theorem 2.

For trajectories designed as a combination of straight lines, transition sections, and circular arcs, there is an inherent limitation on the maximum angle between two legs of the flight plan that is independent of the flight plan geometry.

Theorem 2 implies that, if the flight plan has any two legs such that the angle between the legs is a perfect ${180}^{\circ }$ or within some threshold around it, then the two legs will be treated as a straight line. If, however, the angle is in the range $[max\left({\alpha }_T\right),180)$, i.e., above its limit value, then the angle required by the turn-in and turn-out, i.e., $2{\phi }_{cl}$, will be larger than the $\Delta {\chi }_T$ that is specified in the flight plan. This would lead to a maneuver that is not feasible considering Equation (61), as well as Figure 8, and that $\Delta {\chi }_T^{*}\ge 0$. Hence, such a situation needs to be inhibited as it results in a non-deterministic behavior. This can be achieved by making the desired turn rate dynamic, depending on the required ${\alpha }_T$.

4.2.1. Making the Turn Rate Dynamic

Theorem 1 shows that, in order to decrease the minimum distance between waypoints, which is important for eVTOL aircraft, we need to increase the desired turn rate. Hence, we need to choose the desired turn rate such that it allows us to meet the required minimum distances between waypoints, so long as it is within the performance range of the aircraft.

When we look at Theorem 2 and Equation (63), we notice that there is a limitation on the maximum angle between legs and it not only depends on the performance of the aircraft, but also on the desired turn rate. The higher the turn rate, the lower the maximum angle allowed between two legs. Equation (63) shows that, in order to increase the allowed ${\alpha }_{T,max}$ we need to decrease the desired turn rate. Hence, the initial choice of the desired turn rate, which fulfills the minimum distance requirement assumed during planning, would lead to overly restrictive allowed values for ${\alpha }_T$.

We can overcome this apparent contradiction by making the desired turn rate dynamic, i.e., in the case that the desired ${\alpha }_T$ given by the planned waypoints, exceeds ${\alpha }_{T,max}$ the turn rate will be decreased such that the resulting ${\alpha }_{T,max}$ allows for the desired ${\alpha }_T$. Of course this will again increase the turn distance $d_{turn}$ required to accomplish the maneuver, and it needs to be confirmed that it does not exceed the available distance between the two waypoints. Section 4.2.2 will address this topic. For now, in order to calculate the required turn rate to obtain the desired ${\alpha }_{T,max}$ value, Equation (63) needs to be solved for ${\dot{\chi }}_{T_d}$. Looking at the equation, it is apparent that it cannot be solved analytically, because the variable of interest, i.e., ${\dot{\chi }}_{T_d}$, is both inside and outside an arctan function. This makes the equation transcendental. We therefore need to find an approximation for the arctangent function, which allows us to solve (63) for ${\dot{\chi }}_{T_d}$.

There are several known approximations for the arctangent function with different levels of accuracy. The first step in choosing an approximation is to define the magnitude of the arctan argument. And this is application dependent. Figure 11 shows the magnitude of the arctan function argument, i.e., $V_T{\dot{\chi }}_{T_d}/g_0$, for the typical range of motion of eVTOL aircraft, where the maximum velocity is 30 m/s and the maximum turn rate is 12 $\mathrm{d}\mathrm{e}\mathrm{g}/\mathrm{s}$ based on experimental data. As we can see, the argument remains strictly positive, and reaches a maximum of $0.8$.

We also have to keep in mind that Equation (63) has a product between the ${\dot{\chi }}_{T_d}$ outside the arctan function and inside. Therefore, the equation that needs to be solved after substituting the arctan will have an order $n+1$, where n is the order of the approximation. This makes the Taylor series approximations of order $n>2$ unlikely candidates as the resulting equation becomes too complex.

Figure 11. Typical arctan argument for eVTOL aircraft (see Equation (63)).

Figure 11. Typical arctan argument for eVTOL aircraft (see Equation (63)).

Another promising arctan approximation is:

$$t\left(x\right)=\frac{x}{1+b{x}^2}$$

(65)

where b is a tuning parameter. This function is a particularly good candidate because it has a low order and has the same shape as the arctan function. It is, however, overly complicated for our use case, as the argument of the function is strictly positive and no larger than $0.8$, and the approximation leads to a third-order equation whose solution is not trivial. Therefore, we instead choose to approximate it with a line fitted to minimize the mean squared error. Figure 12 shows a comparison between the actual function and the approximation, where $b_0=0.89813$ is the tuning parameter.

Substituting the approximation into Equation (63), and replacing $max\left({\alpha }_T\right)$ with the desired angle ${\alpha }_{T,d}$, we can solve for ${\dot{\chi }}_d$:

$${\dot{\chi }}_{T_d}=\left(\begin{array}{c}\frac{\sqrt{g_0{p}_{cmd}\left(g_0{p}_{cmd}{T_p}^2+\left(\pi -{\alpha }_{T,d}\right)V_T{b}_o\right)}-T_p{g}_0{p}_{cmd}}{V_T{b}_o}\\ -\frac{\sqrt{g_0{p}_{cmd}\left(g_0{p}_{cmd}{T_p}^2+\left(\pi -{\alpha }_{T,d}\right)V_T{b}_o\right)}+T_p{g}_0{p}_{cmd}}{V_T{b}_o}\end{array}\right)$$

(66)

The ${\dot{\chi }}_d$ computed here is used for the planning of the maneuver as shown in Section 3, and is defined to be strictly positive. Therefore, we take the solution which is positive:

$${\dot{\chi }}_{T_d}=\frac{\sqrt{g_0{p}_{cmd}\left(g_0{p}_{cmd}{T_p}^2+\left(\pi -{\alpha }_{T,d}\right)V_T{b}_o\right)}-T_p{g}_0{p}_{cmd}}{V_T{b}_o}$$

(67)

Proof. 

Validity of Equation (67): Since the denominator $V_T{b}_o$ is positive, we can show that (67) is positive by showing that

$$\sqrt{g_0{p}_{cmd}\left(g_0{p}_{cmd}{T_p}^2+\left(\pi -{\alpha }_{T,d}\right)V_T{b}_o\right)}>T_p{g}_0{p}_{cmd}$$

(68)

Because both terms under the square root are greater than zero: $T_p^2{g}_0^2{p}_{cmd}^2>0$ as well as $g_0{p}_{cmd}\left(\pi -{\alpha }_{T,d}\right)V_T{b}_0>0$, the square root of their sum will also be greater than zero: $\sqrt{T_p^2{g}_0^2{p}_{cmd}^2+g_0{p}_{cmd}\left(\pi -{\alpha }_{T,d}\right)V_T{b}_0}>0$. Because, additionally, $T_p{g}_0{p}_{cmd}>0$ holds, the inequality in (68) can be proven to hold true by showing that the inequality holds for the squared values of both sides:

$$T_p^2{g}_0^2{p}_{cmd}^2+g_0{p}_{cmd}\left(\pi -{\alpha }_{T,d}\right)V_T{b}_0>{T_p}^2{g_0}^2{p}_{cmd}^2$$

(69)

which is the case because the left part of the inequality is a sum of the right part and some other positive expression. Therefore, the left part is always greater than the right part, making the expression under (67) strictly positive.

Another thing to consider is that Equation (67) gives the expression for the maximum ${\dot{\chi }}_d$ that we can have in order to be able to perform a flyby with a specified or desired angle ${\alpha }_{T,d}$ between the legs. We can take a value that is $10\%$ lower as a safety margin. We then obtain the final expression for ${\dot{\chi }}_d$:

$${\dot{\chi }}_{T_d}=0.9\frac{\sqrt{g_0{p}_{cmd}\left(g_0{p}_{cmd}{T_p}^2+\left(\pi -{\alpha }_{T,d}\right)V_T{b}_o\right)}-T_p{g}_0{p}_{cmd}}{V_T{b}_o}$$

(70)

□

4.2.2. Integrating the Dynamic Turn Rate

Now that we have derived the expression for the desired turn rate as a function of the desired angle between legs, we need to integrate it into the system as a whole, while still applying Theorem 1, i.e., making sure that the minimum distance required for the flyby, which results from the applied turn rate, does not exceed the available distance between the planned waypoints. We can perform this as shown in Figure 13. Note that the algorithm in Figure 13 only terminates when the trajectory generation module terminates. Therefore, there is no termination link depicted.

When there is an upcoming flyby maneuver, we check whether the design-fixed turn rate results in a $max\left({\alpha }_T\right)$ that is large enough to encompass the required ${\alpha }_{T,d}$ of the flyby. If this is the case, then we continue with the maneuver planning.

If, however, $max\left({\alpha }_T\right)<{\alpha }_{T,d}$, then we use Equation (70) to calculate a new ${\dot{\chi }}_{T_d}$, which is then used only for the upcoming maneuver. As mentioned earlier, decreasing the desired turn rate increases the turn distance needed for the maneuver. This is an inherent limitation of the trajectory-generation concept, which cannot handle arbitrarily small distances between waypoints. We therefore need to ensure that the maneuver remains feasible when we decrease the turn rate. This is performed in the mission management when the flight plan is uploaded. The mission management, including the feasibility checks, is described in detail in [28].

Figure 13. Integrating dynamic turn rate.

Figure 13. Integrating dynamic turn rate.

In this section we conducted an analysis of the flyby geometry and exposed some limitations to the construction of trajectories as a combination of the straight lines, transition segments, and circular arcs, with a fixed design turn rate. For eVTOL aircraft, the applications typically require short distances between waypoints. We showed that this distance is limited by the distance required to make turns. One way to decrease the turn distance is to increase the turn rate. For fixed-wing aircraft, it is general practice to use a standard turn rate of $3\mathrm{d}\mathrm{e}\mathrm{g}/\mathrm{s}$. However, eVTOL aircraft can achieve much higher turn rates. Therefore we recommended using the highest turn rate achievable by the aircraft. Another thing that the analysis revealed is that the higher turn rate will introduce a limitation to the maximum angle allowed between two legs in the flight plan. We suggested making the turn rate dynamic in order to alleviate this limitation. A more detailed discussion of this analysis and implications is presented in the next section.

5. Consolidation and Results

In Section 2 we introduced the preliminaries of trajectory generation, which can be summarized as follows:

Following the proof in [14], it is common practice to build up trajectories for flying vehicles moving by combining straight lines and circular arcs, to make the so-called Dubins paths [29,30]. One problem that arises with Dubins paths is that connecting straight lines and circular arcs lead to a step in the curvature of the trajectory. As shown in Section 3.1, the step in curvature can be explained physically as an instantaneous step in the required kinematic bank angle of the aircraft. One solution to the curvature step problem is to use transition segments. These segments are typically high-order polynomials such as Bezier curves [31], or other types of cubic splines [32]. In this paper we used the clothoid approach presented in [11].

Regardless of the type of transition segment used, it is needed to accomplish one task; that is, to build up the bank angle of the aircraft. Equation (35) shows an approximation of the time it would take an aircraft to build up a certain bank angle, given the performance of the aircraft, as reflected through the roll time constant ($T_p$) and the roll rate command ($p_{cmd}$). The roll rate command described here is usually the maximum roll rate that the aircraft can achieve. The geometrical length of the transition segment, which is the clothoid in our case, is such that the initial point has the curvature of 0, and the end point has the curvature of the maneuver circle whose radius is calculated with Equation (34).

In Section 3.1 we go through the derivation of the connection between the time it takes the aircraft to bank, and the clothoid shaping parameter A. This results in Equation (41). And, as shown in Figure 5, the shaping parameter A controls the rate of growth of the clothoid. This means that we have established a direct relationship between the aircraft performance and the physical space occupied by the transition segment. This relationship is important because the implications apply regardless of the type of transition segment used. One of those implications is the turn distance as introduced in Section 3.1.1.

5.1. Turn Distance

The expression of Equation (50) shows that, in addition to the angle between the legs, the turn distance is also, and mainly, affected by the speed, V, turn rate, ${\dot{\chi }}_{T_d}$ and $\Delta x_{cl,f}$, and $\Delta y_{cl,f}$, which indirectly represent the aircraft’s performance as they are the distances occupied by the transition segment.

The turn distance is important because it controls how far apart two successive waypoints need to be. Because, as shown in Section 4.1, we need to account for the turn distance at the previous waypoint, and the turn distance at the upcoming waypoint. This defines the minimum distance that can be between waypoints. This is particularly important for eVTOL aircraft as they are typically intended to fly short distances. And what we see here is that there is an inherently lower limit on the distance that can be between waypoints. And this minimum distance is mainly dependent on the performance of the aircraft.

To see this implication in detail, let us take an example of a hypothetical aircraft configuration defined by:

• Roll rate command: $p_{cmd}={30}^{\circ }$
• Roll time constant: $T_p=0.5$ s
• Design turn rate: ${\dot{\chi }}_{T_d}={10}^{\circ }/\mathrm{s}$

The value of the design turn rate ${\dot{\chi }}_{T_d}$ is chosen to be close to the maximum achievable by the aircraft following Theorem 1 which says that the turn distance $d_{turn}$ can be decreased by increasing the design turn rate.

We can start by looking at how $d_{turn}$ varies with the angle between the legs and the velocity.

In Figure 14 we can see that there is a linear relationship between $d_{turn}$ and the velocity V. The higher the velocity, the more distance we need to turn, and by implication, the more distance we need between waypoints. The same Figure also shows a nonlinear dependency of $d_{turn}$ on the angle between legs ${\alpha }_T$, shown by the different color. This is better shown in Figure 15.

Note that the angle ${\alpha }_T$ was chosen to be only up to ${177}^{\circ }$. This is because there is a margin of $3^{\circ }$ around the ${180}^{\circ }$ line to account for small errors when building a flight plan and placing waypoints in a straight line.

What we see in Figure 15 is that, for values of ${\alpha }_T$ lower than a certain threshold, which can be put around 30 to 40 degrees, there is a a slightly sharp increase in the $d_{turn}$. However, as ${\alpha }_T$ start increasing, the turn distance stays almost constant, reflecting the same thing seen in Figure 14. Therefore, in the implementation of the trajectory generation system, the angle between legs of the flight plan for a flyby is limited to ${30}^{\circ }$ to avoid the region of large turn distances.

Another thing that we notice both in Figure 14 and Figure 15 is that the turn distance is heavily influenced by the velocity. So, a possible way to decrease $d_{turn}$ would be to decrease the velocity. But the velocity is given as part of the flight plan as an attribute of each leg. Although this possibility is certainly feasible, it begs several questions, such as, at what point does the aircraft need to start slowing down? If the velocity is to be reduced right before the turn, then the problem becomes vaguely defined, because, as we see in Equation (50), the turn distance does depend on the velocity as well. The turn will start when the aircraft reaches the ${\overrightarrow{r}}^{clp1}$ point from Figure 7. This point is calculated as the point at $d_{turn}$ distance from the $TO$ waypoint. Therefore, the velocity is very much needed in the planning of maneuvers, which happens before the maneuver occurs.

Figure 14. Influence of V and ${\alpha }_T$ on $d_{turn}$.

Figure 14. Influence of V and ${\alpha }_T$ on $d_{turn}$.

Another problem with reducing the velocity is the energy loss that will occur while decelerating and accelerating the aircraft again. Especially if the reason driving the need to reduce the speed is to reduce the turn distance, which will reduce the minimum distance between waypoints. It certainly means that there is a short distance between waypoints, and we might need to start decelerating again as soon as we finish the acceleration at the previous waypoint. We then end up in a vicious circle of accelerations and decelerations, which also defeats the whole purpose of having a velocity assigned to the leg in the flight plan.

It is worth noting that, although the argument above seems to show all the reasons not to try to change the speed assigned in the flight plan in an online trajectory generation system, it is certainly not intended to deter the reader from considering this approach. It is an approach that might work, if considered thoroughly, in addition to a consideration of the regulations over speed changes in the type of airspace where the system is intended to operate.

The solution presented in this paper, to reduce turn distances, and by implication the minimum distance between waypoints, is to simplify increasing the turn rate, as reflected in Theorem 1. The influence of the turn rate is shown in Figure 16 where the black arrow show indicates the increasing velocity.

As we see in Figure 16, the turn distance can be decreased considerably by increasing the turn rate ${\dot{\chi }}_{T_d}$. But this effect is only considerable up to a certain value where the influence on $d_{turn}$ becomes constant. Additionally, in the same Figure, we see the inconsiderable effect of the angle between the legs in the flight plan, and we see the $d_{turn}$ lines being shifted upwards as the velocity increases.

Here, it is also worth mentioning that Theorem 1 is not meant to mean that we need to go to the maximum achievable turn rate. We need to consider Equation (33) which shows that, for different velocities, the turn rate can result in different bank angles. We therefore need to consider the maximum allowed bank angle for the specific type of aircraft and type of intended operations when deciding on the the design turn rate to set in the trajectory generation system.

Figure 15. Influence of ${\alpha }_T$ and V on $d_{turn}$.

Figure 15. Influence of ${\alpha }_T$ and V on $d_{turn}$.

Another implication of the flyby geometrical construction is the maximum angle allowed between legs ${\alpha }_T$, as presented in Section 4.2.

Figure 16. Influence of ${\dot{\chi }}_{T_d}$, ${\alpha }_T$ and V on $d_{turn}$.

Figure 16. Influence of ${\dot{\chi }}_{T_d}$, ${\alpha }_T$ and V on $d_{turn}$.

5.2. Maximum Angle between Legs

As shown in Section 4.2, when we go away from Dubins paths and add transition segments to connect straight lines and circular arcs, there is a limitation that gets introduced to the angle between the legs of the flight that the system is able to handle with deterministic behavior. This limitation comes directly from the geometrical construction of the maneuver, and is influenced directly by the performance of the aircraft. As we see in Figure 8, there is a course angle change ${\varphi }_{cl}$ incurred while flying the turn-in transition segment, and another ${\varphi }_{cl}$ during the turn-out.

As we established earlier, the transition segment depends on the performance of the aircraft. Therefore, it is evident that the course angle incurred during the transition segment also depends on the aircraft performance. This is reflected in Equation (62), where we only see the turn rate, the roll time constant, the roll rate command, and the velocity. This means that, for the typical trajectory generation systems that are designed with a fixed turn rate ${\dot{\chi }}_{T_d}$, the ${\varphi }_{cl}$ angle will always be the same for a given velocity, regardless of the total course angle change $\Delta \chi$ required at one waypoint or another. The implication is that Equation (61) needs to always hold true, such that we are always able to make a turn-in, reach a nominal bank angle, and then a turn-out. For the user of the system, it means that there is an inherent maximum allowed angle ${\alpha }_T$ between two legs in the flight plan for the maneuver to be feasible.

Unsurprisingly, the solution to mitigate this limitation is to reduce the turn rate as shown in Section 4.2.1. This solution works because, as Equation (33) shows, the turn rate defines the bank angle. And we are simply saying that we do not always want to go to the maximum bank angle (defined by the design turn rate). Sometimes, when the angle between the legs is large, i.e., the change in course that we need to perform is small, we can carry out the maneuver with a smaller bank angle, and by implication, a smaller turn rate. By reducing the turn rate, the maneuver circle becomes larger (see Equation (34)). As the circle becomes larger, the curvature becomes smaller (see Equation (19)). To satisfy Equation (38) with the new curvature, the arc length of the clothoid decreases, which also, via Equation (31), decreases the course angle change incurred during the transition segment (This can also be seen graphically in Figure 8). And this makes Equation (61) valid again, making the maneuver feasible.

The expression for the turn rate as a function of the angle ${\alpha }_T$ is shown in Equation (70). It is worth noting that this expression is not universally applicable. As shown in Section 4.2.1, Equation (70) is obtained by solving Equation (63) for ${\dot{\chi }}_{T_d}$. To reach the solution, the arctan function was approximated with a straight line. This approximation was accepted as being valid because, for the use case in this paper, the argument of the arctan function, i.e., $V_T{\dot{\chi }}_{T_d}/g_0$, reaches a maximum of $0.8$ as shown in Figure 11. Therefore, the straight line approximation was valid (see (Figure 12)). For another application, this approximation might need to be replaced with one that is more accurate, and the expression in Equation (70) would be invalid. It would then need to be obtained again by solving Equation (63) with the new arctan approximation.

To see the effects of the other variables on the maximum allowed angle between legs, we can start by looking at Figure 17.

We see in Figure 17 that, as we decrease ${\dot{\chi }}_{T_d}$, the lines of $max\left({\alpha }_T\right)$ start converging towards ${180}^{\circ }$ for all the velocities shown in Section 4.2.1. We can see this clearly also in Figure 18, where $max\left({\alpha }_T\right)$ becomes practically constant, around ${180}^{\circ }$, over the change in velocity, for lower turn rates.

Hence, if the ${\alpha }_T$ of the two given legs of the flight plan is higher than our $max\left({\alpha }_T\right)$ with the design ${\dot{\chi }}_{T_d}$, then we can decrease ${\dot{\chi }}_{T_d}$, such that the resulting $max\left({\alpha }_T\right)$ is larger or equal to the desired ${\alpha }_T$. This is shown in Figure 13.

Although this seems to solve all our problems, there is a catch. Looking at Figure 16, we see that decreasing the turn rate results in a higher turn distance $d_{turn}$. And we run the risk of requesting a ${\dot{\chi }}_{T_d}$ which results in a $d_{turn}$ higher than the distance between the waypoints. For the category of trajectory generation systems considered here there is just no way around this apparent contradiction. Therefore, it is mandatory to have a trajectory verification module that includes all the calculations presented in this paper among other things, in order to check the feasibility of the flight plan when uploaded. The flight plan will be sent to the trajectory generation system only if it is feasible [28].

Figure 17. Influence of ${\dot{\chi }}_{T_d}$, and V on $max\left({\alpha }_T\right)$.

Figure 17. Influence of ${\dot{\chi }}_{T_d}$, and V on $max\left({\alpha }_T\right)$.

Figure 18. Influence of V and ${\dot{\chi }}_{T_d}$ on $max\left({\alpha }_T\right)$.

Figure 18. Influence of V and ${\dot{\chi }}_{T_d}$ on $max\left({\alpha }_T\right)$.

5.3. Simulation Results

5.3.1. Testing Environment Overview

With the equations from Section 3, an online trajectory generation system was implemented in Simulink^® following strict guidelines, allowing us to obtain a code-compliant software according to the DO-178/DO-331 [33]. The clothoid-based lateral implementation of the trajectory generation system has been successfully flight tested on a fixed wing aircraft [13].

As input, the system takes a flight plan consisting of waypoints defining the route to be followed, and generates commands that are sent to a trajectory control system. Each waypoint in the flight plan is defined by its location (longitude, latitude, and altitude), and the speed that the aircraft needs to reach during the leg leading to the said waypoint. For the results shown in this section, a point mass model was implemented to simulate the aircraft and inner loop dynamics, in order to test the functionalities of the trajectory generation system in a neutral environment without the overhead of multiple other systems. An overview of the test harness is shown in Figure 19.

As seen in Figure 19, the flight-plan-handling module is responsible for the selection of the specific waypoint in the flight plan to be sent to the flight path construction module (FPC). The FPC module in turn contains a waypoint buffer which always contains the three waypoints necessary for the maneuver construction as discussed in Section 3. After all the geometry is calculated, the commands to be sent to the trajectory controller are computed. This includes the course angle ${\chi }_T$ and its derivatives, the climb angle ${\gamma }_T$ and its derivatives, and the error dynamics, which are obtained as the difference between the position of the aircraft and the position of the trajectory reference point. The calculation of the trajectory reference point is detailed in [15]. The calculations for the relative kinematics and derivatives of the commands are covered in [12,34].

Figure 19. Trajectory generation testing module implementation overview.

Figure 19. Trajectory generation testing module implementation overview.

5.3.2. Point Mass Model

The first step in the definition of the point mass model is to define the kinematic speed $V_T$ of the aircraft. This is feedforwarded directly from the attributes of the waypoints by the trajectory generation module (TrajGen). The kinematic speed is then fed through a first order filter that simulates the velocity dynamics, and allows to obtain the derivative of the velocity as shown in Figure 20.

Now, we can use the ${\chi }_T$ and ${\gamma }_T$ commands from the TrajGen module to transform the kinematic speed into its North East Down (NED) component as shown in Equation (71).

$$\overrightarrow{V_O}=\left[\begin{array}{c}{u}_T\\ v_T\\ w_T\end{array}\right]=V_T·\left[\begin{array}{c}cos{\chi }_T·cos{\gamma }_T\\ sin{\chi }_T·cos{\gamma }_T\\ -sin{\gamma }_T\end{array}\right]$$

(71)

The NED velocities represent a change in position, which can then be propagated into a change in the WGS84 coordinates (longitude, latitude, and height). To achieve that, we start by calculating the different radii of curvature of the WGS84 spheroid as shown in Figure 21 where the different line formats and colors are intended to improve readability and highlight the different variables in the figure.

• The prime vertical radius of curvature, $N_{\phi }$ is:

  $$N_{\phi }=\frac{a}{1-e^2{sin}^2\phi }$$

  (72)
• The meridian radius of curvature, $M_{\phi }$ is:

  $$M_{\phi }=N_{\phi }·\frac{1-e^2}{1-e^{2}si{n}^2\phi }$$

  (73)

  where:

• $\phi$ is the latitude;
• $a=6,378,137$m is the length of the semi-major axis;
• $b=6,356,752.3142$m is the length of the semi-minor axis;
• e is the first eccentricity, calculated as:

  $$e^2=2f-f^2$$

  (74)

  where the flattening f is:

  $$f=\frac{a-b}{a}$$

  (75)

The derivation of the geodetic equations above is given in [35].

Using the NED velocity components and the radii of the earth, we find the first derivative of the geodetic coordinates as:

$$\left[\begin{array}{c}\dot{\lambda }\\ \dot{\phi }\\ \dot{h}\end{array}\right]=\left[\begin{array}{c}\frac{v_T}{\left(N_{\phi }+h\right)cos\phi }\\ \frac{u_T}{M_{\phi }+h}\\ -w_T\end{array}\right]$$

(76)

where h is the altitude.

We can then integrate the derivatives of Equation (76) starting from an initial position, to obtain the current position. Note that the integrators are reset, frozen, or flipped accordingly, to account for the spheroid shape of the earth.

Figure 21. WGS84 spheroid properties of a point above the surface.

Figure 21. WGS84 spheroid properties of a point above the surface.

5.3.3. Results

For the experiments, a simple flight plan containing three waypoints was defined, similar to Figure 6.

In Section 2 we discussed the generalities of the trajectory generation, and we saw that steps in the curvature of the trajectory occur when using the traditional Dubin’s paths. And accordingly, we discussed the use of the clothoid segment to mitigate this problem. To illustrate this one more time, we simulated the flyby with two TrajGen systems. One based on Dubin’s paths, and another based on the clothoid transition.

Figure 22a shows the turn rate command generated by the TrajGen system during a flyby. For the Dubin’s path-based system (no clothoid), we see an instanteneous change in the turn rate. On the other hand, when using the clothoid transition, we obtain a linearly changing turn rate, representing the linearly changing curvature of the trajectory. Another observation that we can make is that it takes significantly less time for the Dubin’s path system to finish the turn. However, this turn remains non-physical as an aircraft cannot have a discontinuous bank angle behavior, which would lead to unavoidable errors in the tracking behavior of the aircraft. The discontinuity in the turn rate is seen clearly in Figure 22b, where we see that the derivative of the turn rate for Dubin’s paths is infinite.

Further, in Section 5.1 we show the relationship between $d_{turn}$ and $V_T$ showing that an increase in velocity leads to an increase in turn distance. This is well illustrated by the equidistant lines in Figure 23 which move further away from the intersection between the legs as the velocity increases. But, as stated earlier, in this work we did not look at the possibility to adjust the velocity in order to obtain a desired turn distance. Although this is certainly an alternative to the method presented in this paper, we made the basic assumption that the velocity is part of the flight plan uploaded by the user and therefore cannot be changed. Another approach that will be explored in future work is to present the user with an alternative velocity at the moment of the flight plan upload, in case the distance between waypoints is not enough to make turns given the current velocity attributes of the waypoints.

Figure 22. Trajectory-generation commands during flyby.

Figure 22. Trajectory-generation commands during flyby.

One solution presented in the paper to deal with limited distances between waypoints is to simply increase the desired turn rate, given that eVTOL aircraft can achieve much higher turn rates that fixed wing aircraft.

Figure 23. Influence of the velocity on the flight path.

Figure 23. Influence of the velocity on the flight path.

Figure 24 gives us a visual representation of the influence of the turn rate on the turn distance as shown in Figure 16. We can increase the desired turn rate to initially drastically reduce $d_{turn}$ but the effectiveness of increasing $\dot{\chi }$ decreases simultaneously. In the limit, the increase in turn rate has no effect on the decrease in turn distance anymore. This can be seen by the decrease in distance between the flight paths in Figure 24.

Finally, the simulation results will help us build a more visual representation of the results shown at the beginning of Section 5. In Section 5.1, the turn distance $d_{turn}$ was said to be influenced by the angle between the legs ${\alpha }_T$, the velocity $V_T$, and the desired turn rate ${\dot{\chi }}_{T_d}$.

In Figure 25, we can clearly see the influence of the angle between the legs in the flight plan. The smaller ${\alpha }_T$ is, the earlier the turn needs to start, and this results in a higher $d_{turn}$, which in turn would limit the minimum distance between the two waypoints.

All the results shown so far have been looking at the trajectory generation in the horizontal plane. The trajectory in the vertical plane is implemented using the equations shown in Section 3.2.

In Figure 26 we see a straight line in the horizontal plane, and only a climb is performed. This is an example of a case where the waypoints are placed in a straight line, but at different altitudes. At the intersections of the horizontal straight line, and the vertical straight line, a vertical flyby is performed, as seen in Figure 27, in order to maintain a smooth trajectory. The vertical flyby is achieved by the polynomial shown in Section 3.2.

Although the horizontal and vertical planes are developed separately, they are combined in order to be able to generate 3D trajectories. This is illustrated in Figure 28 where the waypoints, indicated in red, are not placed in a straight line, and are placed at different altitudes. The resulting trajectory is a combination of a horizontal flyby with a vertical one, and flying straight lines in both the horizontal and vertical planes in between waypoints. This means that, when two waypoints are placed at different altitudes, they will be connected by a line of constant climb as seen in Figure 26. The feasibility of the climb is not considered in the trajectory generation module. The feasibility gateway [28] verifies the altitudes of the waypoints when they are uploaded, to make sure the climb performance of the aircraft allows us to achieve required climb rate.

These 3D trajectory capabilities are not limited to the simple maneuvers shown so far. The waypoints can be placed at strategic locations in order to achieve much more complicated trajectories as seen in Figure 29. This trajectory is made up of several waypoints with the flyby transition. Circular trajectories can be generated by simply placing flyby waypoints close to each other.

The “start” and “finish” points in Figure 29 do not represent the take-off and landing. As mentioned before, the derivatives of the commands are calculated by making a prediction of the aircraft’s position and trajectory in the future. This requires the use of course and climb angles. These angles are not defined when the aircraft is at a stop. Therefore, the TrajGen module presented in this work can be activated only when the aircraft is already in motion. The takeoff and initial acceleration, and landing can be performed by another module such as the one presented in [36,37]. For the results shown in Figure 29, the trajectory was initialized with a nonzero speed. In addition, it is worth noting that the trajectory was generated for a more agile aircraft, represented by the inner loop dynamics with a roll time constant $T_p=0.6\mathrm{s}$ and a maximum achievable roll rate of $p_{cmd}=50\mathrm{d}\mathrm{e}\mathrm{g}/\mathrm{s}$.

As mentioned earlier, the trajectory is represented by the commands that are generated and sent to the trajectory controller. For the complex 3D trajectory showcased, the commands are shown in Figure 30 and Figure 31.

Looking at the turn rate, ${\dot{\chi }}_{T,cmd}$ in Figure 30a, the first thing we notice is that we have a small region where it is positive, followed by a larger region where it goes into the negative values. This simply represents the direction of the turn as we see in Figure 29, where the first circular maneuver is going to the right, while the following three circular maneuvers are going to the left. The next thing we see is that we do not have a constant turn rate. Rather, it goes to some steady state value and rapidly comes back to 0. This is also reflected in Figure 30a, where the course ${\chi }_{T,cmd}$ exhibits a smooth staircase-like behavior when rising or falling. This is because, as mentioned earlier, the trajectory is made up of several waypoints with the flyby transition. And a flyby is made of a turn-in section, circular arc, and turn-out section, to go back to a straight line. By placing several flyby waypoints close to each other, we are able to simulate circular maneuvers within a certain degree of accuracy. Additionally, we see four apparent jumps in ${\chi }_{T,cmd}$. These happen during the turns, and are simply the angle wrapping between 0 and 360 degrees.

Figure 29. 3D trajectory.

Figure 29. 3D trajectory.

Figure 30. Horizontal plane commands during maneuver of Figure 29.

Figure 30. Horizontal plane commands during maneuver of Figure 29.

Figure 31. Vertical plane commands during maneuver of Figure 29.

Figure 31. Vertical plane commands during maneuver of Figure 29.

Finally, we notice that, for most of the turns, the turn rate ${\dot{\chi }}_{T,cmd}$ goes to a steady state value of just a little over $10$ deg/s. This is the design turn rate. However, for some of the turns, ${\dot{\chi }}_{T,cmd}$ goes to a value that seems to be random. These are turns where the angle between the legs is larger than the limit (see Equation (64)), and the dynamic turn rate formula shown in Section 4.2.1 is triggered, effectively adjusting the desired ${\dot{\chi }}_{T,d}$.

A closer look at the internal moding is shown in Figure 32. In the first sub-figure, we see the moment when the system realizes that ${\dot{\chi }}_{T,d}$ is too high to perform the upcoming turn, and it raises the reduce ${\dot{\chi }}_{T,d}$ flag. We also know that the system is unable to perform the turn because the turn distance $d_{turn}$ in the second sub-figure is zero, meaning that the turn mode will not be activated. We will simply have a jump in the ${\chi }_{T,cmd}$ command when we reach the waypoint. Next, we see in the third sub-figure that the desired turn rate is reduced from $0.2$ to $0.1$, which effectively makes the turn possible, and the $d_{turn}$ goes to a non-zero value.

In Figure 31, we see the commands in the vertical plane, which are given simultaneously with those in the horizontal plane to achieve the 3D trajectory in Figure 29. The values of ${\gamma }_{T,cmd}$ and ${\dot{\gamma }}_{T,cmd}$ are relatively small compared to the other commands because the maneuver is performed with a small change in altitude, with a maximum altitude of 20 m.

Figure 32. Dynamic turn rate moding.

Figure 32. Dynamic turn rate moding.

6. Conclusions

This paper presents a deterministic trajectory-generation module based on clothoid in the horizontal plane and a high-order polynomial in the vertical plane for eVTOL applications. The clothoid-based approach used in the horizontal plane was also discussed in [38], where three suitable functions for generating deterministic trajectories are presented: cubic splines, trigonometric splines, and clothoids. These functions are widely used for path planning due to their simplicity in constructing smooth trajectories.

The trajectory is constructed by combining straight line, clothoid, and circular arc segments. The clothoid is used here as a transition segment between the straight line and clothoid, to provide continuous curvature in the trajectory. We conducted an analysis of the geometrical construction of flyby maneuvers, where a turn is made by a turn-in segment (clothoid), a circular arc segment, and a turn-out segment (clothoid). By using this three-segment approach for the flyby, we show that there will be an inherent limitation to the maximum angle allowed between legs in the flight plan. We show that this limitation is also linked to the use of a constant turn design turn rate used for the maneuvers. Therefore, the solution presented is to make the turn rate dynamic in order to accommodate larger angles between legs.

Our analysis also looks at the turn distance, which affects the required minimum distance between any two waypoints in the flight plan. We show the connection between the turn rate and the turn distance and propose to use a high turn rate for eVTOL in order to minimize the required minimum distance between waypoints. Therefore, the design turn rate is set to the maximum achievable by the aircraft, and is dynamically changed if the angle between the legs is too large. This highlights the need for a feasibility gateway, which checks the flight plan at the moment of upload, to make sure that the distance between waypoint is enough for the maneuvers, and if the angle between the legs is too large, that the resulting turn rate, and correspondingly turn distance would still be feasible.

The analysis presented in this paper is carried out through an analytical derivation of the equations governing the trajectory generation construct, and is confirmed via the simulation of a point mass model and the derived equations.

In addition to the clothoid-based horizontal plane trajectory, we present a polynomial-based approach for the vertical plane. We use a ninth-order polynomial where the expressions of the coefficients of the polynomial are calculated offline by considering the continuity requirements in the trajectory. This allows us to use a high-order polynomial without the high computational requirements typically associated with the optimization of the polynomial. We then combined the horizontal and vertical planes in order to generate 3D trajectories.

The concepts presented in the paper were implemented in Simulink^® R2021a, together with a point mass model assuming perfect tracking of the commands generated by the TrajGen module. The results of our analysis are shown in the paper, along with simulation results.

7. Recommendations for Future Work

The clothoid segment used in this work allowed us to only achieve a second-degree continuity in the control commands. In the future, this function should be replaced with a higher degree of continuity without requiring an excessive amount of computational power.

In addition, as shown in the paper, the clothoid only uses the linearly estimated time for the duration of a bank angle maneuver. This is an oversimplification, and it does not account for the aircraft dynamics. Future works should consider possible ways of accounting for all the dynamics. This point also stands for the polynomial used in the vertical plane.

Another thing is that this paper only considers the flyby maneuver. However, one of the goals of this endeavor is to fully support the maneuvers defined in the ARINC 424 standard [39]. Therefore, future works are planned to extend the supported maneuvers of the system.

Finally, this paper only considers the cruise flight mode for eVTOLs. Future works should consider the extended maneuver envelope of eVTOLs in order to fully explore their capabilities.