Position Constraints in Aircraft AeroSandbox Trajectory Optimization

[michgoer]

I am encountering inconsistent results when applying position constraints in my AeroSandbox model during an aircraft optimization task for the Design Build Fly (DBF) 2023 contest. Despite following what I believe to be the correct approach, the constraints do not always enforce the intended trajectory through the specified waypoints, leading to significant variations between optimization runs. Additionally, the results are highly sensitive to changes in time steps. Occasionally, at the last waypoint, the aircraft slows down and makes unexpected turns, which seems to be an attempt to prolong the flight duration.

My question is whether I am implementing the position constraints correctly, and if not, how a correct implementation should look. Below, I have included the optimization code and the result plotting code:

Optimization Code:

import aerosandbox as asb                                   # AeroSandbox for aircraft design and analysis
import aerosandbox.numpy as np                              # AeroSandbox's numpy variant
import aerosandbox.tools.pretty_plots as p                  # Pretty plots tools from AeroSandbox
import base64                                               # For encoding binary data into ASCII text
import cadquery as cq                                       # For 3D parametric modeling
import casadi as cas                                        # Symbolic framework for numeric optimization
import json                                                 # For handling JSON data
import matplotlib                                           # Base library for plotting
import matplotlib.pyplot as plt                             # For plotting
import neuralfoil as nf                                     # For rapid aerodynamic analysis of airfoils
import os                                                   # For interacting with the operating system
import pandas as pd                                         # For data manipulation and analysis
import re                                                   # For regular expression operations
import scipy                                                # For scientific and technical computing
import scipy.constants as sc                                # For physical and mathematical constants
from scipy.interpolate import interp1d                      # For 1D interpolation
from IPython.display import Image, display, HTML            # For displaying images and HTML in Jupyter notebooks
from matplotlib.colors import LinearSegmentedColormap       # For custom color maps in Matplotlib
from matplotlib.ticker import AutoMinorLocator              # For finer control over plot ticks
from scipy.interpolate import CubicSpline                   # For cubic spline interpolation
from aerosandbox.tools.string_formatting import eng_string  # For engineering number formatting
from aerosandbox.tools import units as u                    # For unit conversions

opti = asb.Opti()

# ==================
# Airfoil Optimization Setup
# ==================
# Initial airfoil guess based on the SD7062 profile
initial_guess_airfoil = asb.KulfanAirfoil("sd7062")
initial_guess_airfoil.name = "Dragonfly Airfoil (SD7062)"

# Optimized airfoil definition with variable parameters
optimized_airfoil = asb.KulfanAirfoil(
    name="Optimized",
    lower_weights=opti.variable(
        init_guess=initial_guess_airfoil.lower_weights,
        lower_bound=-0.5,
        upper_bound=0.25,
    ),
    upper_weights=opti.variable(
        init_guess=initial_guess_airfoil.upper_weights,
        lower_bound=-0.25,
        upper_bound=0.5,
    ),
    leading_edge_weight=opti.variable(
        init_guess=initial_guess_airfoil.leading_edge_weight,
        lower_bound=-1,
        upper_bound=1,
    ),
    TE_thickness=0,  # Trailing edge thickness set to 0
)

# Create an Airfoil object with these coordinates
wing_airfoil = optimized_airfoil

# Wing twist variable
wing_twist = opti.variable(init_guess=0, lower_bound=-10, upper_bound=10)

# Define airfoil for the tail
tail_airfoil = asb.Airfoil("naca0009")

# Antenna Dimensions
antenna_length = opti.variable(init_guess=0.8, lower_bound=0.05, upper_bound=0.95)  # Length in meters to fit into the prescribed box
antenna_diameter = 0.021336  # Outer diameter in meters (0.840 inches)
num_sections = 30  # Number of sections to model the antenna

# Define the airplane model
airplaneDFY = asb.Airplane(
    name="DBF23 Dragonfly",
    wings=[
        # Main Wing Definition
        asb.Wing(
            name="Main Wing",
            symmetric=True,
            xsecs=[
                # Wing Cross Sections: Root, Mid, and Tip
                asb.WingXSec(xyz_le=[-0.025, 0, 0], chord=0.40, twist=wing_twist, airfoil=wing_airfoil),
                asb.WingXSec(xyz_le=[-0.025, 0.5, 0], chord=0.40, twist=wing_twist, airfoil=wing_airfoil),
                asb.WingXSec(xyz_le=[-0.025, 1, 0], chord=0.40, twist=wing_twist, airfoil=wing_airfoil),
            ]        
        ),
        asb.Wing(
            name="Horizontal Stabilizer",
            symmetric=True,
            xsecs=[
                asb.WingXSec(  # root
                    xyz_le=[0, 0, 0],
                    chord=0.2,
                    twist=0,
                    airfoil=tail_airfoil,
                ),
                asb.WingXSec(  # tip
                    xyz_le=[0, 0.3, 0],
                    chord=0.2,
                    twist=0,
                    airfoil=tail_airfoil
                )
            ]
        ).translate([0.923, 0, 0]),
        asb.Wing(
            name="Vertical Stabilizer",
            symmetric=False,
            xsecs=[
                asb.WingXSec(
                    xyz_le=[0, 0, 0],
                    chord=0.304,
                    twist=0,
                    airfoil=tail_airfoil,
                ),
                asb.WingXSec(
                    xyz_le=[0.131, 0, 0.35],
                    chord=0.173,
                    twist=0,
                    airfoil=tail_airfoil
                )
            ]
        ).translate([0.82, 0, 0])
    ],
    fuselages=[
        # Fuselage Tail Boom Definition
        asb.Fuselage(
            name="FuselageTAILBOOM",
            xsecs=[
                asb.FuselageXSec(
                    xyz_c=[0.0602 * xi, 0, 0],
                    width=0.025,
                    height=0.025,
                    shape=7
                )
                for xi in range(10)
            ]
        ).translate([0.4, 0, -0.0125]),

        # Fuselage Body Definition
        asb.Fuselage(
            name="FuselageBODY",
            xsecs=[
                asb.FuselageXSec(
                    xyz_c=[0.042 * xi, 0, 0],
                    width=0.1,
                    height=0.1,
                    shape=7
                )
                for xi in range(10)
            ]
        ).translate([-0.25, 0, -0.05]),

        # Fuselage Nose 
        asb.Fuselage(
            name="FuselageNOSE",
            xsecs=[
                asb.FuselageXSec(
                    xyz_c=[0.22 * xi, 0, 0],
                    width=-0.1 * np.sqrt(xi),  # Adjust the width to be less pointy
                    height=-0.1 * np.sqrt(xi),  # Adjust the height to be less pointy
                    shape=7
                )
                for xi in np.linspace(0, 1, 40)
            ]
        ).translate([-0.47, 0, -0.05]),

        # Transition Section from Fuselage Body to Tail Boom
        asb.Fuselage(
            name="FuselageREAR",
            xsecs=[
                asb.FuselageXSec(
                    # Adjust the center point of each cross-section to slope the bottom surface
                    xyz_c=[0.27 * xi, 0, 0.036 * xi],  # Move the center point up for the slope
                    width=np.interp(xi, [0, 1], [0.1, 0.025]),
                    height=np.interp(xi, [0, 1], [0.1, 0.025]),
                    shape=7
                )
                for xi in np.linspace(0, 1, 40)
            ]
        ).translate([0.13, 0, -0.05]),

        # Antenna modeled with rectangular cross-sections in the x-direction
        asb.Fuselage(
            name="Antenna",
            xsecs=[
                asb.FuselageXSec(
                    xyz_c=[0.075 + xi * antenna_diameter, 1, antenna_length / 2],  # Position in x, y, z
                    width=antenna_diameter * np.sin(np.pi * xi),  # Width of the cross-section
                    height=antenna_length,  # Height of the cross-section
                    shape=1000  # Shape parameter for rectangular shape
                )
                for xi in np.linspace(0, 1, num_sections)
            ]
        ),
    ],
    s_ref= 1.85
)

# Add optimization constraints
opti.subject_to([
    # Local thickness constraints for structural integrity
    optimized_airfoil.local_thickness(x_over_c=0.33) >= 0.128,
    optimized_airfoil.local_thickness(x_over_c=0.90) >= 0.014,
    # Trailing edge angle for aerodynamic efficiency
    optimized_airfoil.TE_angle() >= 6.03,
    # Surface curvature constraints to ensure smooth airfoil shape
    optimized_airfoil.lower_weights[0] < -0.05,
    optimized_airfoil.upper_weights[0] > 0.05,
    optimized_airfoil.local_thickness() > 0
])

# Function to control the 'wiggliness' of the optimized airfoil
get_wiggliness = lambda af: sum([
    np.sum(np.diff(np.diff(array)) ** 2)
    for array in [af.lower_weights, af.upper_weights]
])

# Constraint to minimize 'wiggliness' for manufacturing and performance
opti.subject_to(
    get_wiggliness(optimized_airfoil) < 4 * get_wiggliness(initial_guess_airfoil)
)

# Define time, with horizon length unknown
time = np.linspace(
    0,
    opti.variable(
        init_guess=32, 
        lower_bound=2, 
        upper_bound=40),
    80
)

# Define the number of points in the trajectory
N = np.length(time) 

# Define the total time of the trajectory, to be later used as the objective function
total_time = time[-1]

# Constants specific to the DBF 2023 Mission 2
mission_time = 10 * 60  # s
payload_min_weight_fraction = 0.30  # Dimensionless

# Constants specific to DBF 2023 Mission 3
mission3_laps = 3

# Other aircraft variables
total_aircraft_weight = opti.variable(init_guess=30, log_transform=True, upper_bound=300)  # Aircraft weight in Newtons
total_aircraft_weight_kg = total_aircraft_weight / 9.81  # Aircraft weight in kg
payload_weight = total_aircraft_weight_kg - 6

# Define bounds for altitude and calculate atmospheric properties
altitude = opti.variable(init_guess=np.linspace(0, 40, N), lower_bound=30, upper_bound=40)

# Creation of the dynamics instance for 3D motion of the aircraft
dyn = asb.DynamicsPointMass3DSpeedGammaTrack(
    speed=opti.variable(init_guess=20, n_vars=N, lower_bound=2, upper_bound=50),
    track=opti.variable(init_guess=np.linspace(0, np.pi, N)),
    mass_props=asb.MassProperties(mass=total_aircraft_weight_kg),
    x_e=opti.variable(init_guess=np.linspace(-201, 201, N)),
    y_e=opti.variable(
        init_guess=np.linspace(-40, 40, N),
        lower_bound=-20,
        upper_bound=40,
    ), 
    z_e=-altitude,
    gamma=opti.variable(
        init_guess=np.zeros(N),
        lower_bound=-np.pi/2.1,
        upper_bound=np.pi/2.1,
    ),
    alpha=opti.variable(init_guess=1, n_vars=N, lower_bound=-15, upper_bound=15),
    beta=np.zeros(N),
    bank=opti.variable(
        init_guess=np.zeros(N),
        lower_bound=-np.pi/3,
        upper_bound=np.pi/3,
    ),
)

# Define the operating point for the dynamics
op_point = asb.OperatingPoint(
    atmosphere=asb.Atmosphere(altitude=altitude + 728), #Tucson, Arizona
)

# Define the start 360 waypoints to optimize the turn radius
turn_start_360_x = opti.variable(init_guess=-50, lower_bound=-80, upper_bound=-30)
turn_start_360_y = opti.variable(init_guess=20, lower_bound=17, upper_bound=23)

# Define the end 360 waypoints to optimize the turn radius
turn_end_360_x = opti.variable(init_guess=-10, lower_bound=-28, upper_bound=20)
turn_end_360_y = opti.variable(init_guess=20, lower_bound=17, upper_bound=23)

# Define straight 500 ft y-coordinate as an optimization variable
straight_500_ft_y = opti.variable(init_guess=5, lower_bound=0, upper_bound=10)

# Define return straight end y-coordinate as an optimization variable
return_straight_end_y = opti.variable(init_guess=5, lower_bound=0, upper_bound=10)

# Define waypoints to complete one lap
waypoints = {
    'start': (0, 0),
    'straight_500_ft': (152.4, straight_500_ft_y),  # Straight 500 ft waypoint
    '360_turn_start': (turn_start_360_x, turn_start_360_y),  # Start of 360 turn
    '360_turn_end': (turn_end_360_x, turn_end_360_y),  # End of 360 turn
    'return_straight_end': (-152.4, return_straight_end_y)  # End of return straight
}

# Apply constraints for each waypoint to be reached at the calculated index
for i, (x, y) in enumerate(waypoints.values()):  # Iterate through waypoints
    idx = min(i * (N // (len(waypoints))), N)
    opti.subject_to([
        dyn.x_e[idx] == x,  
        dyn.y_e[idx] == y
    ])

# Apply the gravitational force
dyn.add_gravity_force(g=9.81)

# Define the airplane model
airplane = airplaneDFY

# Define the aerodynamic model and compute aerodynamic forces
aero = asb.AeroBuildup(airplane=airplaneDFY, op_point=dyn.op_point).run()

# Apply forces
dyn.add_force(
    *aero["F_b"], axes="body"
)

# Define variable thrust
thrust = opti.variable(init_guess=30, n_vars=N, lower_bound=0, upper_bound=50)

# Apply variable thrust force in the x-direction
dyn.add_force(
    Fx=thrust, axes="body"
)

# Initial and final state constraints
opti.subject_to([
    dyn.track[0] == 0,
    dyn.gamma[0] == 0,
    dyn.bank[0] == 0,
    dyn.x_e[-1] == 0,
    dyn.y_e[-1] == 0,
])

# Add constraints to the altitude
opti.subject_to(
    dyn.altitude / 1000 > 0,
)

number_of_laps = mission_time / total_time

# ==================
# Objective Function: Maximize Mission Score
# ==================
# Maximum scores based on the winning team's statistics
max_payload_weight_m2 = 7.5  # kg
max_laps_flown_m2 = 12
max_score_m2 = max_payload_weight_m2 * max_laps_flown_m2  # Maximum payload weight * laps flown

max_antenna_length_m3 = 0.97  # meters
max_mission_time_m3 = 3 * 79 # seconds
max_score_m3 = max_antenna_length / max_mission_time_m3  # Maximum antenna length / mission time

mission3_time = total_time * mission3_laps

# Mission score calculation
normalized_score_m2 = ((payload_weight * number_of_laps) / max_score_m2)

# Normalize the score for mission 3
normalized_score_m3 = (antenna_length / mission3_time) / max_score_m3

# Finalize the problem by applying the dynamics constraints to the optimization problem
dyn.constrain_derivatives(opti, time)

# Minimize total distance/average speed (which is equivalent to total time)
opti.minimize(-(normalized_score_m3 + normalized_score_m2))

# Defining callback function to print the values of the variables at each iteration for debugging
def callback(i):
    time_history = opti.debug.value(total_time)
    x_history = opti.debug.value(dyn.x_e)
    y_history = opti.debug.value(dyn.y_e)
    z_history = opti.debug.value(dyn.z_e)
    thrust_history = opti.debug.value(thrust)
    print(f"Thrust: {thrust_history}")
    print(f"X: {x_history}")
    print(f"Y: {y_history}")
    print(f"Z: {z_history}")
    print(f"Time: {time_history}")
    print(f"    ")

# Solve the optimization problem
sol = opti.solve(
    callback=callback,
    max_iter=1000,
)

# Extract the optimized values of the variables
dyn = sol(dyn)
optimized_airfoil = sol(optimized_airfoil)
optimized_twist = sol(wing_twist)
antenna_length = sol(antenna_length)

# Output results
print(f"Optimized total time: {sol.value(time)} s")
print(f"Flight path x-coordinates: {sol.value(dyn.x_e)}")
print(f"Flight path y-coordinates: {sol.value(dyn.y_e)}")
print(f"Optimized altitude: {sol.value(dyn.altitude)} m")
print(f"Optimized track angles: {sol.value(dyn.track)}")
print(f"Optimized flying speed: {sol.value(dyn.speed)} m/s")
print(f"Normalized mission score: {sol.value(normalized_score_m2)}")
print(f"Total aircraft weight: {sol.value(total_aircraft_weight)} N")
print(f"Total aircraft weight: {sol.value(total_aircraft_weight_kg)} kg")
print(f"Payload weight: {sol.value(payload_weight)} kg")
print(f"Number of laps: {sol.value(number_of_laps)}")
print(f"Optimized thrust: {sol.value(thrust)} N")
print(f"Optimized wing twist: {sol.value(wing_twist)}")

Result Plotting Code:

# Plotting airfoils and their aerodynamic performance
fig, ax = plt.subplots(figsize=(10, 6))

# Airfoils to plot and their respective colors
airfoils_and_colors = {
    "Dragonfly Airfoil (SD7062)"           : (initial_guess_airfoil, "dimgray"),
    "NeuralFoil-Optimized"    : (optimized_airfoil, "blue"), 
}

# Plot airfoil shapes and aerodynamic polars
for i, (name, (af, color)) in enumerate(airfoils_and_colors.items()):
    color = p.adjust_lightness(color, 1)
    # Plot airfoil shapes
    ax.fill(
        af.x(), af.y(),
        facecolor=(*color, 0.09),
        edgecolor=(*color, 0.6),
        linewidth=1,
        label=name,
        linestyle=(3 * i, (7, 2)),
        zorder=4 if "NeuralFoil" in name else 3,
    )

# Finalize plot layout and show
ax.legend(fontsize=11, loc="lower right", ncol=len(airfoils_and_colors)//2)
ax.set_title("Airfoil Shapes")
ax.set_xlabel("xc")
ax.set_ylabel("yc")
ax.axis('equal')
plt.savefig("airfoil_shape_comparison.png")
p.show_plot("Comparison of NeuralFoil-Optimized-, and Dragonfly Airfoil (SD7062)", legend=False)

# Create an Airfoil object with these coordinates
wing_airfoil = optimized_airfoil

# Define airfoil for the tail
tail_airfoil = asb.Airfoil("naca0009")

# Define the airplane model
airplaneDFY = asb.Airplane(
    name="DBF23 Dragonfly",
    wings=[
        # Main Wing Definition
        asb.Wing(
            name="Main Wing",
            symmetric=True,
            xsecs=[
                # Wing Cross Sections: Root, Mid, and Tip
                asb.WingXSec(xyz_le=[-0.025, 0, 0], chord=0.40, twist=optimized_twist, airfoil=wing_airfoil),
                asb.WingXSec(xyz_le=[-0.025, 0.5, 0], chord=0.40, twist=optimized_twist, airfoil=wing_airfoil),
                asb.WingXSec(xyz_le=[-0.025, 1, 0], chord=0.40, twist=optimized_twist, airfoil=wing_airfoil),
            ]        
        ),
        asb.Wing(
            name="Horizontal Stabilizer",
            symmetric=True,
            xsecs=[
                asb.WingXSec(  # root
                    xyz_le=[0, 0, 0],
                    chord=0.2,
                    twist=0,
                    airfoil=tail_airfoil,
                ),
                asb.WingXSec(  # tip
                    xyz_le=[0, 0.3, 0],
                    chord=0.2,
                    twist=0,
                    airfoil=tail_airfoil
                )
            ]
        ).translate([0.923, 0, 0]),
        asb.Wing(
            name="Vertical Stabilizer",
            symmetric=False,
            xsecs=[
                asb.WingXSec(
                    xyz_le=[0, 0, 0],
                    chord=0.304,
                    twist=0,
                    airfoil=tail_airfoil,
                ),
                asb.WingXSec(
                    xyz_le=[0.131, 0, 0.35],
                    chord=0.173,
                    twist=0,
                    airfoil=tail_airfoil
                )
            ]
        ).translate([0.82, 0, 0])
    ],
    fuselages=[
        # Fuselage Tail Boom Definition
        asb.Fuselage(
            name="FuselageTAILBOOM",
            xsecs=[
                asb.FuselageXSec(
                    xyz_c=[0.0602 * xi, 0, 0],
                    width=0.025,
                    height=0.025,
                    shape=7
                )
                for xi in range(10)
            ]
        ).translate([0.4, 0, -0.0125]),

        # Fuselage Body Definition
        asb.Fuselage(
            name="FuselageBODY",
            xsecs=[
                asb.FuselageXSec(
                    xyz_c=[0.042 * xi, 0, 0],
                    width=0.1,
                    height=0.1,
                    shape=7
                )
                for xi in range(10)
            ]
        ).translate([-0.25, 0, -0.05]),

        # Fuselage Nose
        asb.Fuselage(
            name="FuselageNOSE",
            xsecs=[
                asb.FuselageXSec(
                    xyz_c=[0.22 * xi, 0, 0],
                    width=-0.1 * np.sqrt(xi),  # Adjust the width to be less pointy
                    height=-0.1 * np.sqrt(xi),  # Adjust the height to be less pointy
                    shape=7
                )
                for xi in np.linspace(0, 1, 40)
            ]
        ).translate([-0.47, 0, -0.05]),

        # Transition Section from Fuselage Body to Tail Boom
        asb.Fuselage(
            name="FuselageREAR",
            xsecs=[
                asb.FuselageXSec(
                    # Adjust the center point of each cross-section to slope the bottom surface
                    xyz_c=[0.27 * xi, 0, 0.036 * xi],  # Move the center point up for the slope
                    width=np.interp(xi, [0, 1], [0.1, 0.025]),
                    height=np.interp(xi, [0, 1], [0.1, 0.025]),
                    shape=7
                )
                for xi in np.linspace(0, 1, 40)
            ]
        ).translate([0.13, 0, -0.05]),
        # Antenna modeled with rectangular cross-sections in the x-direction
        asb.Fuselage(
            name="Antenna",
            xsecs=[
                asb.FuselageXSec(
                    xyz_c=[0.075 + xi * antenna_diameter, 1, antenna_length / 2],  # Position in x, y, z
                    width=antenna_diameter * np.sin(np.pi * xi),  # Width of the cross-section
                    height=antenna_length,  # Height of the cross-section
                    shape=1000  # Shape parameter for rectangular shape
                )
                for xi in np.linspace(0, 1, num_sections)
            ]
        ),
    ],
    s_ref= 1.85
)

airplaneDFY.draw_three_view()

plotter = dyn.draw(
    vehicle_model=airplaneDFY,
    scale_vehicle_model=7,  # Draw the vehicle at x scale
    show=False,
    n_vehicles_to_draw=30
)

plotter.camera.elevation = 45 
plotter.camera.azimuth = 45 
plotter.camera.zoom(0.9)
plotter.show(jupyter_backend="static")