Skip to content

Instantly share code, notes, and snippets.

@Kenneth-T-Moore

Kenneth-T-Moore/gist:a4df46c7e5491bf6e1b5ab8d9e4c2006

Created November 13, 2020 22:26
Show Gist options
  • Save Kenneth-T-Moore/a4df46c7e5491bf6e1b5ab8d9e4c2006 to your computer and use it in GitHub Desktop.
Save Kenneth-T-Moore/a4df46c7e5491bf6e1b5ab8d9e4c2006 to your computer and use it in GitHub Desktop.
Download ZIP
Example: complex step to compute derivs around euler integration of openmdao problem.
Raw
gistfile1.txt
import numpy as np
from scipy.optimize import minimize
import openmdao.api as om
class DynamicPressureComp(om.ExplicitComponent):
def setup(self):
self.add_input(name='rho', val=1.0, units='kg/m**3',
desc='atmospheric density')
self.add_input(name='v', val=1.0, units='m/s',
desc='air-relative velocity')
self.add_output(name='q', val=1.0, units='N/m**2',
desc='dynamic pressure')
self.declare_partials(of='q', wrt='rho')
self.declare_partials(of='q', wrt='v')
def compute(self, inputs, outputs):
outputs['q'] = 0.5 * inputs['rho'] * inputs['v'] ** 2
def compute_partials(self, inputs, partials):
partials['q', 'rho'] = 0.5 * inputs['v'] ** 2
partials['q', 'v'] = inputs['rho'] * inputs['v']
class LiftDragForceComp(om.ExplicitComponent):
"""
Compute the aerodynamic forces on the vehicle in the wind axis frame
(lift, drag, cross) force.
"""
def initialize(self):
self.options.declare('num_nodes', types=int)
def setup(self):
self.add_input(name='CL', val=0.0,
desc='lift coefficient')
self.add_input(name='CD', val=0.0,
desc='drag coefficient')
self.add_input(name='q', val=0.0, units='N/m**2',
desc='dynamic pressure')
self.add_input(name='S', val=0.0, units='m**2',
desc='aerodynamic reference area')
self.add_output(name='f_lift', shape=(1, ), units='N',
desc='aerodynamic lift force')
self.add_output(name='f_drag', shape=(1, ), units='N',
desc='aerodynamic drag force')
self.declare_partials(of='f_lift', wrt=['q', 'S', 'CL'])
self.declare_partials(of='f_drag', wrt=['q', 'S', 'CD'])
def compute(self, inputs, outputs):
q = inputs['q']
S = inputs['S']
CL = inputs['CL']
CD = inputs['CD']
qS = q * S
outputs['f_lift'] = qS * CL
outputs['f_drag'] = qS * CD
def compute_partials(self, inputs, partials):
q = inputs['q']
S = inputs['S']
CL = inputs['CL']
CD = inputs['CD']
qS = q * S
partials['f_lift', 'q'] = S * CL
partials['f_lift', 'S'] = q * CL
partials['f_lift', 'CL'] = qS
partials['f_drag', 'q'] = S * CD
partials['f_drag', 'S'] = q * CD
partials['f_drag', 'CD'] = qS
class FlightPathEOM2D(om.ExplicitComponent):
"""
Computes the position and velocity equations of motion using a 2D flight path
parameterization of states per equations 4.42 - 4.46 of _[1].
References
----------
.. [1] Bryson, Arthur Earl. Dynamic optimization. Vol. 1. Prentice Hall, p.172, 1999.
"""
def initialize(self):
self.options.declare('num_nodes', types=int)
def setup(self):
self.add_input(name='m', val=1.0, units='kg',
desc='aircraft mass')
self.add_input(name='v', val=1.0, units='m/s',
desc='aircraft velocity magnitude')
self.add_input(name='T', val=0.0, units='N',
desc='thrust')
self.add_input(name='alpha', val=0.0, units='rad',
desc='angle of attack')
self.add_input(name='L', val=0.0, units='N',
desc='lift force')
self.add_input(name='D', val=0.0, units='N',
desc='drag force')
self.add_input(name='gam', val=0.0, units='rad',
desc='flight path angle')
self.add_output(name='v_dot', val=0.0, units='m/s**2',
desc='rate of change of velocity magnitude')
self.add_output(name='gam_dot', val=0.0, units='rad/s',
desc='rate of change of flight path angle')
self.add_output(name='h_dot', val=0.0, units='m/s',
desc='rate of change of altitude')
self.add_output(name='r_dot', val=0.0, units='m/s',
desc='rate of change of range')
self.declare_partials('v_dot', ['T', 'D', 'm', 'gam', 'alpha'])
self.declare_partials('gam_dot', ['T', 'L', 'm', 'gam', 'alpha', 'v'])
self.declare_partials(['h_dot', 'r_dot'], ['gam', 'v'])
def compute(self, inputs, outputs):
g = 9.80665
m = inputs['m']
v = inputs['v']
T = inputs['T']
L = inputs['L']
D = inputs['D']
gam = inputs['gam']
alpha = inputs['alpha']
calpha = np.cos(alpha)
salpha = np.sin(alpha)
cgam = np.cos(gam)
sgam = np.sin(gam)
mv = m * v
outputs['v_dot'] = (T * calpha - D) / m - g * sgam
outputs['gam_dot'] = (T * salpha + L) / mv - (g / v) * cgam
outputs['h_dot'] = v * sgam
outputs['r_dot'] = v * cgam
def compute_partials(self, inputs, partials):
g = 9.80665
m = inputs['m']
v = inputs['v']
T = inputs['T']
L = inputs['L']
D = inputs['D']
gam = inputs['gam']
alpha = inputs['alpha']
calpha = np.cos(alpha)
salpha = np.sin(alpha)
cgam = np.cos(gam)
sgam = np.sin(gam)
mv = m * v
partials['v_dot', 'T'] = calpha / m
partials['v_dot', 'D'] = -1.0 / m
partials['v_dot', 'm'] = (D - T * calpha) / (m**2)
partials['v_dot', 'gam'] = -g * cgam
partials['v_dot', 'alpha'] = -T * salpha / m
partials['gam_dot', 'T'] = salpha / mv
partials['gam_dot', 'L'] = 1.0 / mv
partials['gam_dot', 'm'] = -(L + T * salpha) / (m * mv)
partials['gam_dot', 'gam'] = g * sgam / v
partials['gam_dot', 'alpha'] = T * calpha / mv
partials['gam_dot', 'v'] = g * cgam / v**2 - (L + T * salpha) / (v * mv)
partials['h_dot', 'gam'] = v * cgam
partials['h_dot', 'v'] = sgam
partials['r_dot', 'gam'] = -v * sgam
partials['r_dot', 'v'] = cgam
class CannonballODE(om.Group):
def setup(self):
self.add_subsystem(name='dynamic_pressure',
subsys=DynamicPressureComp(),
promotes=['*'])
self.add_subsystem(name='aero',
subsys=LiftDragForceComp(),
promotes_inputs=['*'])
self.add_subsystem(name='eom',
subsys=FlightPathEOM2D(),
promotes=['*'])
self.connect('aero.f_drag', 'D')
self.connect('aero.f_lift', 'L')
prob = om.Problem(model=CannonballODE())
prob.setup(force_alloc_complex=True)
# Set constants
prob.set_val('CL', 0.0) # Lift Coefficient
prob.set_val('CD', 0.0) # Drag Coefficient
prob.set_val('S', 0.25 * np.pi, units='ft**2') # Wetted Area (1 ft diameter ball)
prob.set_val('rho', 1.225) # Atmospheric Density
prob.set_val('m', 5.5) # Cannonball Mass
prob.set_val('alpha', 0.0) # Angle of Attack (Not Applicable)
prob.set_val('T', 0.0) # Thrust (Not Applicable)
def eval_cannonball_range(gam_init, complex_step=False):
"""
Compute distance given initial speed and angle of cannonball.
Parameters
----------
gam_init : float
Initial cannonball firing angle in degrees.
complex_step : bool
Set to True to perform complex step.
Returns
-------
float
Negative of range in m.
"""
dt = 0.1 # Time step
h_init = 1.0 # Height of cannon.
v_init = 100.0 # Initial cannonball velocity.
h_target = 0.0 #
v = v_init
gam = gam_init
h = h_init
r = 0.0
t = 0.0
if complex_step:
prob.set_complex_step_mode(True)
while h > h_target:
# Set values
prob.set_val('v', v)
prob.set_val('gam', gam, units='deg')
# Run the model
prob.run_model()
# Extract rates
v_dot = prob.get_val('v_dot')
gam_dot = prob.get_val('gam_dot', units='deg/s')
h_dot = prob.get_val('h_dot')
r_dot = prob.get_val('r_dot')
h_last = h
r_last = r
# Euler Integration
v = v + dt * v_dot
gam = gam + dt * gam_dot
h = h + dt * h_dot
r = r + dt * r_dot
t += dt
#print(v, gam, h, r)
# Linear interpolation between last two points to get the landing point accurate.
r_final = r_last + (r - r_last) * h_last / (h_last - h)
if complex_step:
prob.set_complex_step_mode(False)
print(r_final, gam_init)
return -r_final
def gradient_cannonball_range(gam_init):
"""
Uses complex step to compute gradient of range wrt initial angle.
Parameters
----------
gam_init : float
Initial cannonball firing angle in degrees.
Returns
-------
float
Derivative of range wrt initial angle in m/deg.
"""
step = 1.0e-12
dr_dgam = eval_cannonball_range(gam_init + step * 1j, complex_step=True)
return dr_dgam.imag / step
#print(eval_cannonball_range(45.0))
#print(gradient_cannonball_range(40.0))
result = minimize(eval_cannonball_range, 27.0,
method='SLSQP',
jac=gradient_cannonball_range)
print(result)
print('done')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment