Getting elevator deflection in real time, for a optimized trajectory.

[sumitkumarsahoo2017]

How I can get the elevator deflection as a function of time by using DynamicsRigidBody3DBodyEuler for endurance optimization. Can anyone help me with this code I am very new to coding and aerosandbox.

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

g = 3.72

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.813537984559229



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


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=-500,
    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
)

time = np.arange(0, 120, 0.5)


opti = asb.Opti()
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 * u_e, scale=1e3),
        y_e=opti.variable(init_guess=time * v_e, scale=1e3),
        z_e=opti.variable(init_guess=time * w_e, scale=1e3),
        u_b=opti.variable(init_guess=dyn_init.u_b * np.ones_like(time), scale=1e2),
        v_b=opti.variable(init_guess=dyn_init.v_b * np.ones_like(time), scale=1e2),
        w_b=opti.variable(init_guess=dyn_init.w_b * np.ones_like(time), scale=1e2),
        phi=opti.variable(init_guess=dyn_init.phi * np.ones_like(time), scale=1e-1),
        theta=opti.variable(init_guess=dyn_init.theta * np.ones_like(time), scale=1e-1),
        psi=opti.variable(init_guess=dyn_init.psi * np.ones_like(time), scale=1e-1),
        p=opti.variable(init_guess=np.zeros_like(time), scale=1e-2),
        q=opti.variable(init_guess=np.zeros_like(time), scale=1e-2),
        r=opti.variable(init_guess=np.zeros_like(time), scale=1e-2),
        elevator_deflection = opti.variable(init_guess=0, lower_bound=-25, upper_bound=25),
    )
)

elevator_deflection = opti.variable(init_guess=0, lower_bound=-25, upper_bound=25)
airplane = airplane.with_control_deflections(
    control_surface_deflection_mappings={
        "elevator": elevator_deflection,
    }
)
aero = asb.AeroBuildup(
    airplane=airplane,
    op_point=dyn.op_point,
    xyz_ref=dyn.mass_props.xyz_cg
).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"]>0,
    aero["CL"]<1.5,
    dyn.op_point.mach()<0.8,
    dyn.op_point.velocity < 86,
    dyn.theta < np.pi/2,
])  

endurance = ?
opti.maximize(endurance)
sol = opti.solve(max_iter=4000, behavior_on_failure="return_last")
dyn = sol(dyn)

```