This thesis presents a solution that enables aerial robots to reason about surrounding wind flow fields in real time using on board sensors and embedded flight hardware. The core novelty of this research is the fusion of range measurements with sparse in situ wind measurements to predict surrounding flow fields. We aim to address two fundamental questions: first, the sufficiency of topographical data for accurate wind prediction in dense urban environments; and second, the utility of learned wind models for motion planning with an emphasis on energy efficiency and obstacle avoidance. Drawing on tools from deep learning, fluid mechanics, and optimal control, we establish a framework for local wind prediction using navigational LiDAR, and then incorporate local wind model priors into a receding-horizon optimal controller to study how local wind knowledge affects energy use and robustness during autonomous navigation. Through simulated demonstrations in diverse urban wind scenarios we evaluate the predictive capabilities of the wind predictor, and quantify improvements to autonomous urban navigation in terms of crash rates and energy consumption when local wind information is integrated into the motion planning. Sub-scale free flight experiments in an open-air wind tunnel demonstrate that these algorithms can run in real time on an embedded flight computer with sufficient bandwidth for stable control of a small aerial robot. Philosophically, this thesis contributes a new paradigm for localized wind inference and motion planning in unknown windy environments. By enabling robots to rapidly assess local wind conditions without prior environmental knowledge, this research accelerates the introduction of aerial robots into increasingly challenging environments.
Figures & tables
Figure 2.1 : Free body diagram of a rotary-wing UAV subjected to forces and torques generated by gravity and aerodynamic interactions with the air. The blue, red, and black arrows indicate coordinate frames, wrenches, and positional vectors, respectively.
Figure 2.2 : Empirically-driven power curves for a rotary-wing UAV as a function of ground speed and wind speed. On the left, power curves are plotted for varying headwinds (-) and tailwinds (+) with the max range indicated for each. On the right, the optimal range airspeeds and ground speeds are plotted versus the wind speed.
Parameter
Unit
Range (min-max)
m
kg
0.375–0.9375
cDx
N-(m/s) -2
0–1( 10−3 )
cDy
N-(m/s) -2
0–1( 10−3 )
cDz
N-(m/s) -2
0–2( 10−2 )
kd
N-rad-m-s -2
0–1.19( 10−3 )
kz
N-rad-m-s -2
0–2.32( 10−3 )
Table 3.1 : Randomized UAV and wind parameters for the Monte Carlo evaluation of the wind filter. The symbols are consistent with those in Chapter 2 .
Figure 3.1 : Monte Carlo evaluation of the Wind UKF over 50 simulations; each instance has randomized mass, drag coefficients, and average wind magnitudes.
Figure 3.2 : A simulated instance of the unscented Kalman filter estimating the local wind velocity vector for a UAV subject to Dryden wind gusts.
Figure 4.1 : Method overview for predicting local wind vector fields in real time without a map of the environment. Buildings are sensed using range sensors and fed through a deep neural network alongside sparse wind measurements to predict the components of the time-averaged wind flow field.
Map
CPU Time Per Timestep (s)
Real Time Factor
Training Map
( 500×500 m)
0.378
0.66
Evaluation Map
( 1,200×1,200 m)
1.606
0.16
Prediction Domain
( 25×25 m)
0.019
12.78
Table 4.1 : Run times for our fluid solver on varying grid sizes compared to our network on an Intel Core i5 @ 2.60 GHz.
Figure 4.2 : (a) Strouhal numbers for isolated rectangles of various shapes and onset flow angles compared with experimental data [ 113 , 112 ] . (b) Strouhal numbers for various aspect ratios compared with experimental data [ 113 , 112 ] . (c) Qualitative comparisons of our fluid solver against commercial solutions for isolated (bottom) and grouped (top) rectangles.
Figure 4.3 : An overview of the data acquisition, training, and evaluation of a wind flow decoder model. Thousands of wind simulations on small maps are densely sampled (left), and the labeled pairs are used to train the model via supervised learning (middle). The model is evaluated on randomly sampled locations in much larger maps (right) to measure generalization.
Figure 4.4 : Because of the large spatial scale of urban flow features, local croppings of urban flow fields can actually be thought of as evolving on a much lower dimensional manifold composed of elementary flows. The 30×30 m local flow fields on the left can be reproduced by unique linear combinations of the principal components from the right, some of which mimic fundamental flows in fluid mechanics such as uniform flow, vortex pairs, and shear flow.
Figure 4.5 : Cumulative explained variance as a function of the number of components retained from principal component analysis. As the region size grows, the local flow fields become more complex and therefore more components are required to recover the training data statistics.
Figure 4.6 : A global reconstruction of an urban wind flow field time-averaged over 50 seconds, generated by translating the robot in and around an evaluation map and overlapping 25×25 m predictions: (a) ground truth velocity magnitude; (b) reconstructed velocity magnitude; (c) velocity magnitude error; (d) velocity direction error. Points of interest for comparison are marked in dashed boxes highlighting strong shears (black), city wakes (red), and intersecting (pink) channel flows.
Figure 4.7 : Density scatter plots comparing between the predicted and actual spatial gradients of the wind flow field at the robot’s location for a model with H=25 m given both local and upstream wind measurements.
Figure 4.8 : The pointwise absolute error distributions for velocity magnitude (left) and direction (right) for H=25 m grid size models with different feature inputs, compared to baselines that assume a uniform flow field. Dots indicate the mean.
Figure 4.9 : The average absolute error in velocity magnitude (left) and direction (right) for models with varying prediction grid sizes ( H ), reported as the ratio of the prediction area to the maximum LiDAR perception region.
Figure 5.1 : Energy-aware autonomous navigation through an urban environment with dynamic wind conditions. The controller adapts in real time using on board wind flow predictions, balancing energy consumption and obstacle avoidance to reach the goal destination more robustly and with improved energy efficiency.
Figure 5.2 : Wind-aware navigation through unknown windy environments. A LiDAR Scan measures distances to nearby obstacles. The wind prediction module synthesizes a measurement of the local wind with the LiDAR scan to predict the time-averaged wind flow field in a local region, Lt , around the UAV. The stochastic receding horizon controller (MPPI Controller) produces a low frequency acceleration reference, which is then tracked with a higher frequency lower level controller (Lower-Level Controller). The simulation environment includes full rigid body dynamics and unsteady wind flow fields.
Figure 5.3 : Time-averaged wind flow fields for 4 out of the 16 wind scenarios used for experiments on one of the 16 maps. Each scenario is labeled with one of the four possible start (green circle) and goal (red square) waypoints.
Parameter
Value
Unit
Nτ
100
-
Δt
0.25
s
Nk
187
-
fplan
10
Hz
λ
0.79
-
σu
5.15
m/s 2
Table 5.1 : Hyperparameters for the simulated nominal MPPI controller.
Figure 5.4 : Average percent reduction in energy consumption with one standard deviation regions for EA-MPPI+W relative to MPPI against different energy cost weights (left) and wind prediction region sizes (right). Trials where either EA-MPPI+W or MPPI failed are not included.
Ablation
Avg. Energy (kJ)
Cmd. Acc. (m/s 2 )
Time to Goal (s)
Distance Traveled (m)
Failure Rate ( collision , hover )* (%)
MPPI
13.61
20.95 ± 1.48
28.29 ± 4.51
329.54 ± 11.29
36.3 (22.4,14.0)
MPPI+W
14.16
21.35 ± 1.85
28.41 ± 4.52
329.94 ± 11.98
28.8 (18.9,9.9)
EA-MPPI
11.07
14.15 ± 1.04
38.31 ± 4.90
329.15 ± 9.07
18.7 (2.5,16.1)
EA-MPPI+W
10.81
14.05 ± 0.62
37.94 ± 5.63
334.18 ± 11.66
15.8 (2.6,13.2)
Table 5.2 : Performance metrics of different ablations of the MPPI controller across simulated scenarios.
Figure 5.5 : The best trial with respect to energy reduction over the baseline from the simulated wind scenarios. The blue lines correspond to MPPI whereas the green lines correspond to EA-MPPI+W . The wind flow field is seen averaged over time in the top right alongside the paths each ablation took. Airspeed and instantaneous power consumption are shown in the bottom two plots.
Figure 5.6 : The worst trial with respect to energy reduction over the baseline from the simulated wind scenarios. The blue lines correspond to MPPI whereas the green lines correspond to EA-MPPI+W . The wind flow field is seen averaged over time in the top right alongside the paths each ablation took. Airspeed and instantaneous power consumption are shown in the bottom two plots.
Figure 5.7 : Case studies demonstrating the qualitative differences between the local receding horizon (MPPI and EA-MPPI+W) and the global graph search (GS and GS+W) methods for navigation in windy urban environments. The EA-MPPI+W method has parameters QE=3.0 and H=90 . The background contours visualize the time-averaged wind magnitude, although in simulation the winds vary over time.
Figure 6.1 : An overview of the experimental design, in which a small UAV is tasked with tracking a dynamic goal position through a small-scale city-like environment subjected to wind generated by the Ames WindShaper fan array. The left and right sides of the diagram describe the hardware and software components, respectively. The software components are split between two off board computers: a Beelink Mini PC handles the LiDAR simulation, networking, and lower level control, while a flight-ready Nvidia Jetson AGX Xavier runs the core wind prediction and planning algorithms.
Figure 6.2 : The WindShaper facility at NASA Ames Research Center has a 9 ft by 6 ft open air wind tunnel with 1,134 individually programmable fans capable of producing winds up to 16 m/s. The wind tunnel is positioned inside a large motion capture space which is used to track the pose of objects with centimeter precision. On the right, the six experiment configurations are visualized with the sub-scale buildings highlighted in red. See Section B.1 for pictures of each obstacle with dimensions.
Figure 6.3 : A modified palm-sized Crazyflie 2.1-Brushless UAV designed and manufactured by Bitcraze. The UAV has a rotor diameter of 55 mm and is modified with IR reflective markers, a dedicated sensor for high frequency power measurements, and a custom whisker-based airflow sensor to measure wind.
Vehicle
Mass (kg)
Rotor Diameter (m)
Wingspan (m)
Max Speed (m/s)
Crazyflie 2.1-Brushless*
0.057
0.055
-
-
Joby S4 Prototype**
2000†
1.8
10.6
89.4
Table 6.1 : A rough comparison between the sub-scale flyer (Crazyflie 2.1-Brushless) and the Joby S4 prototype aircraft.
Figure 6.4 : The power spectral density (PSD) of the wake velocity behind the 16×16 inch obstacle, as measured by the WindProbe positioned roughly 0.4 m downstream from the obstacle’s leeward side while the WindShaper operated at 20% power. For additional context, Kolmogorov’s - 5/3 turbulence scaling law [ 130 ] is plotted for the higher frequencies.
Figure 6.5 : Identification of relationships between thrust, motor speed, input voltage, and commanded pulse width modulation (PWM) for a single rotor of the Crazyflie 2.1-Brushless. The bench top experimental setup is seen in the inset.
Figure 6.6 : Identification of the induced drag coefficient, kd , for the Crazyflie via linear regression, following the procedure outlined in [ 25 ] .
Figure 6.7 : Data used for the identification of the thrust kη , inflow kz , and translational lift kh coefficients for the Crazyflie, following the procedure outlined in [ 25 ] .
Parameter
Variable
Value
Unit
Thrust Coefficient
kη
4.26(10−8)
(rad/s)2N
Induced Drag Coefficient
kd
5.09(10−6)
m−radN−s2
Inflow Coefficient
kz
1.19(10−5)
m−radN−s2
Translational Lift Coefficient
kh
1.11(10−3)
m2Ns2
Table 6.2 : Lumped parameter aerodynamic coefficients identified using the procedure in [ 25 ] . The magnitudes of each coefficient are indicated in the parenthesis.
Figure 6.8 : Identification of a cubic power consumption model as a function of airspeed for the Crazyflie. The data is annotated with the hover power and max range and endurance speeds (assuming no wind). The max range speed was determined via the range tangent line, which extends down to the point (0, 0).
Figure 6.9 : The root mean square error (RMSE) between the CMU Whisker and Wind UKF across all the experiments. In the left plot, the numbers correspond to the experiment IDs in Table B.1 .
Figure 6.10 : Time-averaged wind flow field reconstructions at z=0.75 for experiment 14 of the Flow Traversal series. The CMU Whisker and Wind UKF measurements were collected from the UAV, while the WindProbe was collected during a separate sweep.
Figure 6.11 : Comparisons of the instantaneous wind estimates between the UKF and CMU Whisker for Flow Traversal experiment 14. The estimate provided by the WindProbe represents a sampling of the mean and standard deviation of the flow field from the bins corresponding to the UAV’s trajectory.
Figure 6.12 : Local time-averaged velocity magnitude comparisons between reconstructions using the WindProbe and the instantaneous predictions from the wind decoder network running on an Nvidia Jetson AGX Xavier. The green sphere represents the desired position, the triad is the pose of the UAV, and the green and magenta arrows extending off screen are wind estimates from CMU Whisker and Wind UKF , respectively. The wind decoder network takes as input the Wind UKF wind estimate as well as the simulated point cloud data (red dots) to sense topography. The top 25 forecasted trajectories from the wind-aware motion planner can also be seen originating the UAV’s position.
Parameter
Value
Unit
Nτ
1000
-
Δt
0.04
s
Nk
25
-
fplan
25 (target)
Hz
λ
5.0
-
σu
N/A
m/s 2
Table 6.3 : Hyperparameters for the MPPI controller implemented on an Nvidia Jetson AGX Xavier. Note that because the actuator sampling scheme changed with the hardware implementation, the associated parameters below were not used.
Category
Routine
Time (ms)
% Contribution
GPU Compute
update_mppi
25.5
36.4%
Cost Calculation
compute_cost
2.5
-
Rollouts
(remainder)
23.0
-
Memory Transfer
Tensor.cpu()
44.4
63.4%
Overhead
numpy/torch operations
0.1
0.2%
Total
70.0 ms
100%
Table 6.4 : Breakdown of computational latency for a single MPPI planning iteration on the Nvidia Jetson AGX Xavier. Data is averaged over 25 planning cycles, and the raw profiling log is available in Appendix B .
Figure 6.13 : Time series measurements of the distance to the nearest obstacle for the wand and UAV while running the MPPI motion planner. The UAV tracks the wand adequately when it is far enough away from obstacles. However, the obstacle avoidance term in the MPPI cost function prevents the UAV from tracking the wand when it gets too close to an obstacle.
Figure 6.14 : Heat maps of the standard deviation of the lateral force commands and the average body rotation rates of the UAV for experiment 11 , demonstrating regions of unsteadiness in urban-like winds.
Appendix figures & tables8 assets
Supplementary material from the paper’s appendix.
Appendix
Robot
Nominal Mass ( kg )
Rotor Diameter ( m )
Acc Uncertainty ( m/s2 )
Drag Coefficient ( Ns/m )
X/Y Wind Uncertainty ( m/s )
Falcon 4
4.2
0.381
0.0105 †
0.55*
0.82
Crazyflie 2.1
0.030
0.045
0.0105 †
0.01*
0.28
Ingenuity
1.61
1.61
0.0105 †
0.03**
6.03
Appendix
Table A.1 : Example wind uncertainties associated with three robots with varying size, mass, drag coefficients, and operating environments.
Figure B.1 : The four obstacles used during experiments in the WindShaper facility. The “slim” (top) and “wide” (bottom) obstacles were constructed from recycled 1 inch thick plywood and fitted with motion capture marker arrays for pose tracking and tuft grids for flow visualization.
Figure B.2 : Characterization of the average velocity of the bulk wind flow generated by the WindShaper with varying fan power levels and fan speeds. Velocity measurements were taken at three downstream locations to confirm that the wind maintains a constant speed throughout the test section.
Figure B.3 : Validation of the WindShaper characterization on two experiments. On the left, the WindShaper receives a sequence of discrete increasing power commands, whereas on the right the command is a continuous sinusoid. The predicted wind from the fan speeds (magenta lines) is in good agreement with the data collected from the WindProbe .
Figure B.4 : The whisker sensor calibration procedure: A) calibration curves for each of the sensor and wind axes; B) oblique view of the whisker sensor rigidly mounted on top of the anemometer; C) side view of the whisker experiencing maximum deflection due to the wind; D) side view of the whisker with a WindShaper fan array in the background.
Exp. ID
Experiment Series
# Trials (SE3/MPPI)
Building Config.
Fan Power (%)
A1
Aero. Characterization
1 (1/0)
B0
0-50
1-2
Flow Traversal
4 (2/2)
B0
10, 30
3-5
Flow Traversal
6 (3/3)
B1
10, 25, 35
6-8
Flow Traversal
7 (4/3)
B2
10, 25, 35
9-11
Flow Traversal
6 (3/3)
B3
10, 25, 35
12-14
Flow Traversal
6 (3/3)
B4
10, 25, 35
Appendix
Table B.1 : A complete list of the experiments conducted at the NASA Ames WindShaper facility. The building configurations are shown in Figure 6.2 . The number of trials per experiment are noted and broken down into which controller was running.
Figure B.5 : Part 1 of the flow field reconstructions and unsteadiness metrics for the Flow Traversal experimental series.
Figure B.6 : Part 2 of the flow field reconstructions and unsteadiness metrics for the Flow Traversal experimental series.
Aerial robotics studies the design, perception, control, and coordination of flying robots for navigation, inspection, transport, and physical interaction.