Cannonball Implicit Duration - Dymos
Cannonball Implicit Duration - Dymos
Cannonball Implicit Duration - Dymos
This example demonstrates determining the duration of a phase using a nonlinear solver to satisfy an endpoint condition. As in
the multi-phase cannonball problem, we will simultaneously find the optimal design for the cannonball (its radius) and the
optimal flight profile (its launch angle). However, here we will use a single phase that terminates at the condition
h(tf ) = 0 (91)
This component sits upstream of the trajectory model and feeds its outputs to the trajectory as parameters.
import numpy as np
from scipy.interpolate import interp1d
import matplotlib.pyplot as plt
import openmdao.api as om
import dymos as dm
from dymos.models.atmosphere.atmos_1976 import USatm1976Data
#############################################
# Component for the design part of the model
#############################################
class CannonballSizeComp(om.ExplicitComponent):
"""
Compute the area and mass of a cannonball with a given radius and density.
Notes
-----
This component is not vectorized with 'num_nodes' as is the usual way
with Dymos, but is instead intended to compute a scalar mass and reference
area from scalar radius and density inputs. This component does not reside
in the ODE but instead its outputs are connected to the trajectory via
input design parameters.
"""
def setup(self):
self.add_input(name='radius', val=1.0, desc='cannonball radius', units='m')
self.add_input(name='dens', val=7870., desc='cannonball density', units='kg/m**3')
self.declare_partials(of='mass', wrt='dens')
self.declare_partials(of='mass', wrt='radius')
self.declare_partials(of='S', wrt='radius')
class CannonballODE(om.ExplicitComponent):
"""
Cannonball ODE assuming flat earth and accounting for air resistance
"""
def initialize(self):
self.options.declare('num_nodes', types=int)
def setup(self):
nn = self.options['num_nodes']
# static parameters
self.add_input('m', units='kg')
self.add_input('S', units='m**2')
# 0.5 good assumption for a sphere
self.add_input('CD', 0.5)
# state rates
self.add_output('v_dot', shape=nn, units='m/s**2', tags=['dymos.state_rate_source:v'])
self.add_output('gam_dot', shape=nn, units='rad/s', tags=['dymos.state_rate_source:gam'])
self.add_output('h_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:h'])
self.add_output('r_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:r'])
self.add_output('ke', shape=nn, units='J')
gam = inputs['gam']
v = inputs['v']
h = inputs['h']
m = inputs['m']
S = inputs['S']
CD = inputs['CD']
q = 0.5*rho*inputs['v']**2
qS = q * S
D = qS * CD
cgam = np.cos(gam)
sgam = np.sin(gam)
outputs['v_dot'] = - D/m-GRAVITY*sgam
outputs['gam_dot'] = -(GRAVITY/v)*cgam
outputs['h_dot'] = v*sgam
outputs['r_dot'] = v*cgam
outputs['ke'] = 0.5*m*v**2
return guess
p = om.Problem(name='cannonball_implicit_duration', model=om.Group())
p.driver = om.pyOptSparseDriver()
p.driver.options['optimizer'] = 'IPOPT'
p.driver.declare_coloring()
p.driver.opt_settings['derivative_test'] = 'first-order'
p.driver.opt_settings['mu_strategy'] = 'monotone'
p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'
p.driver.opt_settings['bound_mult_init_method'] = 'mu-based'
p.driver.opt_settings['mu_init'] = 0.01
p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'
p.set_solver_print(level=0, depth=99)
p.model.add_subsystem('size_comp', CannonballSizeComp(),
promotes_inputs=['radius', 'dens'])
p.model.set_input_defaults('dens', val=7.87, units='g/cm**3')
p.model.add_design_var('radius', lower=0.01, upper=0.10,
ref0=0.01, ref=0.10, units='m')
phase.add_boundary_constraint('ke', loc='initial',
upper=400000, lower=0, ref=100000)
# Prior to setup, this phase has the default NonlinearRunOnce and LinearRunOnce solvers.
# Because of the implicit behavior introduced by set_duration_balance, it will automatically
# use a NewtonSolver with an Armijo Goldstein linesearch.
# In this case that linesearch causes extrapolation of the atmospheric model into invalid regions.
# To account for this, we explicitly assign a NewtonSolver with a different linesearch.
phase.nonlinear_solver = om.NewtonSolver(iprint=0, solve_subsystems=True)
phase.nonlinear_solver.linesearch = None
p.model.connect('size_comp.mass', 'traj.phase.parameters:m')
p.model.connect('size_comp.S', 'traj.phase.parameters:S')
#############################################
# Set constants and initial guesses
#############################################
p.set_val('radius', 0.05, units='m')
p.set_val('dens', 7.87, units='g/cm**3')
p.set_val('traj.phase.parameters:CD', 0.5)
guess = initial_guess(t_dur=30.0)
p.set_val('traj.phase.t_initial', 0.0)
p.set_val('traj.phase.t_duration', 30.0)
#####################################################
# Run the optimization and final explicit simulation
#####################################################
dm.run_problem(p, simulate=True, make_plots=True)
/usr/share/miniconda/envs/test/lib/python3.10/site-packages/dymos/transcriptions/transcription_base.py:448:
OpenMDAOWarning:'traj.phases.phase' <class Phase>
Non-default solvers are required
implicit duration: True
solved segments: True
input initial: False
Setting `traj.phases.phase.linear_solver = om.DirectSolver(iprint=2)`
Explicitly set traj.phases.phase.linear_solver to override.
Set `traj.phases.phase.options["auto_solvers"] = False` to disable this behavior.
Full total jacobian for problem 'cannonball_implicit_duration' was computed 3 times, taking 0.2783104929999354
seconds.
Total jacobian shape: (98, 100)
Solution:
--------------------------------------------------------------------------------
Total Time: 17.3794
User Objective Time : 14.6487
User Sensitivity Time : 2.4103
Interface Time : 0.1811
Opt Solver Time: 0.1392
Calls to Objective Function : 239
Calls to Sens Function : 55
Objectives
Index Name Value
0 traj.phase.states:r -3.182962E+03
Exit Status
Inform Description
0 Solve Succeeded
--------------------------------------------------------------------------------
False
Results of cannonball_implicit_duration
Creation Date: 2024-03-29 20:22:47
Plot phases:
phase