VLM for trajectory optimization

[cdhainaut]

Hi,

I am trying to use vortex lattice methods for optimal control problems. I successfully used the AeroBuildup module for optimizing a trajectory, but optimizing the design with this module doesn't make sense because my design isn't to adapted to it. When trying to switch to vortex lattice method with a number of variables > 1 using the following code:

import aerosandbox.numpy as np
import aerosandbox as asb

opti = asb.Opti()

n_vars = 2
alpha = opti.variable(n_vars=n_vars)
beta = np.linspace(0,2,n_vars)

vlm_results = asb.VortexLatticeMethod(
    airplane=airplane,
    op_point=asb.OperatingPoint(
        atmosphere=asb.Atmosphere(),
        alpha = alpha,
        beta = beta
    )
).run()

opti.subject_to([
    vlm_results["F_g"][1] = 100
])

sol = opti.solve()

I get this error:

Traceback (most recent call last):
  File "AeroSandbox/aerosandbox/aerodynamics/aero_3D/vortex_lattice_method.py", line 197, in run
    steady_freestream_velocity = self.op_point.compute_freestream_velocity_geometry_axes()  # Direction the wind is GOING TO, in geometry axes coordinates
  File "AeroSandbox/aerosandbox/performance/operating_point.py", line 235, in compute_freestream_velocity_geometry_axes
    return self.compute_freestream_direction_geometry_axes() * self.velocity
  File "AeroSandbox/aerosandbox/performance/operating_point.py", line 231, in compute_freestream_direction_geometry_axes
    return self.compute_rotation_matrix_wind_to_geometry() @ np.array([-1, 0, 0])
  File "AeroSandbox/aerosandbox/performance/operating_point.py", line 212, in compute_rotation_matrix_wind_to_geometry
    alpha_rotation = np.rotation_matrix_3D(
  File "AeroSandbox/aerosandbox/numpy/rotations.py", line 104, in rotation_matrix_3D
    return array(rot)
  File "/AeroSandbox/aerosandbox/numpy/array.py", line 31, in array
    return _cas.vertcat(
  File "/home/charles/miniconda3/envs/py39/lib/python3.9/site-packages/casadi/casadi.py", line 466, in vertcat
    def vertcat(*args): return _vertcat(args)
  File "/home/charles/miniconda3/envs/py39/lib/python3.9/site-packages/casadi/casadi.py", line 20204, in _vertcat
    return _casadi._vertcat(*args)
NotImplementedError: Wrong number or type of arguments for overloaded function '_vertcat'.
  Possible prototypes are:
    _vertcat([Sparsity])
    _vertcat([DM])
    _vertcat([SX])
    _vertcat([MX])
  You have: '(([MX|int],DM,[MX|int]))'

I guess this compute_rotation_matrix_wind_to_geometry does not support mixed types. Is it a feature that is planned to be available any time soon ?

Thank you

[peterdsharpe]

Hey there,

This is a good question! Apologies for the slow response.

Essentially, the reason why is that asb.VortexLatticeMethod is not yet vectorized across multiple simultaneous OperatingPoints.

In short, there are two quite-simple issues here that prevent this from happening as-is:

1. The first issue is just a practical one, that doesn't preclude this, but makes it slow: the VLM will need to solve a linear system at each collocation point (say, 200 of them). We could naively implement this as 200 linear solves, but this is wasteful, since the A-matricies (aerodynamic influence coefficients) are identical, and only the b-vectors (essentially the freestream flow through each panel) is changing. So, it would be much faster (almost 200x faster) to do an LU-factorization, save the factorization, and then back-substitute for each b-vector. We can do this, but we'd have to rewrite the linear solve portion of the VLM module (around 30 lines of code edits, probably).

2. The second one is that asb.VortexLatticeMethod still seems to use a few older functions, like compute_rotation_matrix_wind_to_geometry - which is not vectorized, since that would require vectorizing across matrices, and CasADi doesn't efficiently support >2D arrays. But, we're in luck here! In the time since last editing asb.VortexLatticeMethod, we've written a function .convert_axes() which is a method of both asb.OperatingPoint objects and asb.Dynamics_... objects. This function is a matrix-free, fully-vectorized aerospace axis conversion method, which means it would be a drop-in replacement here of that old rotation-matrix function. This would definitely be something good to do (a pure upgrade, as it adds vectorization and costs almost nothing), and not very much work (around 10 lines of code edits) - but my plate is a bit full at the moment.

So, tl;dr is that this would be possible with not too much work! Tackle both these things (~40 lines of code), and asb.VortexLatticeMethod would be vectorized across operating points. I will try to keep this on my radar; in the meantime, PRs welcome here!

[cdhainaut]

Hi Peter, Thank you for answering. For your second point, are you really sure its only a matter of 10 code lines ? Looking at the code, freestream_velocity consists in a Nx3 matrix, where N is the number of collocation points.

freestream_velocities = np.add(
    wide(steady_freestream_velocity) , rotation_freestream_velocities
)
# Nx3, represents the freestream velocity at each panel collocation point (c)

As CasADi MX not support n > 2 dimensions, how would you shape the matrix so that it accounts for this third "time" dimension (or the number of steps of the trajectory) ?

The only quick workaround I see to solve this issue is to consider constant spanwise freestream velocity, so that the first dimension of the array is the number of steps of the simulation/optimization. But this can be limiting in some cases. What do you think ?

Last, I am happy to work on it and share it as a PR when I'll do it. I am just interested in your point of view on this

Thank you

[peterdsharpe]

Easiest way here would be to split up the x, y, and z components here. Generally, you want to split across a) axes with known size (3), and b) small axes, as the loop overhead is minimized. The x/y/z axis meets both these criteria. So you would have:

freestream_velocities_x = ...
freestream_velocities_y = ...
freestream_velocities_z = ...

and then:

freestream_dot_normal = (
    freestream_velocities_x * normal_x +
    freestream_velocities_y * normal_y +
    freestream_velocities_z * normal_z +
)

or, to do the exact same thing a shorthand, more extensible way:

freestream_velocities = [
    freestream_velocities_x,
    freestream_velocities_y,
    freestream_velocities_z
]

normal = [
    normal_x,
    ...
]

import aerosandbox.numpy as np

freestream_dot_normal = np.dot(freestream_velocities, normal, manual=True)

where the manual=True argument will expand and perform the dot product across the list.

---

In retrospect, you're right - it definitely would be more than a 10-line code edit! Perhaps I was too optimistic haha :)

It's doable, just with a bit more work than I had anticipated! Maybe a 200-line code edit?