is Opti() not supported for optimization of CST parameters of airfoils

[Zcaic]

Hello, is Aerosandbox.Opti() not supported for optimization of CST parameters of airfoils?
I user opti.varibale to define get_kulfan_coordinates()'s upper_weights and lower_weights,
then use get_aero_from_neuralfoil() to get CD, then opti.minimize(CD),
but it doesn't seem to work.

[Zcaic]

I found ""np.argmin()" doesn't support casadi.

[peterdsharpe]

Hey @Zcaic !

This is an area of active work-in-progress that I am extremely excited for! :)

Previously, with XFoil-driven aerodynamics, optimizing on airfoils through the asb.Opti() interface was not possible (as XFoil was in the loop, which is a black box). Instead, one could use the separate airfoil optimizer code here to "pre-optimize" the airfoils before overall aircraft optimization - but this airfoil optimization was done with non-differentiable (Nelder-Mead simplex via scipy.optimize) optimization that was separate from asb.Opti(), and quite slow.

Now, with NeuralFoil, optimizing CST parameters of airfoils (i.e., airfoil shape optimization as part of a larger aircraft design optimization problem) should in theory be possible!

However, the nuts and bolts of the code haven't yet been fully implemented - I imagine this will happen in the next few weeks.

In the time between, you may have luck with the following steps:

1. Define an asb.Opti() optimization environment.
2. Define CST parameters as optimization variables (ideally, with some sensible lower and upper bounds). Maybe try keeping the TE_thickness parameter set to 0 for the first attempts.
3. Plug those CST parameters directly into neuralfoil.get_aero_from_kulfan_parameters() This function is yet-to-be-documented, but the kulfan_parameters keyword argument is a dictionary of the form described here.
4. Set an optimization objective/constraints (perhaps, go for max L/D?) and solve.

Critically, Step 3 bypasses the conversion to/from a coordinate array, eliminating the argmin needed to identify the LE index.

This approach may need some massaging to work correctly, but in theory it should work. If you give this a shot, please let me know how it goes! I'm extremely excited to see how NeuralFoil does when stress-tested against airfoil aerodynamic shape optimization.

When I have time to take a stab at this and get it working, I will post a Jupyter notebook in AeroSandbox/tutorials/ and update this issue here.

[Zcaic]

it is my test

import aerosandbox as asb
import aerosandbox.numpy as np
from aerosandbox.geometry.airfoil.airfoil_families import get_kulfan_coordinates,get_kulfan_parameters
from neuralfoil import get_aero_from_kulfan_parameters
import logging
logging.basicConfig(level=logging.INFO,format='%(levelname)s - %(asctime)s - %(message)s...')

def track_opt(iter):
    logging.info(f"optimizer state at iteration {iter}")
    track_lp.append(opti.debug.value(lp))
    track_up.append(opti.debug.value(up))

oriairfoil=asb.Airfoil(name='rae2822')
initcst=get_kulfan_parameters(coordinates=oriairfoil.coordinates)

opti = asb.Opti()

lp = opti.variable(init_guess=initcst['lower_weights'], upper_bound=-0.03, lower_bound=-0.5)
up = opti.variable(init_guess=initcst['upper_weights'], upper_bound=0.5, lower_bound=0.05)

mycst = dict(lower_weights=lp, upper_weights=up, TE_thickness=0, leading_edge_weight=0)
rsaero = get_aero_from_kulfan_parameters(kulfan_parameters=mycst, alpha=2, Re=6.98886e+06)

ck = -rsaero["CL"] / rsaero["CD"]
opti.minimize(ck)

track_lp=[]
track_up=[]
sol = opti.solve(callback=track_opt,verbose=False)
print(track_lp[-1])
print(track_up[-1])
print(sol.stats()['iterations']['obj'])

the history of obj is:
[-71.97528142783199, -72.13468982455441, -67.62159392513874, -79.97649494975234, -122.52266687532232, -133.51252039590796, -180.17954523397108, -183.93608137653388, -184.81119028819847, -210.3904160677403, -210.82175351449737, 
-232.7252212728867, -243.22206547147832, -251.63075279570074, -254.38170153061793, -255.48376530501085, -255.55609826991125, -255.59643616947838, -255.5969481213508, -255.59695565203143]

it seems to work sucessfully
then I use optimal solution (lp and up) to verify:

lp = np.array([-0.05640185, -0.03, -0.03, -0.03, -0.03, -0.03, -0.03, -0.03])  # the optimal lp
up = np.array([0.13740123, 0.29148949, 0.29450211, 0.27986522, 0.5, 0.5, 0.17367956, 0.3624671]) # the optimal up

myairfoil=asb.Airfoil(coordinates=get_kulfan_coordinates(lower_weights=lp,upper_weights=up))
rsaero=myairfoil.get_aero_from_neuralfoil(alpha=2,Re=6.98886e+06)
ck=rsaero['CL']/rsaero['CD']
print(f"neuralfoil: {ck}")
xf=asb.XFoil(airfoil=myairfoil,Re=6.98886e+06)
rsaero=xf.alpha(2)
ck=rsaero['CL']/rsaero['CD']
print(f"xfoil: {ck}")

output:
neuralfoil: [255.59443472]
xfoil: [248.25301205]

it works well although the CL/CD has a little difference between neuralfoil and xfoil

[peterdsharpe]

This is so awesome! 😄 <3% difference in L/D between NeuralFoil and XFoil after optimization is excellent to hear, and better than I was expecting. You may try supplying model_size="xxxlarge" to the analysis method call, which may tighten this gap even more.

Great work here!

[Zcaic]

@peterdsharpe hello, Is it possible to use an external optimizer?
for example, I use aerosandbox to get objective and gradient values, and pass to the external gradient optimizer.

[Zcaic]

I use code below to get objective and gradient：

import casadi as ca
import aerosandbox as asb
opti=asb.Opti()
twi=opti.variable()
...
ap=asb.AeroBuildup(
    airplane=myairplane,
    op_point=asb.OperatingPoint(
        velocity=100,
        alpha=twi
    )
).run()
ck=ap['CL']/ap['CD']
f=ca.Function('ck',[twi],[ck])
dfdt=ca.Function('dfdt',[tw0],[ca.gradient(ck,twi)])
print(f(0))
print(dfdt(0))

it can work as I expect.

[peterdsharpe]

Nice work! This is exactly the right way to access the gradient.

This actually reminds me that I should put together an AeroSandbox tutorial on how to connect AeroSandbox to NLOPT, for users who want to try other optimizers (e.g., gradient-free or global). Perhaps even making this a method of asb.Opti would be useful!

[peterdsharpe]

By the way @Zcaic ,

Check out these two brand-new tutorials which show AeroSandbox airfoil optimization through a few examples:

• Airfoil Analysis
• Airfoil Optimization

[Zcaic]

Oh, it's a so amazing work! @peterdsharpe
I have learned these two tutorials Airfoil Analysis and Airfoil Optimization.
For optimization with many parameters, It saves a lot of time than using evolutionary algorithms.

[Zcaic]

But I am curious about what happens about optimization at high Mach. I also tried to optimize at Ma=0.8, Re=1.8e7, the optimizer ran 1000 steps and warn that optimization failed. The optimized airfoil also looks a little strange.

[peterdsharpe]

I'm glad to hear that the new tutorials are helpful!

I have also noticed some weird behavior in terms of transonic airfoil optimization - it seems to behave worse than subsonic airfoil optimization.

I think that this is because of a few reasons:

1. When I was writing the loss function (i.e., objective function) for NeuralFoil's training, I weighted the importance of accurate CL, CD, and CM prediction more heavily than $C_{p,\min}$ prediction. In transonic cases, accurately estimating $C_{p,\min}$ is crucial in order to get the wave drag right. So, this inaccuracy is sometimes driving the optimizer in the wrong direction.
2. I want to embed more physics-based features into NeuralFoil (basically, feature engineering), which I think will make the neural network more accurate and parameter-efficient.
3. The training data for NeuralFoil consists of randomly-generated airfoils (with some steering to bias it towards "typical" airfoils), which are then fed into XFoil to generate training data. In cases of weird airfoils, XFoil will fail to converge, which means that the weird airfoil is not included into the training data. I want to rework the training data for NeuralFoil so that, instead of dropping the sample if XFoil cannot converge, it progressively regularizes the airfoil (and penalizes its performance accordingly) to still include the weird airfoil. This will cause optimization using NeuralFoil to have more "optimization pressure" driving it away from weird airfoil designs.

I won't have time for the next few weeks to work on these improvements (busy with other projects), but I hope to get around to it within the next few months. This is definitely an area where contributions would be welcome!

[peterdsharpe]

(Converting this issue to a discussion, as the initial question has been resolved and now we're exploring other fun topics)

[Zcaic]

I did some validation and I think it might be some problem with xfoil when calculating high Reynolds and high Mach
This is my optimization result using Aerosandbox:

and aerosandbox shows CL/CD=40.15 at Ma=0.8,Re=1.8e7,alpha=4deg, then I export optimized airfoil to xflr5, the CL/CD in xflr5 is about 37, at last I verify with CFD, the CL/CD only about 1.3.
It shows that NeuralFoil get a good result close to xfoil (xflr5), but maybe xfoil has difficulties to calculating high Reynolds and high Mach.

[Zcaic]

@peterdsharpe , hello!
Do you think it is a good way that I want to run DOE airfoil cfd with about tens of thousands of grid cells (each case doesn't take long ) using SU2 and then add the cfd result to NeuralFoil training set.
Because I have a optizimed_airfoil optimized_airfoil.txt from NACA 0012 using Aerosandbox. And I found its CL/CD at Re=1.4e7, Mach=0.6 has a big difference between xflr5 (xfoil, NeuralFoil) and CFD (STARCCM+, Fluent, SU2).
So I want to add CFD results to NeuralFoil.