Scaling Fuselage Geometry during Optimization

[joyo44]

Hi Peter and ASB Users,

I'm trying to scale a Fuselage during optimization and got an error "Implicit conversion of symbolic CasADi type to numeric matrix not supported".

Here is the code segment:

FuselageLength = np.abs(x_nose) + x_tail

afFuseScaled = asb.Airfoil(
    name=Fuselage_Airfoil,
    coordinates=np.stack((afFuse.x()*FuselageLength, afFuse.y()*Fuselage_Diameter), axis=1)
    )        
    
fuselage=asb.Fuselage(
            name="Fuselage",
            color='b',
            xsecs=[
                asb.FuselageXSec(
                    xyz_c=[x_nose + FuselageLength * xi, 0.0, np.tan(-Fuselage_Ifuse/rad) * xi + Fuselage_Zpos],
                    radius=(0.5 * afFuseScaled.local_thickness(x_over_c=xi))
                )
                for xi in np.cosspace(0, 1, Fuselage_Xpanel_resolution)]
            )

and here is the related error message:

Traceback (most recent call last):
  File "C:\Users\LENOVO\anaconda3\envs\asb3\Lib\site-packages\casadi\casadi.py", line 16407, in __array__
    return self.full()
           ^^^^^^^^^
  File "C:\Users\LENOVO\anaconda3\envs\asb3\Lib\site-packages\casadi\casadi.py", line 13990, in <lambda>
    __getattr__ = lambda self, name: _swig_getattr(self, MX, name)
                                     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\LENOVO\anaconda3\envs\asb3\Lib\site-packages\casadi\casadi.py", line 83, in _swig_getattr
    raise AttributeError("'%s' object has no attribute '%s'" % (class_type.__name__, name))
AttributeError: 'MX' object has no attribute 'full'

During handling of the above exception, another exception occurred:

Traceback (most recent call last):
  File "E:\Home_T420\APAwAX_V550\ASB_Optimization_PY\scripts\ASB_FlyingWingOptimization_Neu_WL.py", line 1869, in <module>
    xsecs=[
          ^
  File "E:\Home_T420\APAwAX_V550\ASB_Optimization_PY\scripts\ASB_FlyingWingOptimization_Neu_WL.py", line 1872, in <listcomp>
    radius=(0.5 * afFuseScaled.local_thickness(x_over_c=xi))
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\LENOVO\anaconda3\envs\asb3\Lib\site-packages\aerosandbox\geometry\airfoil\airfoil.py", line 821, in local_thickness
    upper = self.upper_coordinates()[::-1]
            ^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\LENOVO\anaconda3\envs\asb3\Lib\site-packages\aerosandbox\geometry\airfoil\airfoil.py", line 974, in upper_coordinates
    return self.coordinates[:self.LE_index() + 1, :]
                             ^^^^^^^^^^^^^^^
  File "C:\Users\LENOVO\anaconda3\envs\asb3\Lib\site-packages\aerosandbox\geometry\airfoil\airfoil.py", line 952, in LE_index
    return int(np.argmin(self.x()))
               ^^^^^^^^^^^^^^^^^^^
  File "<__array_function__ internals>", line 200, in argmin
  File "C:\Users\LENOVO\anaconda3\envs\asb3\Lib\site-packages\numpy\core\fromnumeric.py", line 1338, in argmin
    return _wrapfunc(a, 'argmin', axis=axis, out=out, **kwds)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\LENOVO\anaconda3\envs\asb3\Lib\site-packages\numpy\core\fromnumeric.py", line 54, in _wrapfunc
    return _wrapit(obj, method, *args, **kwds)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\LENOVO\anaconda3\envs\asb3\Lib\site-packages\numpy\core\fromnumeric.py", line 43, in _wrapit
    result = getattr(asarray(obj), method)(*args, **kwds)
                     ^^^^^^^^^^^^
  File "C:\Users\LENOVO\anaconda3\envs\asb3\Lib\site-packages\casadi\casadi.py", line 16415, in __array__
    raise Exception("Implicit conversion of symbolic CasADi type to numeric matrix not supported.\n"
Exception: Implicit conversion of symbolic CasADi type to numeric matrix not supported.
This may occur when you pass a CasADi object to a numpy function.
Use an equivalent CasADi function instead of that numpy function.

I tried the direct scale method of the Airfoil Class but that failed also.

Is there mistake in my approach or is there a workaround for this problem?

Thanks for any help.

John

[peterdsharpe]

Hi,

This issue is occuring because the coordinates of afFuseScaled are the downstream result of an optimization variable, not a NumPy array. If you do this, certain operations on that Airfoil are no longer possible to automatic-differentiate through. One such disallowed operation is evaluating the local_thickness of that airfoil.

To understand why, consider:

1. To evaluate the thickness of the airfoil, we first need to identify which points form the upper surface, and which points form the lower surface. To do this, we need to identify the index of the leading edge in the coordinates array.
2. A coordinates array has no explicit information about which index is the leading edge, other than by looking at the point locations themselves (i.e., looking at values).
3. The problem with looking at values is that now we're doing value-dependent code execution, which is a no-go with many automatic differentiation libraries (e.g., JAX, CasADi), since the computational graph is no longer static.

To solve this, you should either:
a) Make sure Airfoil.coordinates is not a downstream result of an optimization variable, or
b) Use the asb.KulfanAirfoil type, which incorporates intrinsic explicit information about where the leading edge is.

---

Separately, it's best practice to keep airfoils scaled in normalized coordinates, rather than rescaling them into physical coordinates (as the code snippet above does). Ultimately, this is the real root cause of this issue, presuming that afFuse has static coordinates. In this case, I think you can fix the issue just by moving that scaling to physical coordinates outside of the Airfoil class instantiation.

[joyo44]

Hi Peter,

thanks for your quick reply. Now I understand the problem.

John

[cdhainaut]

Hi,

The np.argmin function is causing the issue. Like Peter said, argmin is not something that is compatible in an automatic-differentiation context.

If you are optimizing a section with a given topology the LE_index shouldn't change. So you can override the function with an integer value you should know in advance.

I already had similar issue in the past and solved it that way

Bests

[peterdsharpe]

^ This is also a clever solution, and should work!

Just be careful with this, since if you later change the airfoil and forget about this workaround, overwriting LE_index will cause it to silently fail.

(Maybe put some big all-caps warning comments around the place in the code where the airfoil is instantiated?)

[cdhainaut]

^ This is also a clever solution, and should work!

Just be careful with this, since if you later change the airfoil and forget about this workaround, overwriting LE_index will cause it to silently fail.

(Maybe put some big all-caps warning comments around the place in the code where the airfoil is instantiated?)

Yes true... This could be an interesting option to add to the airfoil class actually. Weather to automatically find it or impose a certain index

[joyo44]

Hi,

I tried KulfanAirfoil following the description:

        # From Class KulfanAirfoil:
        # -------------------------
        # if scale_x != 1:
            # raise ValueError("\n".join([
                # "Scaling a KulfanAirfoil in the x-direction is not supported due to inherent limitations of the",
                # "Kulfan parameterization. If you need to scale in the x-direction: "
                # "\t- Convert to an Airfoil first (`KulfanAirfoil.to_airfoil()`)"
                # "\t- Scale it (`Airfoil.scale()`)"
                # "\t- Convert back to a KulfanAirfoil (`Airfoil.to_kulfan_airfoil()`)"
            # ]))


        FuselageLength = np.abs(x_nose) + x_tail    # opti.variables

        afFuse       = asb.KulfanAirfoil(name=Fuselage_Airfoil).normalize()
        FuseAirfoil  = afFuse.to_airfoil().scale(scale_x=FuselageLength, scale_y=Fuselage_Diameter)
        afFuseScaled = FuseAirfoil.to_kulfan_airfoil()
        
        fuselage=asb.Fuselage(
                    name="Fuselage",
                    color='b',
                    xsecs=[
                        asb.FuselageXSec(
                            xyz_c=[x_nose + FuselageLength * xi, 0.0, np.tan(-Fuselage_Ifuse/rad) * xi + Fuselage_Zpos],
                            radius=(0.5 * afFuseScaled.local_thickness(x_over_c=xi))
                        )
                        for xi in np.cosspace(0, 1, Fuselage_Xpanel_resolution)]
                    )

but again got a casadi error :

Traceback (most recent call last):
  File "E:\Home_T420\APAwAX_V550\ASB_Optimization_PY\scripts\ASB_FlyingWingOptimization_Neu_WL.py", line 1863, in <module>
    FuseAirfoil = afFuse.to_airfoil().scale(scale_x=FuselageLength, scale_y=1.0)
                  ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "C:\Users\LENOVO\anaconda3\envs\asb3\Lib\site-packages\aerosandbox\geometry\airfoil\airfoil.py", line 1423, in scale
    if scale_x < 0:
       ^^^^^^^^^^^
  File "C:\Users\LENOVO\anaconda3\envs\asb3\Lib\site-packages\casadi\casadi.py", line 14029, in __bool__
    return _casadi.MX___bool__(self, *args)
           ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
RuntimeError: .../casadi/core/mx_node.cpp:182: Can only determine truth value of a numeric MX.

I'm still struggling with the casadi types.

[cdhainaut]

When you use to_airfoil() function, you convert a KulfanAirfoil to an Airfoil object. And in the Airfoil scaling function, you have a boolean condition (line 1423) on "scale_x" which happens to be a symbolic variable in your case.

If/Else conditions is only meaningful for numeric variables. Not for symbolic expressions that are built from CasADi variables.

You need to see writing code with CasADi as if you had no way to guess numeric values of the variables you are using... The only things you can do from CasADi variables is using basic maths functions (addition, multiplication, log, sum, exp, cos, etc...).