Cannonball Example with Euler Integration and an External Optimizer#
This example will show you how to do the following:
Use an OpenMDAO Problem inside of a loop to perform a broader calculation (in this case, simple integration.)
Optimize a Problem using an external optimizer with a functional interface.
Perform a complex step across an OpenMDAO Problem.
In the example, we want to find the optimal angle to fire a cannon to maximize the distance it travels down range. We already know that a 45 degree angle will maximize the range in ideal conditions, but the model we will use also includes aerodynamic effects, so we will solve for the firing angle in the presence of a small amount of drag force.
Model#
A very general set of dynamics for a simplified aircraft are given by these differential equations:
A very general set of dynamics for a simplified aircraft are given by these differential equations:
We will use these for the cannonball dynamics, though angle of attack, lift, and thrust will all be
zero. We implement these equations into an OpenMDAO ExplicitComponent whose outputs are the
time rates of change of the states which are velocity ‘v’, flight path angle ‘gam’, altitude
‘h’ and range ‘r’.
import numpy as np
import openmdao.api as om
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')
def setup_partials(self):
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
We also need to compute the aerodynamic forces L and D, which are computed from the coefficients of lift and drag as well as the dynamic pressure. This was implemented in two components:
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
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
Finally, we put it all together in our top model. Given any cannonball flight path angle, and velocity we can compute the rates of change of all states by running this model.
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')
Time Integration#
Given an initial angle and velocity for the cannonball, one way we can compute the range is by integrating the equations of motion that are provided by OpenMDAO. A simple way to do this is to use the Euler integration. We can perform this by choosing a time step, running the Problem at the starting location to compute the state rates, and then use Euler’s method to compute the new cannonball state. We can do this sequentially until the height turns negative, which means we have hit the ground sometime between this time step and the previous one. Finally, we use linear interpolation between the final point and the previous one to find the location where the height is zero.
Note that Euler integration is a first-order method, so the accuracy will be linearly proportional to the step size. We use a ‘dt’ of 0.1 seconds here, which isn’t particularly accurate, but runs quickly. In practice, you might use a higher order integration method here (e.g., Runge-Kutta.)
import numpy as np
from scipy.optimize import minimize
def eval_cannonball_range(gam_init, prob, complex_step=False):
"""
Compute distance given initial speed and angle of cannonball.
Parameters
----------
gam_init : float
Initial cannonball firing angle in degrees.
prob : <Problem>
OpenMDAO problem that contains the equations of motion.
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(f"Distance: {r_final}, Time: {t}, Angle: {gam_init}")
return -r_final
Here, we have placed the integration code inside of a function that takes the initial angle as an argument and returns the negative of the computed range. This is the objective we wish to optimize with an external optimizer which requires a function that it can call to evaluate the objective.
Providing Derivatives for an External Optimizer#
The optimizer also allows you to specify a function to evaluate the gradient. If you do not provide one, it will use finite difference. We can improve the accuracy by performing complex step instead. OpenMDAO allows you to run a model in complex mode. When the mode is enabled on the Problem, you can use ‘set_val’ to set complex values and ‘get_val’ to retrieve them. In the code for the objective evaluation above, we turn this feature on by calling prob.set_complex_step_mode(True). Likewise, it is important to turn it off when not needed, and it should only be used when you are performing a complex step to compute the derivatives.
The following function computes the total derivatives that the external optimizer needs.
def gradient_cannonball_range(gam_init, prob):
"""
Uses complex step to compute gradient of range wrt initial angle.
Parameters
----------
gam_init : float
Initial cannonball firing angle in degrees.
prob : <Problem>
OpenMDAO problem that contains the equations of motion.
Returns
-------
float
Derivative of range wrt initial angle in m/deg.
"""
step = 1.0e-14
dr_dgam = eval_cannonball_range(gam_init + step * 1j, prob, complex_step=True)
return dr_dgam.imag / step
Running the Optimization#
Now we can put everything together and run the optimization. Our optimizer is scipy.minimize for this example.
import numpy as np
from scipy.optimize import minimize
def eval_cannonball_range(gam_init, prob, complex_step=False):
"""
Compute distance given initial speed and angle of cannonball.
Parameters
----------
gam_init : float
Initial cannonball firing angle in degrees.
prob : <Problem>
OpenMDAO problem that contains the equations of motion.
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(f"Distance: {r_final}, Time: {t}, Angle: {gam_init}")
return -r_final
def gradient_cannonball_range(gam_init, prob):
"""
Uses complex step to compute gradient of range wrt initial angle.
Parameters
----------
gam_init : float
Initial cannonball firing angle in degrees.
prob : <Problem>
OpenMDAO problem that contains the equations of motion.
Returns
-------
float
Derivative of range wrt initial angle in m/deg.
"""
step = 1.0e-14
dr_dgam = eval_cannonball_range(gam_init + step * 1j, prob, complex_step=True)
return dr_dgam.imag / step
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.05) # 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)
result = minimize(eval_cannonball_range, 27.0,
method='SLSQP',
jac=gradient_cannonball_range,
args=(prob))
print(result['x'])
[42.3810579]