Minimize power | ASB + Neuralfoil

[kjayawar]

Hi Peter

Thanks heaps for these awesome libraries. I am new to both Aerosandbox and NeuralFoil.
I tried modifying the wing optimization tutorial by incorparting CL,CD values taken from NeuralFoil
When I try to optimize drag, it works but when I try to optimize for power the optimizer fails.

import aerosandbox as asb
import aerosandbox.numpy as np

from neuralfoil import get_aero_from_kulfan_parameters
from aerosandbox.geometry.airfoil.airfoil_families import get_kulfan_parameters

airfoil = asb.Airfoil("dae11").repanel().normalize()
kulfan_parameters = get_kulfan_parameters(
    airfoil.coordinates,
    n_weights_per_side=8
)

def airfoil_clcd(alpha, re):
	x = get_aero_from_kulfan_parameters(alpha=alpha, Re=re, kulfan_parameters=kulfan_parameters)
	return x['CL'], x['CD']

# Let's define some constants.
density = 1.225  # Note: AeroSandbox uses base SI units (or derived units thereof) everywhere, with no exceptions.
viscosity = 1.81e-5  # So the units here are kg/(m*s).
weight = 100  # 100 N, corresponding to the weight of around 10 kg.
CDAo = 0.04
e = .85

# Time to start constructing our problem!
opti = asb.Opti()

aspect_ratio = opti.variable(init_guess=10)  # as an experiment, let's leave this unbounded.
wing_area = opti.variable(init_guess=1)
alpha = opti.variable(init_guess=4)

# We don't know the best airspeed to fly at, so let's make that an unknown, too.
airspeed = opti.variable(init_guess=30, lower_bound=0) # See the following comment block

span = (aspect_ratio * wing_area) ** 0.5  # Equivalent to saying: AR = b^2 / S
chord = span / aspect_ratio

dynamic_pressure = 0.5 * density * airspeed ** 2  # Definition of dynamic pressure

Re = density * airspeed * chord / viscosity  # Definition of Reynolds number

CL, Cd = airfoil_clcd(alpha, Re)
CD_profile = Cd + CDAo/wing_area
CD_induced = CL ** 2 / (np.pi * e * aspect_ratio)  # Textbook equation for induced drag, assuming optimal lift distribution.
CD = CD_profile + CD_induced
drag = dynamic_pressure * wing_area * CD
lift = CL * dynamic_pressure * wing_area  # This must be true in steady level flight
power = drag * airspeed

opti.subject_to([
    lift == weight
])

opti.minimize(power)

sol = opti.solve()

AR = sol.value(aspect_ratio)
WA = sol.value(wing_area)
AL = sol.value(alpha)
AS = sol.value(airspeed)

print(f"AR = {AR}")
print(f"WA = {WA}  m^2")
print(f"AL = {alpha}")
print(f"AS = {AS}")

I tried changing log_transformations on and off for variables but was not able to get the solution to converge.

RuntimeError: Error in Opti::solve [OptiNode] at .../casadi/core/optistack.cpp:157:
.../casadi/core/optistack_internal.cpp:998: Assertion "return_success(accept_limit)" failed:
Solver failed. You may use opti.debug.value to investigate the latest values of variables. return_status is 'Diverging_Iterates'

Could someone please point me to the right direction as of how this has to be changed in order to get sensible values ?

Thanks
Kanishka.

[peterdsharpe]

Hey Kanishka,

Thanks for writing in, and for including the code snippet!

The "diverging iterates" error means that the problem appears to be unbounded - some variable is being driven towards infinity.

To determine which it is, we can modify our opti.solve() command to instead say:

opti.solve(behavior_on_failure="return_last")

This will make it so that the optimizer will ignore any RuntimeError that is raised (such as the diverging iterates error) and return the last value of the variables anyway. The reason this is not the default behavior (and is instead "opt-in" with the behavior_on_failure keyword argument) is to ensure that users are 100% aware that the result is not a converged solution.

With that line of code changed, we can re-run the code. On my machine, that produces:

AR = 106891.89476263986
WA = 2.4438014646453908e+20  m^2
AL = 22382.562930647404
AS = 5.260276056086263e-07

The issue here is apparent - airspeed is going to zero. Some more digging (use sol(x) to inspect various variables) shows that this is because wing area is going to infinity. Mathematically, this makes sense. If there's no penalty for making wing area huge (which, in the real world, would be weight), then a power-minimization problem wants you to fly at an infinitesimally-slow speed with an infinitely-large wing.

So, the reason AeroSandbox isn't giving a solution is because no optimum (at least, with finite values) exists to the problem as-written if the objective is power minimization - mathematically, the problem is not well-posed. Two possible strategies (among many possible ones) to make the problem well-posed are:

1. Add some kind of weight scaling that accounts for the growing weight of a large wing. Purely as an illustrative example (this equation is totally made up), one could overwrite the weight definition to read:

weight = 100 * wing_area ** 1.5

2. Alternatively, add an upper bound on the wing area (this is the "wrong" approach in that it doesn't capture the design intent, but it will make the problem mathematically well-posed). As another crude example:

wing_area = opti.variable(init_guess=1, upper_bound=10)

---

Unrelated, but you probably also want to add sensible bounds constraints on alpha. NeuralFoil will handle any value of alpha just fine (e.g., if you ask for airfoil aerodynamics at $\alpha = 370\degree$, it's essentially modulo'd back to $\alpha = 10\degree$), but this inherently creates multiple optima - 370 degrees is just as good as 10 degrees. Hence, by adding bounds constraints, the design space is more likely to be unimodal, and the gradient-based optimizer will tend to be happier.

---

Separately, you probably also want to add a CL >= 0 constraint in order to make convergence more reliable. If the optimizer gets itself into a state where CL is negative, the reversed gradient information will just drive the optimizer try to fly slower and slower to restore lift == weight. This is, of course, impossible so long as CL is negative - so it can drive you to weird parts of the design space.

[kjayawar]

Thanks heaps for the great answer Peter. Much appreciated :)