Developing a Generalized Trajectory Modeling

Framework for Small UAS Performance in the

Presence of Wind

Sarah N. D’Souza∗

NASA Ames Research Center, Moffett Field, CA, 94035

In the context of UAS Traffic Management operations, where varied vehicle types and

uncertain winds are expected, there is a need for a trajectory simulation that trades a

vehicle-centric approach with a systems level approach to model expected vehicle performance with an uncertain operational environment. This paper focuses on the first phase

of this trajectory prediction development. The goal is to understand the behavior of this

new framework as applied to a vehicle targeting a specific terminal altitude. Specifically,

a generalized six degree of freedom trajectory model was built to identify vehicle performance in the presence of wind. Generalization of this model was achieved by reducing the

number of simplifying assumptions and using vehicle performance parameters, that were

non-specific to the control system, to drive the control solution. The inner loop control

strategy behind this algorithm is to find the forces and torques required to guide the vehicle to an intended altitude within vehicle-specific force and attitude constraints in the

presence of winds. The control solution was optimized via the Artificial Bee Colony genetic

optimization method. It was successfully demonstrated that this framework can utilize vehicle performance parameters, that were non-specific to the control system, to optimize

a control solution for targeting an operational altitude. Additionally, vertical and lateral

path deviations from the nominal case were observed in the wind case results. This forms

the basis for quantifying the spatial and temporal requirements a vehicle will need for a

given operation.

I. Introduction

Small unmanned aircraft systems1 (sUAS) are becoming increasingly popular for commercial and governmental applications, e.g. fighting wildfires, agricultural monitoring, and surveillance.2 In order to

enable safe and efficient operations of these sUAS, NASA is conducting research into the technologies and

systems that may be required for UAS Traffic Management (UTM).3 The overarching principles of the UTM

research prototype are, 1) to provide flexibility where possible and structure where necessary, and 2) to

apply a risk-based approach where geographical needs and use cases determine the airspace performance

requirements. Implementation of these principles will require a quantitative understanding of a vehicle’s

operational capability (sense and avoid systems, etc.) and airspace demand. Addressing these aspects will

require a sUAS trajectory prediction model that is generalized to accommodate multiple vehicle types and

airspace environments (wind, terrain, etc.). In the absence of generalization, trajectory modeling would require knowledge proprietary information such as control systems and methods of a particular vehicle. Given

the large number and various types of sUAS expected to be operating in the airspace, requesting and implementing this information would prohibitive due to the constant modifications that would be required to

accommodate the various control laws and the potential legal barriers to acquiring proprietary information

from all operators.

In this paper, a trajectory prediction framework is developed to meet two primary requirements: 1) the

simulation must generalize vehicle and control models and 2) it must incorporate wind disturbances. There

is a prevalence of traffic management and path planning literature that focuses on simplified equations and

kinematic modeling. These methods handle multiple different fixed wing vehicles, reliably provide a solution,

∗Aerospace Flight Systems Engineer, Systems Analysis Office, NASA Ames Research Center, and AIAA Young Professional.

1 of 16

American Institute of Aeronautics and Astronautics

Downloaded by NASA AMES RESEARCH CENTER on February 15, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2017-0447 

 AIAA Information Systems-AIAA Infotech @ Aerospace 

 9 - 13 January 2017, Grapevine, Texas 

 AIAA 2017-0447 

 This material is declared a work of the U.S. Government and is not subject to copyright protection in the United States. 

 AIAA SciTech Forum 
Downloaded by NASA AMES RESEARCH CENTER on February 15, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2017-0447 

and incorporate some aspect of wind effects. NASA’s Air Traffic Management trajectory prediction algorithms use these methods, with various fidelity,4–6 to enable trajectory based operations. Path planners for

autonomous UAS use a decomposition approach, where assumptions are made to derive simplified kinematic

and algebraic equations.7 This allows for numerous, fast-time solutions to be generated and evaluated for

a specific vehicle. McGee and Hedrick8, 9 demonstrated that a mapping of the kinematic equations can be

used on fixed wing vehicles to account for the effects of wind. Many of these methods make simplifications

that narrow modeling of the vehicle performance envelope, which preclude non-linearities caused by wind

disturbances as the winds become more pronounced. Additionally, these trajectory algorithms do not have

to contend with a non-homogenous fleet of aircraft.

There are also numerous control system models for multi-rotor UAS and, in some models, disturbances

are implemented to determine vehicle performance. In Pounds et al.10 the flight dynamics of a quadrotor are

modeled using point mass, rigid body dynamics with a Proportional, Integral, Derivative (PID) controller,

that accounts for the effects of flapping, roll, and pitch damping, but no winds. In Hoffman et al.11 a

point mass model is implemented with a PID controller for path tracking and attitude control. In Raza and

Etele12 and Waslander and Wang13 the flight dynamics of a quadrotor, in the presence of wind gusts, is

modeled using a similar methodology to Hoffman et al.11 Additionally, a comprehensive survey of rotorcraft

UAS control methods14 identified PID controllers as one of the most widely used controllers, where its

capability is typically expanded via gain scheduling. In D’Souza et al.15 a gain-scheduled, PID controller

was implemented in a generalized trajectory prediction model and it was demonstrated that using a PID

controller as means to model vehicle performance in the presence of a sustained wind proved difficult because

of the inherent non-linearities of that exist for wind compensated vehicle dynamics.

In the context of UTM operations, where varied vehicle types and uncertain winds are expected, there

is a need for a trajectory simulation that trades a vehicle-centric approach with a systems level approach to

model expected vehicle performance with uncertainty in the context of an operational environment. This

paper focuses on the first phase of this trajectory prediction development. The goal is to understand the

behavior of this new framework as applied to a vehicle targeting a specific terminal altitude. Specifically,

a generalized six degree of freedom (6DOF) trajectory model was built to identify vehicle performance in

the presence of wind. Generalization of this model was achieved by: 1) limiting the number of simplifying

assumptions applied to the formulation to maintain fidelity, and 2) using vehicle performance parameters,

that were non-specific to the control system, to drive the control solution. The control solution was optimized

via the Artificial Bee Colony (ABC) genetic optimization method.16, 17 This method was selected because

it is a gradient-free search, immune to problems related to local optima. This characteristic is important

because stable muti-rotor control solutions lie in a small solution space and the simulation should not be at

risk of not converging to a solution. In Section II the trajectory simulation and optimization framework are

detailed, in Section III a nominal trajectory case and wind disturbance case are detailed, and Section IV

discusses future work in the context of UTM research.

II. Methodology

A 6DOF simulation was developed based on the equations provided in Beard18 and Zipfel19 to enable

modeling of fixed wing and multirotor vehicle dynamics. The work in this paper, however, focuses on a

multirotor test case as a starting point. The outer loop trajectory solution is driven by an inner loop

controls solution that includes a multi-segment step/linear force function and a constant/step function for

torques. The optimal inner loop control profile was found by wrapping the trajectory simulation in a genetic

optimization code. The objective function and constraints were chosen to minimize altitude error and ensure

that position rate and vehicle airspeed did not exceed defined limits. The effect of wind was incorporated

into the simulation as aerodynamic drag. This simulation provides the framework to incorporate vehicle

performance parameters and wind in order to quantify deviations of the vehicle from its nominal path.

The inner loop control strategy behind this algorithm is to find the forces and torques required to guide the

vehicle to an intended altitude within vehicle-specific force and attitude constraints in the presence of winds.

This strategy allows for generalized parameters to be used in the trajectory prediction, without knowing

exact details of the control system. Instead of requiring gains or knowledge of the control methodology

from a specific manufacturer, the simulation input parameters are based on quantities that define the vehicle

model and performance limits. The following list of parameters is only limited by parameters available for

the selected test case, but flexibility exists in this framework to include other parameters:

2 of 16

American Institute of Aeronautics and Astronautics
Downloaded by NASA AMES RESEARCH CENTER on February 15, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2017-0447 

• vehicle mass

• vehicle dimensions

• pitch and roll angle

• pitch, roll, and yaw acceleration

• altitude rate

• vehicle airspeed

• commanded forces and torques

A. Test Case

While the intent of this simulation is to model multiple vehicle types, development of a new simulation

requires an initial starting point. The Delta H (Fig. 1)20 quadrotor was selected as an initial test case because vehicle parameters and flight test data were available from flight tests in NASA’s first UTM research

prototype flight demonstration. The data from this flight test included altitude (z), latitude, longitude,

(a) Top View (b) Front View

Figure 1. Delta H drone

inertial velocity components(x, ̇ y, ̇ z ̇), vehicle rotation rate components (φ, θ, ψ), vehicle acceleration components (u, ̇ v, ̇ w ̇), and the quaternion (q0, q1, q2, q3). The flight data was recorded every 0.1 seconds. This flight

data was used to provide the basis for setting constraints initial conditions and terminal conditions in this

simulation. The flight data also provides a basis for validation of the simulation. Note that it is not in the

scope of this paper to complete a full validation but to develop an approach that demonstrates the potential

to generalize trajectory prediction modeling with wind.

The vehicle parameters included in this simulation were total mass, component masses (frame and body),

component dimensions (length, width, and height), CD, Aref , maximum velocity, maximum attitude angles,

and maximum angular accelerations. The component masses and dimensions are used in the calculation of

the moment of inertia matrix. The variables CD and Aref are used to calculate the drag. The maximum

attitude angles and angular accelerations are used as constraints in the simulation.

Some of these parameters were readily available from specification sheets and other online sources, however, some parameters required best estimates. Estimates were made for the component masses, component

dimensions, maximum angular accelerations, CD, and Aref . The total mass, maximum pitch angle, and

maximum roll angle were reported to NASA by the operator of the vehicle. The vehicle parameters are

detailed in Table 1.

The trajectory simulation was set-up to fly a take-off scenario, where the vehicle targets an altitude of

20 meters within 7 seconds. The termination altitude and time were determined from the flight data plotted

in Fig. 2. The mission scenario and vehicle parameters detailed in this section are direct inputs to the

equations of motion (detailed in Sec. II.B) and the optimal control strategy (detailed in Sec. II.C).

3 of 16

American Institute of Aeronautics and Astronautics
Downloaded by NASA AMES RESEARCH CENTER on February 15, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2017-0447 

Table 1. Vehicle Parameters

Parameter Description Value Units

mtotal total vehicle mass 4 kg

mratio ratio of frame mass to total mass 0.077 ND

mf rame vehicle frame mass mtotal ∗ mratio kg

mbody vehicle body mass mtotal − mframe kg

rframe radius length of the frame arms from center 0.057 m

[lbody, wbody, hbody] body (fuselage) dimensions [0.61, 0.23, 0.19] m

CD drag coefficient 0.5 ND

Aref vehicle reference area 0.14 m2

Vmax maximum vehicle airspeed 10 m/s

θmax, φmax maximum pitch and roll 20 deg

Figure 2. Altitude profile from flight data

B. Vehicle Dynamics

In order to capture the non-linearities of vehicle flight dynamics in the presence of wind, the vehicle dynamics

were modeled based on the 6DOF equations of motion derived in Beard18 and Zipfel.19 These equations

are generalized such that a fixed wing or multirotor sUAS could be modeled. The equations of motion are

derived using rigid body dynamics and kinematics. The derivation begins with the definition of the coordinate

systems for each rigid body, the Earth and the vehicle. Given the distance and operating altitudes of small

UAS, it is assumed that the Earth is flat and non-rotating. The Earth coordinate system, whose origin is

at the starting ground location of the vehicle, is an inertial frame where the x-axis aligns with North, the

y-axis aligns with East, and the z-axis points down toward the Earth. The origin of the body frame is set

at the vehicle’s center of mass, with the x-axis pointing out of the nose of the airframe, the y-axis positive

in the direction of the right of the vehicle and the z-axis positive pointing down and out of the bottom of

the aircraft.

The kinematics of the vehicle can be derived by applying a transformation matrix (1) which goes from

the inertial frame to the body frame.

⎡ ⎤

⎢

⎣sinφsinθcosψ − cosφsinψ

cosφcosψ cosθsinψ −sinθ

⎥ T

⎦

B

I = sinφsinθsinψ + cosφcosψ sinφcosθ (1)

cosφsinθcosψ + sinφsinψ cosφsinθsinψ − sinφcosψ cosφcosθ

The translational kinematics are found by multiplying the transpose of (1) by the body frame velocity

vector to find the velocity vector in the inertial frame (2),

⎤⎡ ⎤⎡

x ̇ u

⎢

⎣

⎥

⎦ = (TB

I )

T ⎢

⎣ ⎥ y ̇ v⎦ (2)

z ̇ w

4 of 16

American Institute of Aeronautics and Astronautics
Downloaded by NASA AMES RESEARCH CENTER on February 15, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2017-0447 

where [x, y, z] is the position vector in the inertial frame, [u, v, w] is the vehicle velocity vector in the body

frame, φ is the roll angle, θ is the pitch angle, and ψ is the yaw angle. The angular kinematics of the vehicle

are found using an intermediate transformation matrix to find the angular rates (3) in the intermediate

vehicle frames (see Beard18 for description),

⎤⎡ ⎡ ⎤ ⎤⎡

φ ̇ 1 sinφtanθ cosφtanθ

⎢

⎣

⎥

⎦

p

q

⎥

⎦ (3) ⎢⎣

⎥

⎦ = ⎢ θ

 ̇ ⎣0 cosφ −sinφ

 ̇ sinφ cosφ ψ 0 r cosθ cosθ

where p is the roll rate, q is the pitch rate, and r is the yaw rate in the inertial frame. The next step is to

apply Newton’s first and second law to find the translational and angular equations of motion as a function

of the forces and moments acting on the vehicle. The translational equation is found using Newton’s 1st

Law:

dv m = f (4) dtI

where v is the velocity vector of the vehicle in the body frame, m is the mass of the vehicle, f is the force

vector acting on the vehicle, and the subscript I indicates that the time derivative is in the inertial frame.

The equation of Coriolis is applied to find the inertial time derivative of the velocity vector in terms of the

body frame derivatives and rotations. This results in the following acceleration equations of motion:

⎤⎡ ⎡ ⎤ ⎤⎡

u ̇

⎥

⎦ = ⎢⎣

rv − qw

pw − ru⎥⎦ + ⎢⎣

fX

fY

⎥

⎦ (5) ⎢⎣

1

v ̇ m

w ̇ qu − pv fZ

where [u, ̇ v, ̇ w ̇ ] are the components of the vehicle acceleration in the body frame, fX is the force acting in

the X-axis of the body frame, fY is the force acting in the Y-axis of the body frame, and fZ is the force

acting in the Z-axis of the body frame. Finally, the angular accelerations are found by applying Newton’s

second law:

dhB

dtI

= τ (6)

where h is the angular momentum and τ is the applied torque vector. The equation of Coriolis is again

applied and results in the following angular acceleration equations of motion:

⎡ ⎤

J = ⎢

⎣

JX −JXY −JXZ

−JXY JY −JY Z

−JXZ −JY Z JZ

⎥

⎦ (7)

⎤⎡ ⎡⎛ ⎤ ⎤⎡ ⎤⎡ ⎞

p ̇ 0 −q

⎥

⎦

J ⎢

⎣

p

q

⎥

⎦ + ⎢⎣

τφ

τθ

⎟⎠⎥

⎦

r

⎢

⎣

⎥

⎦ = J−1 ⎢⎣ ⎜ q ̇ ⎝ −r 0 p (8)

r ̇ q −p 0 r τψ

where J is the moments of inertia matrix, τφ is the roll torque, τθ is the pitch torque, and τψ is the yaw

torque. Implementation of equations (2), (3), (5), and (8) with genetic optimization, for a multirotor, requires

additional assumptions and the use of the quaternion. Note that these assumptions and simplifications do

not preclude the modeling of a fixed wing vehicle because the trajectory simulation code was written in

the most generalized form ((2), (3), (5), and (8)). First, it is assumed that the vehicle is symmetric, which

simplifies (7) to a diagonal matrix where JXY = JXZ = JY Z = 0. The sum of the forces on the vehicle, (9),

~ ~ for a multirotor include the force from the motors (Fmotors), the gravitational force (Fg), and aerodynamic

drag ~(D) as a function of relative velocity . Note that all of these forces are in the body frame.

~ ~ ~ ~ f = Fmotors + Fg + D(Vwind) (9)

The quadrotor is treated as a point mass where the force from the motors acts at the center of mass, opposite

the body Z-axis. The drag force (13) is a quadratic with relative velocity and a constant drag coefficient

5 of 16

American Institute of Aeronautics and Astronautics
Downloaded by NASA AMES RESEARCH CENTER on February 15, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2017-0447 

applied for simplicity. A drag coefficient of 0.5 was selected based on an ellipsoidal body21 with a length-towidth ratio of 0.6. Since the wind velocity is in the inertial coordinate system, the transformation matrix

was applied to get wind velocity in the body frame (10),

⎤⎡ ⎤⎡

VX Vx

V~wind,B = ⎢

⎣VY

⎥

⎦ = TB

I

⎢

⎣Vy

⎥

⎦ (10)

VZ wind Vz wind

V~ ~ rel = Vvehicle,B − V~wind,B (11)

V~rel,sqr = V 2 ˆiB + V 2 ˆjB + V 2 k

ˆB (12) rel,X rel,Y rel,Z

⎤⎡

⎢

⎣

V 2

rel,X

V 2

rel,Y

⎥

⎦ (13) 1 D~ = − CDAref ρ 2

V 2

rel,Z

where C ~ D is the drag coefficient, Aref is the reference area of the vehicle, ρ is the air density, Vrel is the

~ relative wind velocity, and Vwind is the three dimensional wind vector. The gravity force is in the inertial

frame, so the transformation matrix was applied to get gravity force in the body frame:

⎤⎡

0

~ = TB Fg I

⎢

⎣

⎥ 0 ⎦ (14)

mg

At this point, the quaternion is introduced because it is well suited for near-Earth simulations and, due

to its relationship to the Euler angles (φ, θ, ψ), avoids singularities.19 For brevity the quaternion differential

equation, transformation matrix and Euler angle relations are included without derivation (see Chapter 4 of

Zipfel19 for an extended description of the quaternion). The four element quaternion differential equations

are a function of the Euler rates of the vehicle:

⎤⎡ ⎡ ⎤ ⎤⎡ −p −q −r

⎢

⎢

⎢

⎣

q ̇

0

q ̇

1

q ̇

2

⎥

⎥

⎥

⎦

1 = 2

⎢

⎢

⎢

⎣

0

⎢

⎢

⎢

⎣

⎥

⎥

⎥

⎦

q0

q1

q2

⎥

⎥

⎥

⎦

(15) p 0 r −q

q −r 0 p

q ̇

3 r q −p 0 q3

The initial quaternion state was calculated using the relationship of the Euler angles to the four element

quaternion:

⎡⎤⎡ ⎤ ψ θ φ cos cos + sin ψ sinθ sin φ

2 2 2

⎢

⎢

⎢

⎣

q0

q1

q2

⎥

⎥

⎥

⎦ =

⎢

⎢

⎢

⎣

cos 2 2 2

⎥

⎥

⎥

⎦

ψ θ φ cos sin φ − sin ψ sinθ cos 2 2 2 cos 2 2 2 (16) ψ φ θ sinθ cos + sin ψ cos sin φ

2 2 2 cos 2 2 2

θ φ ψ q3 sin ψ cos cos − cos sinθ sin φ

2 2 2 2 2 2

Next, using (1) and (16), the transformation matrix from the inertial frame to the body frame can be defined

in terms of the quaternion:

⎡ ⎤ 2 2 2 2 q0 + q1 − q2 − q 2(q1q2 + q0q3) 2(q1q3 − q0q2) 3

TB

I = ⎢

⎣

⎥

⎦ (17) 2 2 2 2

0 − q1 + q 2(q1q2 − q0q3) q 2 − q 2(q2q3 + q0q1) 3

2 2 2 2 2(q1q3 + q0q2) 2(q2q3 − q0q1) q0 − q1 − q2 + q3

Using the quaternion relationships and the sum of the forces acting on the vehicle, the complete set of

6DOF equations of motion are as follows:

⎤⎡ ⎡ ⎤ ⎤⎡ 2 2 2 2

0 + q1 − q2 − q 2(q1q2 + q0q3) 2(q1q3 − q0q2) x 3  ̇ q u

⎢

⎣ ⎥⎦ = ⎢⎣

⎢

⎣

⎥

⎦

⎥

⎦

2 2 2 2

0 − q1 + q y ̇ 2(q 2 − q (18) 1q2 − q0q3) q 2(q2q3 + q0q1) v 3

z ̇ 2(q1q3 + q0q2) 2(q2q3 − q0q1) 2 2 2 2 q0 − q1 − q2 + q3 w

6 of 16

American Institute of Aeronautics and Astronautics
Downloaded by NASA AMES RESEARCH CENTER on February 15, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2017-0447 

⎤⎡ ⎡ ⎡⎤ ⎤ ⎡⎛ ⎤ ⎡ ⎤⎞ V 2

rel,x

⎥

⎦ = ⎢⎣

rv − qw

pw − ru⎥⎦ + ⎢⎣

2g(q1q3 − q0q2)

2g(q2q3 + q0q1) ⎥

⎦ +

u ̇ 0

⎢

⎣

⎢

⎣ ⎜⎝ ⎥⎦ + ⎢⎣ ⎟⎠⎥⎦

1 1

V 2

rel,y F~ v ̇ 0 g − CDAref ρ (19) m 2 2 2 2 2 w ̇ qu − pv g(q0 − q1 − q2 + q V 2

3) −FZ,motor rel,z

⎤⎡ ⎡ ⎡⎤ ⎤ JY −JZ τφ

p ̇ qr JX JX

τθ

JY

τψ

⎢

⎣

⎥

⎦

⎢

⎣

⎥

⎦ + ⎢⎣

⎥

⎦

JZ −JX

JY q ̇ = pr (20)

JX−JY r ̇ JZ pq JZ

⎤⎡ ⎡ ⎤ ⎤⎡ −p −q −r

⎢

⎢

⎢

⎣

q ̇

0

q ̇

1

q ̇

2

⎥

⎥

⎥

⎦

1 = 2

⎢

⎢

⎢

⎣

0

⎢

⎢

⎢

⎣

⎥

⎥

⎥

⎦

q0

q1

q2

⎥

⎥

⎥

⎦

p 0 r −q (21)

q −r 0 p

q ̇

3 r q −p 0 q3

    2(q1q2+q0q3) 2(q2q3+q0q1) ψ = atan θ = asin(−2(q1q3 + q0q2)) φ = atan (22) q2

0+q21−q22−q23 q20−q21−q22+q23

The next step is to determine an optimal control strategy that will drive the trajectory solution such

that the vehicle meets the target objectives.

C. Optimal Control Framework: Design Variables and Objective Function

This section details the strategy for finding an optimal control solution for this problem. A gradient-free

genetic optimizer was utilized to avoid problems related to local optima when searching for an optimal set of

controls. Specifically, an ABC genetic optimization method was implemented. This method was was inspired

by the way bees search for food. In bee colonies, there are three roles: employed, on-looker, and scout. The

employed bees identify food sources and associated quality. This information is shared with the on-looker

bees who search for similar food sources with better quality. Meanwhile, the scout bees search, at random,

for new food sources and identify the quality. As applied to optimization, the solution (or set of design

variables) is the food source and the number of food sources considered are associated with the number of

employed bees. The associated value of the objective function identifies the quality of the solution. The

value of the objective function and the solution is passed to on-looker bees to iterate on the solution in search

of a better objective function value. At the same time, scouts are employed to look at an entirely different

set of solutions and capture solutions where the objective function is improved from the previous solutions.

For additional details regarding this method, reference papers by Abachizadeh et al.16 and Ghiglino et al.17

The optimization process is started by randomly selecting the design variables contained in a solution for

each solution set (or employed bee). Based on the resulting value of the objective function, the probability

of each solution to minimize the objective function is calculated. If the probability is high, the solution is

iterated on during the on-looker phase to minimize the objective function. If the probability is low, then

the scout phase is initiated to generate new solutions. For the full set of solutions, many cycles are run to

converge to a solution that minimizes the objective function. The following section discusses the multirotor

control strategy in the context of the design variables and objective function that drive the optimization.

1. Control Strategy Development

The inner loop control strategy was developed so that details of the vehicle control system are not needed.

In this simulation the control vector is defined as: [FZ,motor, τφ, τθ, τψ]. Identifying the forces and torques is

not a function of the number motors or the motor performance specifications, but is determined based on

the physics of multirotor flight and vehicle performance parameters. The forces and torques that result from

the optimization will be compared to the forces and torques derived from the flight data for verification.

In order to determine the force command, the force acceleration profile from the flight data was plotted

in Fig.3a and examined. This profile shows that the force acceleration initially rises to a value that is

approximately two times gravity. After approximately 1 second the force acceleration has ramped down and

maintains values ranging from 1.5 to 0.5 times gravity. With this understanding and through some trial

and error, a force acceleration function was developed that consisted of 4 segments and 3 switch points in

7 of 16

American Institute of Aeronautics and Astronautics
Downloaded by NASA AMES RESEARCH CENTER on February 15, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2017-0447 

Figure 3. Force acceleration flight data and analytical force profile (notional).

time (Fig. 3b). The force acceleration changes linearly in segments 1 and 3, while force acceleration is held

constant during segments 2 and 4.

Since the vehicle must overcome gravity and maintain altitude, the minimum value of force acceleration

(force) cannot be zero, but must be a non-zero factor of gravity. Therefore, each segment in the force

acceleration function is defined by a varied set of factors. The first factor is applied at the start of segment 1

and changes linearly up to segment 2, where the second factor is applied as a constant. The constant factor

applied over segment 4 provides the values by which the force acceleration in segment 3 changes linearly.

The time switches t1 and t2 are hard coded at 1 second and 2 seconds, respectively, based on the behavior in

the flight data. The time switch t3 becomes a design variable in the optimization. The force design variables

and associated equations are as follows,

h i



 

DVF C4

C2−C1

t (t − t0) + C1 ∗ g, 0.0 ≤ t ≤ t1 1−t0

C2 ∗ g, t1 < t ≤ t2

= C1 C2 t3

⎡

⎢

⎢

⎢

⎢

⎣

⎤

⎥

⎥

⎥

⎥

⎦

(23) FZ,motors

m =  C4−C2

t (t − t2) + C2 ∗ g, t2 < t ≤ t3 3−t2

C4 ∗ g, t3 < t ≤ t4

where Ci is the coefficient assigned to a particular segment, i indicates the segment number, ti is the time

switch that ends a particular segment. The maximum and minimum values (24) for the force acceleration

coefficients were selected such that the force value in segment 1 is always larger than g and the subsequent

force accelerations are between 1.5g and 0.1g to meet the altitude objective. The maximum and minimum

bounds on the time switch will lie within the flight time of the vehicle for the take-off segment.

h i

h

DVF,lower bound = 1.5 0.5 0.1 3[sec]

i (24) DVF,upper bound = 2.5 1.5 1.0 7[sec]

The next set of design variables are the torques due to the commanded attitude angles,

h i

DVτ = τφ τθ τψ (25)

After some test runs of the simulation, it was determined that optimizing multiple values of torque every

second resulted in poor ABC performance. Thus, torque profiles were set as constant values across the

trajectory since the objective is to target an altitude and not a specific x-y position.

The maximum and minimum bounds for each torque were calculated based on simplifying the attitude

acceleration equations (20). Using the assumption that attitude angles are small, provides a set of analytical

8 of 16

American Institute of Aeronautics and Astronautics
Downloaded by NASA AMES RESEARCH CENTER on February 15, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2017-0447 

(a) Roll acceleration (φ ̈) (b) Pitch acceleration (θ

 ̈) (c) Yaw acceleration (ψ ̈)

Figure 4. Euler angle acceleration profiles from flight data.

equations for calculating the upper bound on DVτ , (26).

 ̈ τφ,max = JX ∗ φmax

τθ,max = JY ∗ θ

 ̈

max (26)

 ̈ τψ,max = JZ ∗ ψmax

The maximum attitude accelerations were taken from the flight data, plotted in Fig. 4, and used to find

maximum torques in (26). The minimum values were determined to be the negative of the maximum torque

values found. The final bounds on the torques are as follows:

h i h i

DVτ,lower bound = −τφ,max −τθ,max −τψ,max = −350 −250 −125 N · m

h i h i (27) DVτ,upper bound = τφ,max τθ,max τψ,max = 350 250 125 N · m

To summarize, the design variables that fully define a control solution for this trajectory optimization include

force coefficients, a time switch, and torques (28)

h i

DV = C1 C2 C4 t3 τφ τθ τψ (28)

Since the focus of this test case is to define take-off performance of the vehicle and ensure that the attitude

and vehicle velocity limits are not exceeded, the objective function was defined based on final altitude error

and position rate violations,

2 OBJ = |zdesired − zactual| + jx

2

 ̇max,violations + jy

2

 ̇max,violations + jz

2

 ̇max,violations (29)

where zdesired is the desired altitude, and j indicates the number of violations for the subscript variable

identified. The altitude rate threshold, z ̇max, was set to 10 m/s to ensure that the altitude profile does not

oscillate largely in time, e.g. a parabolic rise to 40m and back down to 20m at the end of the trajectory. The

[  ̇xmax, y ̇max] threshold was set to a value of 0.01 m/s. This was selected to ensure that the inertial velocities

were small and as a result kept lateral path changes in x and y relatively small. While this appears to be a

counterintuitive number, this value was selected because test runs with this value resulted in better genetic

optimization performance. Finally, any solution that resulted in a vehicle airspeed that exceeded Vmax was

assigned a very large objective function value to ensure that the genetic optimization did not iterate on

solutions with infeasible airspeeds.

This control strategy does not rely on knowledge of the motor specifications or vehicle-specific control laws.

But is a function of vehicle performance parameters that identify the limits and constraints of the vehicle. The

design variables and the objective function are now implemented in the simulation to determine a nominal

and wind disturbed trajectory that achieves the altitude objective within these limits and constraints.

III. Results

The results that follow detail the altitude, velocities, attitudes, and control profiles for: 1) a nominal case

with no winds and 2) a wind disturbance case. The wind case exposes the vehicle to a three dimensional

9 of 16

American Institute of Aeronautics and Astronautics
Downloaded by NASA AMES RESEARCH CENTER on February 15, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2017-0447 

wind velocity vector, [Vx, Vy, Vz], where each component is equal to 5.0 m/s. The vehicle is exposed to this

wind for 1 second and the compensating control is initiated with a 0.5 second latency to model response

times for a control system.

A. Nominal Case

For the nominal case, the genetic optimizer found the following optimal control strategy:

h i

DVnom = 1.8295 0.9325 1.015 3.0[sec] 0.004016[N · m] 0.006351[N · m] 0.1490[N · m] (30)

In Fig. 5 the resulting force acceleration profile is plotted against the flight data and shows values that

are within the capability of the vehicle. In Fig. 6 the resulting constant values for each torque are plotted

Figure 5. Nominal force control profile optimized in the simulation.

against the flight data and are within the capability of the vehicle. The maximum force acceleration from

(a) Torque due to roll (b) Torque due to pitch (c) Torque due to yaw

Figure 6. Torques commanded from simulation.

the simulation is not as high as that from the flight data and the torques are close to zero, but not negligible.

Since the vehicle is targeting an altitude and is not attempting to maintain hard lateral path constraints,

torque modulation is not adversely affecting the altitude. Therefore, the maximum force in the first segment

does not require a full 2g0

s to compensate.

In Fig. 7a the simulated altitude profile shows that the vehicle meets the final desired altitude of 20 meters

within 7 seconds and has negligible oscillations. The lack of oscillations in altitude is further supported by

the altitude rate in Fig. 7b, which shows that the altitude rate peaks within the first second of the trajectory

and then maintains an almost constant value. The simulated vehicle velocity components, plotted in Fig. 8,

are below the Vmax limit of 10 m/s. Confirming that the constraints and limits implemented in the objective

function successfully ensured desired vehicle behavior. The simulated roll and pitch angles shown in Fig. 9

10 of 16

American Institute of Aeronautics and Astronautics
Downloaded by NASA AMES RESEARCH CENTER on February 15, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2017-0447 

(a) Altitude (b) Altitude Rate

Figure 7. Altitude and altitude rate from simulation.

Figure 8. Vehicle velocities from simulation.

are below the limit of 20 degrees due in part to the combination of the force acceleration profile with constant

torques and the objective function heuristic.

(a) Roll (b) Pitch

Figure 9. Roll and pitch attitude from simulation.

Finally, the simulated track for this control profile is plotted in Fig. 10. It should be noted that this

nominal case did not replicate the flight data ground track because the goal was to understand the behavior

of this framework for a simplified case, where the vehicle must meet a desired altitude. Having proved that

this simulation methodology works, modifications and improvements can be implemented to model path

tracking performance.

Although this is a simplified case, the resulting nominal case provides both spatial and temporal information for the intended operation (take-off to an altitude of 20 meters). This information becomes important

because it quantifies system level airspace demand and trajectory separation in the context of expected UTM

operations. The next step is to use this nominal case as a point of departure for the wind case, which is

detailed in the following section.

11 of 16

American Institute of Aeronautics and Astronautics
Downloaded by NASA AMES RESEARCH CENTER on February 15, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2017-0447 

Figure 10. Vehicle track from simulation.

B. Wind Case

Vehicle performance in the presence of wind must be identified because it has a direct impact on operational

safety and efficiency within UTM.

A three dimensional wind vector was simulated, where the wind velocity components, [Vx, Vy, Vz], were

set equal to 5.0 m/s, corresponding to a magnitude of 8.7 m/s. The wind disturbance was initiated two

seconds into flight and was terminated at three seconds. Additionally, a 0.5 second latency was implemented

to model the delay between detection of the disturbance and initiation of a control response to the wind.

Thus, the nominal control profile (31) is used for the first 2.5 seconds and the optimizer finds an optimal

control solution for [C4, t3, τφ, τθ, τψ] at time greater than 2.5 seconds, (32).

h i

DVwind,1 = 1.8295 0.9325 1.015 3.0 0.004016 0.006351 0.1490 , 0.0 ≤ t ≤ 2.5[sec] (31)

h i h i

DVwind,2 = C4 t3 τφ τθ τψ = 0.9727 2.5 −0.1675 0.006435 0.9394 , 2.5 < t ≤ 7.0[sec] (32)

In Fig. 11 the optimal force acceleration profile is plotted against the nominal case and shows that a slightly

smaller value of force acceleration was needed to compensate for the wind at time greater than 2.5 seconds.

This implies that the wind aided in the vehicle ascent such that the commanded force could be reduced. In

Figure 11. Force profile optimized to compensate for wind in the simulation.

Fig. 12 the optimal torques are plotted against the flight data and the nominal case. The roll torque was

reduced from the nominal case, the pitch torque is practically unchanged and the yaw torque is significantly

12 of 16

American Institute of Aeronautics and Astronautics
Downloaded by NASA AMES RESEARCH CENTER on February 15, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2017-0447 

higher than the nominal case. In Fig. 14 note that the largest deviations are in the u and v component

vehicle velocities. Since the algorithm is directed to maintain low velocities, most of the lateral torque control

effort is achieved via the roll and yaw torque. The yaw torque became significant because both u ̇ and v ̇ are

a function of yaw rate.

(a) Torque due to roll (b) Torque due to pitch (c) Torque due to yaw

Figure 12. Torques commanded from simulation.

In Fig. 13 the simulated altitude profile meets the final altitude of 20 meters within 7 seconds. The

altitude deviates from the nominal case as the vehicle compensates for the wind disturbance. The simulated

vehicle velocities, plotted in Fig. 14, are below the Vmax limit of 10 m/s, however, there are significant

oscillations. It is known that multirotor vehicles are unstable and it is possible that this optimal control

profile has induced an instability. Understanding the source of this behavior and modifying the control

profile will require further investigation.

(a) Altitude (b) Altitude Rate

Figure 13. Altitude and altitude rate from simulation.

Figure 14. Vehicle velocities from simulation.

The simulated roll and pitch angles plotted in Fig. 15 are below the limit of 20 degrees, however, the

profile is oscillatory at time greater than 2.5 seconds. The roll angle oscillations appear to be damping

out, but the pitch angle oscillations are increasing in magnitude. Again, understanding the source of this

behavior and modifying the control profile will require further investigation.

Finally, the nominal track for this control profile is plotted in Fig. 16. It is expected that the track would

deviate significantly from the nominal track because the current implementation does not guarantee lateral

13 of 16

American Institute of Aeronautics and Astronautics
Downloaded by NASA AMES RESEARCH CENTER on February 15, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2017-0447 

(a) Roll (b) Pitch

Figure 15. Roll and pitch attitude from simulation.

path tracking. However, the limitation on inertial velocities did keep the path deviations from the nominal

to less than 10 meters.

Figure 16. Vehicle track from simulation.

The resulting control strategy successfully compensated for the wind disturbance to meet the target

altitude. However, the control strategy exhibited an instability in the vehicle dynamics. In order for this

simulation to be an effective tool for future UTM operations, an optimal control solution must provide a

valid trajectory prediction but also correctly indicate when the vehicle cannot successfully compensate for

the wind disturbance.

IV. Conclusions and Future Work

A framework for a generalized trajectory prediction model, based on vehicle performance in the presence

of wind, was successfully implemented for targeting an operational altitude. It was found that vehicle

performance parameters, that were non-specific to the control system, could be used as means to determine a

set of optimal controls that guide the vehicle to its target. Additionally, vertical and lateral path deviations

from the nominal case were observed in the wind case results. This forms the basis for quantifying the

spatial and temporal requirements a vehicle will need for a given operation. In order to fully quantify

these requirements, the next steps in the development of this framework include implementation of a hover

condition constraint, targeting a 3D waypoint in time, and including additional phases of flight (cruise,

descent, etc.) in the simulation.

The ultimate goal of this framework it to provide UTM with a validated trajectory prediction tool

that can be easily modified for any mission and vehicle type and incorporates the effect of wind and other

uncertainties on the operation. The final output of the model would be a geo-fence that is defined by the

spatial and temporal path deviations that result from winds and other uncertainties. Figure 17 illustrates

this trajectory based geo-fencing and the current UTM prototype implementation. Dynamic calculation of

geo-fences across vehicle operations will allow UTM services to understand the demand on airspace so that

sense and avoid systems are able to maintain trajectory separation within a specific technological capability.

14 of 16

American Institute of Aeronautics and Astronautics
Downloaded by NASA AMES RESEARCH CENTER on February 15, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2017-0447 

Figure 17. Notional graph of geo-fence around an operation.

Acknowledgements

The author would like to acknowledge colleagues in the Systems Analysis Office and the Modeling and

Simulations Branch at NASA Ames Research Center for their helpful contributions to this work: John

Melton, Jaewoo Jung, and Heather Arneson.

References

1“Small Unmanned Aircraft Systems.” 14 ”CFR” 107.3. 2016. 2Federal Aviation Administration, FAA Aerospace Forecast: Fiscal Years 2016-2036, 2016. 3Kopardekar, P., Rios, J., Prevot, T., Johnson, M., Jung, J., and Robinson, J., “Unmanned Aircraft System Traffic

Management (UTM) Concept of Operations,” in “16th AIAA Aviation Technology, Integration, and Operations Conference,”

Washington, DC, 2016. 4Bilimoria, K. and Sridhar, B., “FACET: Future ATM Concepts Evaluation Tool,” in “3rd USA/Europe Air Traffic

Management R&D Seminar,” Napoli, Italy, 2000. 5Slattery, R. and Zhao, Y., “Trajectory Synthesis for Air Traffic Automation,” Journal of Guidance, Control, and Dynamics, Vol. 20, No. 2, 1997, pp. 232–238. 6Poles, D., “Base of Aircraft Data (BADA) Aircraft Performance Modeling Report,” Tech. Rep. EEC Technical/Scientific

Report No. 2009-009, EUROCONTROL, France, 2009. 7Goerzen, C., Kong, Z., and Mettler, B., “A Survey of Motion Planning Algorithms from the Perspective of Autonomous

UAV Guidance,” Journal of Intelligent and Robotic Systems, Vol. 57, 2010, pp. 65–100. 8McGee, T. and Hedrick, J. K., “Path Planning and Control for Multiple Point Surveillance by an Unmanned Aircraft in

Wind,” in “Proceedings of the 2006 American Control Conference,” Minneapolis, MN, 2006. 9McGee, T. G. and Hedrick, J. K., “Optimal Path Planning with a Kinematic Airplane Model,” Journal of Guidance,

Control, and Dynamics, Vol. 30, No. 2, 2007, pp. 629–633. 10Pounds, P., Mahony, R., and Corke, P., “Modelling and control of a large quadrotor robot,” Control Engineering Practice,

Vol. 18, 2010, pp. 691–699. 11Hoffman, G. M., Waslander, S. L., and Tomlin, C., “Quadrotor Helicopter Trajectory Tracking Control,” in “AIAA

Guidance, Navigation, and Control Conference and Exhibit,” . 12Raza, S. A. and Etele, J., “Simulation Tool for Testing and Validating UAV Autopilots in Wind Gust Environments,” in

“AIAA Atmospheric Flight Mechanics Conference,” . 13Waslander, S. L. and Wang, C., “Wind Disturbance Estimation and Rejection for Quadrotor Position Control,” in “AIAA

Infotech@Aerospace Conference,” . 14Kendoul, F., “Survey of Advances in Guidance, Navigation, and Control of Unmanned Rotorcraft Systems,” Journal of

Field Robotics, Vol. 29, No. 2, 2012, pp. 315–378. 15D’Souza, S., Ishihara, A., Nikaido, B., and Hasseeb, H., “Feasibility of Varying Geo-Fence around an Unmanned Aircraft

Operation based on Vehicle Performance and Wind,” in “2016 IEEE/AIAA 35th Digital Avionics Systems Conference (DASC),”

Sacramento, CA, 2016,

doi:TBD. 16Abachizadeh, M., Yazdi, M., and Yousefi-Koma, A., “Optimal Tuning of PID Controllers Using Artificial Bee Colony

Algorithm,” in “IEEE/ASME International Conference on Advanced Intelligent Mechatronics,” IEEE, Montreal, Canada, 2010. 17Ghiglino, P., Forshaw, J., and Lappas, V., “Online Evolutionary Swarm Algorithm for Self-Tuning Unmanned Flight

Control Laws,” Journal of Guidance, Control, and Dynamics, Vol. 38, No. 4, 2015, pp. 772–782.

15 of 16

American Institute of Aeronautics and Astronautics
Downloaded by NASA AMES RESEARCH CENTER on February 15, 2017 | http://arc.aiaa.org | DOI: 10.2514/6.2017-0447 

18Beard, R., “Quadrotor Dynamics and Control Rev 0.1,” Tech. Rep. Paper 1325, Brigham Young University All Faculty

Publications, 2008. 19Zipfel, P. H., Modeling and Simulation of Aerospace Vehicle Dynamics, Second Edition, American Institute of Aeronautics

and Astronautics, Inc., Washington, DC, 2007. 20Delta H, Delta Drone, http://droneinfini.fr/site/?page id=611, 2015. 21White, F. M., Fluid Mechanics, Fifth Edition, McGraw-Hill, New York, NY, 2003.

16 of 16

American Institute of Aeronautics and Astronautics