NumPy 1.25+ breaks something with NeuralFoil Optimization

[peterdsharpe]

This code runs fine with NumPy 1.24.4, but does not work with NumPy 1.25.0. Will update as investigation continues.

import aerosandbox as asb
import aerosandbox.numpy as np

CL_multipoint_targets = np.array([0.8, 1.0, 1.2, 1.4, 1.5, 1.6])
CL_multipoint_weights = np.array([5, 6, 7, 8, 9, 10])

Re = 500e3 * (CL_multipoint_targets / 1.25) ** -0.5
mach = 0.03

initial_guess_airfoil = asb.KulfanAirfoil("naca0012")
initial_guess_airfoil.name = "Initial Guess (NACA0012)"

opti = asb.Opti()

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,
)

alpha = opti.variable(
    init_guess=np.degrees(CL_multipoint_targets / (2 * np.pi)),
    lower_bound=-5,
    upper_bound=18
)

aero = optimized_airfoil.get_aero_from_neuralfoil(
    alpha=alpha,
    Re=Re,
    mach=mach,
)

opti.subject_to([
    aero["CL"] == CL_multipoint_targets,
    np.diff(alpha) > 0,
    aero["CM"] >= -0.133,
    optimized_airfoil.local_thickness(x_over_c=0.33) >= 0.128,
    optimized_airfoil.local_thickness(x_over_c=0.90) >= 0.014,
    optimized_airfoil.TE_angle() >= 6.03, # Modified from Drela's 6.25 to match DAE-11 case
    optimized_airfoil.lower_weights[0] < -0.05,
    optimized_airfoil.upper_weights[0] > 0.05,
])

opti.subject_to(
    optimized_airfoil.local_thickness() > 0
)

get_wiggliness = lambda af: sum([
    np.sum(np.diff(np.diff(array)) ** 2)
    for array in [af.lower_weights, af.upper_weights]
])

opti.subject_to(
    get_wiggliness(optimized_airfoil) < 2 * get_wiggliness(initial_guess_airfoil)
)

opti.minimize(np.mean(aero["CD"] * CL_multipoint_weights))

sol = opti.solve(
    behavior_on_failure="return_last"
)

optimized_airfoil = sol(optimized_airfoil)
aero = sol(aero)

Update: seems to not be an issue on Linux. Maybe Windows-only?

[peterdsharpe]

Update: resolved. Discrepancy was caused by updates to the NumPy backend from 1.24 -> 1.25, causing slight differences in floating-point operation results. This gets magnified throughout NeuralFoil's neural network, causing a tiny but measurable (on the order of 10 * floating-point epsilon, so around 10^-15) difference in NeuralFoil's outputs (lift, drag). This particular optimization case study happens to be very sensitive to initial conditions, causing it not to converge. Essentially, NeuralFoil's model landscape has a small "wiggle" near the optimum, and this floating-point math difference is enough to trip the optimizer into an infeasible local minimum. For example, changing the initial guess from a NACA0012 to a NACA0013 airfoil resolves the issue.

Fundamentally this is is a one-off issue caused by NeuralFoil's model landscape, not an issue with AeroSandbox, so closing this issue. However, it is useful to know that NumPy 1.24 and 1.25 produce computational results that, while within acceptable floating-point tolerance, are ever-so-slightly different. This may affect repeatability across versions. On the NeuralFoil note, a new version with improved model architecture will be released soon that will improve robustness to initial guesses.

By-product of this is that we can re-enable NumPy 1.25+, allowing AeroSandbox to be installed on Python 3.12+.