adding wind as a function of altitude in simulate.py, where altitude is a decision variable

[sumitkumarsahoo2017]

The code below is giving an error as follows when I am trying to add wind to our consideration which varies with z_e(which is a decision variable). Is there any solution to this problem please share, I will be extremely grateful.

--> 100     for i, h in enumerate(height):
    101 
    102         # Find the index of the nearest lower height
    103         index = 0

[~/miniconda3/envs/SuEnv/lib/python3.12/site-packages/casadi/casadi.py](https://file+.vscode-resource.vscode-cdn.net/d%3A/8TH%20SEM/Martian%20Glider%20Codes/~/miniconda3/envs/SuEnv/lib/python3.12/site-packages/casadi/casadi.py) in ?(self)
  16298     def __iter__(self):
> 16299       raise Exception("""CasADi matrices are not iterable by design.
...

Exception: CasADi matrices are not iterable by design.
                      Did you mean to iterate over m.nz, with m IM/DM/SX?
                      Did you mean to iterate over horzsplit(m,1)/vertsplit(m,1) with m IM/DM/SX/MX?

import aerosandbox as asb
import aerosandbox.numpy as np
import casadi as cas

m_fuselage = 4  # aircraft weight excluding wing [N]

g = 3.72
# opti = asb.Opti()  # Initialize an optimization environment
span = 2.2976561585848687
x2 = 0.1958611108267831
z2 = 0.08523498840079038
chord1 = 0.24964754416884002
chord2 = 0.15366729678999783
twist2 = -5.543260585617984
alpha = 8.016486709052963
speed = 132.2275150337934


translate_tail_x = 0.23427720512140934
translate_tail_y = 1.1188253995687496
vtail_chord_root = 0.1999974438396932
vtail_chord_tip = 0.14999735568673864
vtail_x_tip = 0.1741087959237856
vtail_z_tip = 0.3783308841244814
vtail_twist = 0.5461087414168971
fuselage_width = 0.15
fuselage_nose_x = 0.15128645451938888
fuselage_twist = -0.5105524759927914

x_cg = 0.1178286638578197
mass_glider = 10.81353798455922


Ixx=1.600312206
Iyy=0.133854194
Izz=1.659027376
Ixy=0
Iyz=0
Ixz=0.041791945
mass_properties_of_glider=asb.MassProperties(
            mass=mass_glider,
            x_cg=x_cg, y_cg=0, z_cg=0,
            Ixx=1.600312206, Iyy=0.133854194, Izz=1.659027376,
            Ixy=0, Iyz=0, Ixz=0.041791945,
        )

y_aileron_outboard = translate_tail_y - 0.05
y_aileron = 0.9522717861734794

deflection_aileron = 0
hinge_point_aileron = 0.7

phi = np.pi/4
deflection_rudder = 0
hinge_point_rudder = 0.88

def get_x(y_le):
    return 0+(x2)*(y_le/(span/2))
def get_z(y_le):
    return 0+(z2)*(y_le/(span/2))
def get_twist(y_le):
    return 0+(twist2)*(y_le/(span/2))
def get_chord(y_le):
    return chord1+(chord2-chord1)*(y_le/(span/2))



wing_airfoil = asb.Airfoil("a63a108c")
tail_airfoil = asb.Airfoil("a63a108c")
fuselage_airfoil = asb.Airfoil("sc20012")

wing_1 = asb.Wing(
            name="Main Wing",
            symmetric=True, 
            xsecs=[  
                asb.WingXSec( 
                    xyz_le=[0, 0 , 0], 
                    chord=chord1,
                    twist=0,  # degrees
                    airfoil=wing_airfoil, 
                    control_surfaces=[
                            asb.ControlSurface(
                                name="flap",
                                symmetric=True,
                                hinge_point=0.75,
                                deflection=0,
                            )
                        ]
                
                ),
                asb.WingXSec(
                        xyz_le=[get_x(y_aileron), y_aileron, get_z(y_aileron)],
                        chord=get_chord(y_aileron),
                        twist=get_twist(y_aileron),
                        airfoil=tail_airfoil,
                        control_surfaces=[
                            asb.ControlSurface(
                                name="aileron",
                                symmetric=False,
                                hinge_point=hinge_point_aileron,
                                deflection=deflection_aileron,
                            ),
                            asb.ControlSurface(
                                name="elevator",
                                symmetric=True,
                                hinge_point=hinge_point_aileron,
                                deflection=deflection_aileron,
                            )
                        ]
                    ),

                asb.WingXSec(
                        xyz_le=[get_x(y_aileron_outboard), y_aileron_outboard , get_z(y_aileron_outboard)],
                        chord=get_chord(y_aileron_outboard),
                        twist=get_twist(y_aileron_outboard),
                        airfoil=tail_airfoil,
                    ),
                    
                asb.WingXSec(  # Tip
                    xyz_le=[x2, span / 2 , z2],
                    chord=chord2,
                    twist=twist2,
                    airfoil=wing_airfoil,
                ),
            ]
        )


vtail = asb.Wing(
            name="Vertical Stabilizer",
            symmetric=True,
            xsecs=[
                asb.WingXSec(
                    xyz_le=[0, 0, 0],
                    chord=vtail_chord_root,
                    twist=vtail_twist,
                    airfoil=tail_airfoil,
                    control_surfaces=[
                            asb.ControlSurface(
                                name="rudder",
                                symmetric=True,
                                hinge_point=hinge_point_rudder,
                                deflection=deflection_rudder,
                            )
                        ]
                 
                ),
             
                asb.WingXSec(
                    xyz_le=[vtail_x_tip, 0, vtail_z_tip],
                    chord=vtail_chord_tip,
                    twist=vtail_twist,
                    airfoil=tail_airfoil,
                )
            ]
        )

fuselage_of_glider =  asb.Wing(
            name="fuselage_wing",
            symmetric=True,
            xsecs=[
                asb.WingXSec(
                    xyz_le=[-fuselage_nose_x*(np.cosd(yi)**2), fuselage_width*yi/90, 0],
                    chord=0.24964754416884002+(fuselage_nose_x*(np.cosd(yi)**2)),
                    twist=fuselage_twist,
                    airfoil=fuselage_airfoil,
                )
                for yi in np.linspace(0, 60, 10)
            
            ] + [
                asb.WingXSec(  # Root
                    xyz_le=[0.0, fuselage_width, 0],  # Coordinates of the XSec's leading edge, relative to the wing's leading edge.
                    chord=0.24964754416884002,
                    twist=0,  # degrees
                    airfoil=wing_airfoil
                ),
            ]
        )

xyz_cg = [
    x_cg, 
    0, 
    0
]
airplane = asb.Airplane(
    name="Martian Glider",
    xyz_ref=xyz_cg,  # CG location
    
    wings=[
        wing_1.translate([0, fuselage_width, 0]),
        vtail.translate([translate_tail_x, translate_tail_y+fuselage_width, z2]),
        # vtail,
        fuselage_of_glider,

    ],   
)

alpha_init = alpha
V_init = speed

dyn_init = asb.DynamicsRigidBody3DBodyEuler(
    mass_props=mass_properties_of_glider,
    x_e=0,
    y_e=0,
    z_e=-50,
    u_b=V_init * np.cosd(alpha_init),
    v_b=0,
    w_b=V_init * np.sind(alpha_init),
    phi=0,
    theta=np.radians(alpha_init),
    psi=0,
    p=0,
    q=0,
    r=0
)



opti = asb.Opti()
t_max_init_guess = 200
t_max = opti.variable(init_guess=t_max_init_guess)
time = np.linspace(0, t_max, 100)
time_init = np.linspace(0, t_max_init_guess, 100)

elevator_deflection = opti.variable(init_guess=np.ones_like(time)*(-5), lower_bound=-25, upper_bound=25)
# rudder_deflection = opti.variable(init_guess=np.ones_like(time)*(0), lower_bound=-25, upper_bound=25)
# aileron_deflection = opti.variable(init_guess=np.ones_like(time)*(0), lower_bound=-25, upper_bound=25)

u_e, v_e, w_e = dyn_init.convert_axes(
    dyn_init.u_b, dyn_init.v_b, dyn_init.w_b,
    from_axes="body",
    to_axes="earth"
)


dyn = dyn_init.get_new_instance_with_state(
    dict(
        x_e=opti.variable(init_guess=time_init * u_e, scale=1e3),
        y_e=opti.variable(init_guess=time_init * v_e),
        z_e=opti.variable(init_guess=dyn_init.z_e + time_init * w_e, scale=1e2,lower_bound=-1000 ,upper_bound=-10),
        u_b=opti.variable(init_guess=dyn_init.u_b * np.ones_like(time_init), scale=1e2),
        v_b=opti.variable(init_guess=dyn_init.v_b * np.ones_like(time_init)),
        w_b=opti.variable(init_guess=dyn_init.w_b * np.ones_like(time_init)),
        phi=opti.variable(init_guess=dyn_init.phi * np.ones_like(time_init)),
        theta=opti.variable(init_guess=dyn_init.theta * np.ones_like(time_init)),
        psi=opti.variable(init_guess=dyn_init.psi * np.ones_like(time_init)),
        p=opti.variable(init_guess=np.zeros_like(time_init)),
        q=opti.variable(init_guess=np.zeros_like(time_init)),
        r=opti.variable(init_guess=np.zeros_like(time_init)),
    )
)

def mars_wind_profile(height):
    """
    Computes the wind speed at given heights in the Martian atmosphere.

    Args:
        height: Array-like object containing heights above the surface of Mars, in meters.

    Returns:
        wind_speed: Array of wind speeds at the given heights, in meters per second.
    """
    # Define the wind profile data
    heights = [0, 170,  522.8, 1500]  # Heights corresponding to WH0, WH1, WH2 in meters
    wind_speeds = [5.8, 20.1, 26.0, 28.63844]  # Wind speeds at corresponding heights in m/s
    gradients = [0.08411765, 0.016723356, 0.0027]  # Gradients (dW/dz) at corresponding heights in sec^-1

    # Initialize an empty array to store interpolated wind speeds
    interpolated_wind_speeds = np.zeros_like(height)

    # Linearly interpolate between nearest heights
    for i, h in enumerate(height):
       
        # Find the index of the nearest lower height
        index = 0
        while h > heights[index + 1]:
            index += 1

        # Linear interpolation
        height_lower = heights[index]
        height_upper = heights[index + 1]
        wind_speed_lower = wind_speeds[index]
        wind_speed_upper = wind_speeds[index + 1]
        gradient = gradients[index]

        # Interpolated wind speed using linear approximation
        interpolated_wind_speeds[i] = wind_speed_lower + (h - height_lower) * gradient

    return interpolated_wind_speeds


height = -dyn.z_e

wind_speed = mars_wind_profile(height)
wind_x_b, wind_y_b, wind_z_b = dyn.convert_axes(
    wind_speed, 0, 0,
    from_axes="earth",
    to_axes="body"
)
dyn_for_aero = dyn.get_new_instance_with_state(dict(
    u_b = dyn.u_b + wind_x_b,
    v_b = dyn.v_b + wind_y_b,
    w_b = dyn.w_b + wind_z_b,
))



airplane = airplane.with_control_deflections(
    control_surface_deflection_mappings={
        "elevator": elevator_deflection,
    }
)

llt = asb.AeroBuildup(
    airplane=airplane,
    op_point=dyn_for_aero.op_point,
    xyz_ref=xyz_cg,
)

aero=llt.run_with_stability_derivatives()

dyn.add_gravity_force()
dyn.add_force(*aero["F_b"], axes="body")
dyn.add_moment(*aero["M_b"], axes="body")
dyn.constrain_derivatives(
    opti=opti,
    time=time,
)



opti.subject_to([
    v[0] == v_init
    for v, v_init in zip(dyn.state.values(), dyn_init.state.values())
])
opti.subject_to([
    aero["CL"]>-1,
    aero["CL"]<1.5,
    dyn.op_point.mach()<0.8,
    dyn.op_point.velocity > 86,
    dyn.theta < np.pi/2,
    dyn.phi < np.pi/2,
    t_max > 0

])  

state_derivatives = dyn.state_derivatives()
objective = ((dyn.kinetic_energy[-1]+dyn.potential_energy[-1])-(dyn.kinetic_energy[0]+dyn.potential_energy[0]))
opti.maximize(objective)
sol = opti.solve(max_iter=4000, behavior_on_failure="return_last")

dyn = sol(dyn)

print("elevator_deflection", sol.value(elevator_deflection))
print("dyn.op_point", sol.value(dyn.op_point))
print("energy gain", ((dyn.kinetic_energy[-1]+dyn.potential_energy[-1])-(dyn.kinetic_energy[0]+dyn.potential_energy[0]))/(dyn.kinetic_energy[0]+dyn.potential_energy[0]))


aero = asb.AeroBuildup(
    airplane=airplane,
    op_point=dyn.op_point,
    xyz_ref=dyn.mass_props.xyz_cg,
).run_with_stability_derivatives(alpha=True, beta= True, p=True, q=True, r=True)

import matplotlib.pyplot as plt
import aerosandbox.tools.pretty_plots as p
import matplotlib

vars_to_plot = {
    **dyn.state,
    "alpha"   : dyn.alpha,
    "beta"    : dyn.beta,
    "speed"   : dyn.speed,
    "altitude": dyn.altitude,
    "CL"      : aero["CL"],
    "CD"      : aero["CD"],
    "Cm"      : aero["Cm"],
}


fig, ax = plt.subplots()
plt.plot(dyn.x_e, dyn.altitude, "k")
sc = plt.scatter(dyn.x_e, dyn.altitude, c=dyn.speed, cmap=plt.cm.get_cmap("rainbow"), zorder=4)
plt.axis('equal')
plt.colorbar(sc, label="Airspeed [m/s]")  # Use the returned object from plt.scatter to create the colorbar
p.show_plot("Trajectory using `asb.AeroBuildup` Flight Dynamics", "$x_e$", "$-z_e$")

unique_speeds, counts = np.unique(dyn.speed, return_counts=True)
print("Unique speeds:", unique_speeds)
print("Counts:", counts)

[peterdsharpe]

Hi - the immediate answer is that you can iterate by direct indexing:

for i in range(len(height):
    h = height[i]

But, this is likely not what you want to do here. The Mars atmosphere model you're trying to implement is not $C^1$-continuous. So, if this is implemented as-written, the resulting problem will likely fail to solve.

1. For a longer explanation of why this is, see details in this tutorial.
2. Read through the asb.Atmosphere implementation to see how a differentiable atmosphere can be modeled. In particular, see here.

Separately, direct indexing like this is discouraged because it leads to deep, slow-to-trace computational graphs in the automatic differentiation engine. Instead, try to vectorize or use vector conditionals when able.

[sumitkumarsahoo2017]

"Hello, Peter Sir🙏. I wanted to thank you for your guidance on modifying the wind_speed function. I made the changes you suggested and it was very helpful. Here is the updated code:

from aerosandbox.modeling.interpolation import InterpolatedModel

altitude_knot_points = np.linspace(0, 10e3, 1000)

def wind_speed_calculation(HeightLower, WindSpeedLower, Gradient, h):
    return WindSpeedLower + (h - HeightLower) * Gradient

def wind_speed(height):
    wind = np.zeros_like(height)  # Initialize an array to store wind speed
    heights = [0, 170, 522.8, 5000]  # Heights corresponding to WH0, WH1, WH2 in meters
    wind_speeds = [5.8, 20.1, 26.0, 38.08844]  # Wind speeds at corresponding heights in m/s
    gradients = [0.08411765, 0.016723356, 0.0027, 0.0027]  # Gradients (dW/dz) at corresponding heights in sec^-1
    for i in range(len(heights)):
        wind = np.where(
            height > heights[i],
            wind_speed_calculation(
                HeightLower=heights[i],
                WindSpeedLower=wind_speeds[i],
                h=height,
                Gradient=gradients[i]
            ),
            wind
        )
    return wind

# Calculate wind speed at knot points
wind_knot_points = wind_speed(altitude_knot_points)

# Create interpolated model for wind speed
interpolated_wind = InterpolatedModel(
    x_data_coordinates=altitude_knot_points,
    y_data_structured=wind_knot_points,
)

def wind_differentiable(height):
    return interpolated_wind(height)

height = -dyn.z_e
wind_speed = wind_differentiable(height)

[sumitkumarsahoo2017]

@peterdsharpe could you please confirm if the simulation is correctly accounting for the effect of wind?

dyn_init = asb.DynamicsRigidBody3DBodyEuler(
    mass_props=mass_properties_of_glider,
    x_e=0,
    y_e=0,
    z_e=-1000,
    u_b=V_init * np.cosd(alpha_init),
    v_b=0,
    w_b=V_init * np.sind(alpha_init),
    phi=0,
    theta=np.radians(alpha_init),
    psi=0,
    p=0,
    q=0,
    r=0
)



opti = asb.Opti()
t_max_init_guess = 20
t_max = opti.variable(init_guess=t_max_init_guess, lower_bound=0, upper_bound=600)
# time = np.arange(0, t_max, 1)
# t_max= 10
time = np.linspace(0, t_max, 20)
time_init = np.linspace(0, t_max_init_guess, 20)

elevator_deflection = opti.variable(init_guess=np.ones_like(time)*(9), lower_bound=-25, upper_bound=25)
rudder_deflection = opti.variable(init_guess=np.ones_like(time)*(0), lower_bound=-25, upper_bound=25)
aileron_deflection = opti.variable(init_guess=np.ones_like(time)*(0), lower_bound=-25, upper_bound=25)

u_e, v_e, w_e = dyn_init.convert_axes(
    dyn_init.u_b, dyn_init.v_b, dyn_init.w_b,
    from_axes="body",
    to_axes="earth"
)
dyn = dyn_init.get_new_instance_with_state(
    dict(
        x_e=opti.variable(init_guess=time_init * u_e, scale=1e3),
        y_e=opti.variable(init_guess=time_init * v_e),
        z_e=opti.variable(init_guess=dyn_init.z_e + time_init * w_e, scale=1e2,lower_bound=-5000 ,upper_bound=-10),
        u_b=opti.variable(init_guess=dyn_init.u_b * np.ones_like(time_init), scale=1e2),
        v_b=opti.variable(init_guess=dyn_init.v_b * np.ones_like(time_init)),
        w_b=opti.variable(init_guess=dyn_init.w_b * np.ones_like(time_init)),
        phi=opti.variable(init_guess=dyn_init.phi * np.ones_like(time_init)),
        theta=opti.variable(init_guess=dyn_init.theta * np.ones_like(time_init)),
        psi=opti.variable(init_guess=dyn_init.psi * np.ones_like(time_init)),
        p=opti.variable(init_guess=np.zeros_like(time_init)),
        q=opti.variable(init_guess=np.zeros_like(time_init)),
        r=opti.variable(init_guess=np.zeros_like(time_init)),
    )
)

height = -dyn.z_e

wind_speed = wind_differentiable(height)


wind_x_b, wind_y_b, wind_z_b = dyn.convert_axes(
    wind_speed, 0, 0,
    from_axes="earth",
    to_axes="body"
)
dyn_for_aero = dyn.get_new_instance_with_state(dict(
    u_b = dyn.u_b + wind_x_b,
    v_b = dyn.v_b + wind_y_b,
    w_b = dyn.w_b + wind_z_b,
))



airplane = airplane.with_control_deflections(
    control_surface_deflection_mappings={
        "elevator": elevator_deflection,
        "rudder" : rudder_deflection,
        "aileron" : aileron_deflection
    }
)

llt = asb.AeroBuildup(
    airplane=airplane,
    op_point=dyn_for_aero.op_point,
    xyz_ref=xyz_cg,
)

aero=llt.run_with_stability_derivatives()

dyn.add_gravity_force()
dyn.add_force(*aero["F_b"], axes="body")
dyn.add_moment(*aero["M_b"], axes="body")
dyn.constrain_derivatives(
    opti=opti,
    time=time,
)




opti.subject_to([
    v[0] == v_init
    for v, v_init in zip(dyn.state.values(), dyn_init.state.values())
])
opti.subject_to([
    aero["CL"]>-1,
    aero["CL"]<1.2,
    dyn.op_point.mach()<0.8,
    dyn.op_point.velocity > 86,
    dyn.theta < np.pi/2,
    dyn.phi < np.pi/2,
    t_max>=20,
])  

state_derivatives = dyn.state_derivatives()

# Access the value corresponding to a particular key

objective = ((dyn.kinetic_energy[-1]+dyn.potential_energy[-1])-(dyn.kinetic_energy[0]+dyn.potential_energy[0]))/(dyn.kinetic_energy[0]+dyn.potential_energy[0])*100
opti.maximize(objective)
sol = opti.solve(max_iter=4000, behavior_on_failure="return_last")

dyn = sol(dyn)