Skip to content

Commit

Permalink
Updated parameter introspection and traj parameters (#964)
Browse files Browse the repository at this point in the history
* Significant updates to introspection.
_get_targets_metadata now returns metadata on a per target basis, and determining the appropriate shape of a control, state, or
parameter is done outside of the function.
brachistochrone and some other tests are working but a lot of work is left to be done.
this should allow us to introspect the appropriate default shape and value for parameters and other variables.

* removal of some unused arguments

* cleaning up some ambiguity in values for tests

* fixed a bug in handling of static parameter targets during refinement

* Cleanup of a final few tests.

* Fixed an error where shape introspection wasn't happening on control rate2 targets

* removed checks for OpenMDAO older than the oldest supported version

* Better error messages with specific tests.  Allow corner case where parameter has no targets since it may be used as a rate source.

* pep8 fixes

* cleanup

* updates per Bret's comments
  • Loading branch information
robfalck authored Aug 23, 2023
1 parent 1aca42a commit 99ec21e
Show file tree
Hide file tree
Showing 24 changed files with 1,455 additions and 531 deletions.
Original file line number Diff line number Diff line change
Expand Up @@ -75,10 +75,10 @@ def test_conflicting_static_target(self):

expected_msg = "Invalid parameter in phase `traj.phases.phase0`.\n" \
"Parameter `g` has invalid target(s).\n" \
"User has specified 'static_target = False' for parameter g,but one or more " \
"targets is tagged with 'dymos.static_target': g"
"User has specified 'static_target = False' for parameter `g`,\nbut one or more " \
"targets is tagged with 'dymos.static_target':\n{'g'}"

self.assertEqual(str(e.exception), expected_msg)
self.assertEqual(expected_msg, str(e.exception))

def test_control_to_static_target_fails_gl(self):
""" Tests that control cannot be connected to target tagged as 'dymos.static_target'. """
Expand Down

Large diffs are not rendered by default.

Original file line number Diff line number Diff line change
Expand Up @@ -80,8 +80,6 @@ def compute(self, inputs, outputs):
@use_tempdirs
class TestConnectControlToParameter(unittest.TestCase):

@unittest.skipIf(om_version < (3, 4, 1) or (om_version == (3, 4, 1) and om_dev_version),
'test requires OpenMDAO >= 3.4.1')
@require_pyoptsparse(optimizer='SLSQP')
def test_connect_control_to_parameter(self):
""" Test that the final value of a control in one phase can be connected as the value
Expand Down
220 changes: 220 additions & 0 deletions dymos/examples/cannonball/test/test_traj_param_static_and_dynamic.py
Original file line number Diff line number Diff line change
@@ -0,0 +1,220 @@
import unittest

import numpy as np

import openmdao
import openmdao.api as om

from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse

from dymos.examples.cannonball.cannonball_ode import rho_interp

om_dev_version = openmdao.__version__.endswith('dev')
om_version = tuple(int(s) for s in openmdao.__version__.split('-')[0].split('.'))


GRAVITY = 9.80665 # m/s**2


class CannonballODEVectorCD(om.ExplicitComponent):
"""
Cannonball ODE assuming flat earth and accounting for air resistance
"""

def initialize(self):
self.options.declare('num_nodes', types=int)
self.options.declare('static_gravity', types=bool)

def setup(self):
nn = self.options['num_nodes']

# static parameters
self.add_input('m', units='kg')
self.add_input('S', units='m**2')

if self.options['static_gravity']:
self.add_input('g', units='m/s**2', val=GRAVITY)
else:
self.add_input('g', units='m/s**2', val=GRAVITY * np.ones(nn,))

# This will be used as both a control and a parameter
self.add_input('CD', 0.5, shape=nn)

# time varying inputs
self.add_input('h', units='m', shape=nn)
self.add_input('v', units='m/s', shape=nn)
self.add_input('gam', units='rad', shape=nn)

# state rates
self.add_output('v_dot', shape=nn, units='m/s**2')
self.add_output('gam_dot', shape=nn, units='rad/s')
self.add_output('h_dot', shape=nn, units='m/s')
self.add_output('r_dot', shape=nn, units='m/s')
self.add_output('ke', shape=nn, units='J')

# Ask OpenMDAO to compute the partial derivatives using complex-step
# with a partial coloring algorithm for improved performance
self.declare_partials('*', '*', method='cs')
self.declare_coloring(wrt='*', method='cs')

def compute(self, inputs, outputs):

gam = inputs['gam']
v = inputs['v']
h = inputs['h']
m = inputs['m']
S = inputs['S']
CD = inputs['CD']

# handle complex-step gracefully from the interpolant
if np.iscomplexobj(h):
rho = rho_interp(inputs['h'])
else:
rho = rho_interp(inputs['h']).real

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


@use_tempdirs
class TestTrajParamStaticAndDynamic(unittest.TestCase):

@require_pyoptsparse(optimizer='SLSQP')
def test_traj_param_static_and_dynamic(self):
""" Test that the final value of a control in one phase can be connected as the value
of a parameter in a subsequent phase. """
import openmdao.api as om
from openmdao.utils.assert_utils import assert_near_equal

import dymos as dm
from dymos.examples.cannonball.size_comp import CannonballSizeComp

p = om.Problem(model=om.Group())

p.driver = om.pyOptSparseDriver()
p.driver.options['optimizer'] = 'SLSQP'
p.driver.declare_coloring()

external_params = p.model.add_subsystem('external_params', om.IndepVarComp())

external_params.add_output('radius', val=0.10, units='m')
external_params.add_output('dens', val=7.87, units='g/cm**3')

external_params.add_design_var('radius', lower=0.01, upper=0.10, ref0=0.01, ref=0.10)

p.model.add_subsystem('size_comp', CannonballSizeComp())

traj = p.model.add_subsystem('traj', dm.Trajectory())

transcription = dm.Radau(num_segments=5, order=3, compressed=True)
ascent = dm.Phase(ode_class=CannonballODEVectorCD, transcription=transcription,
ode_init_kwargs={'static_gravity': True})

ascent = traj.add_phase('ascent', ascent)

# All initial states except flight path angle are fixed
# Final flight path angle is fixed (we will set it to zero so that the phase ends at apogee)

ascent.set_time_options(fix_initial=True, duration_bounds=(1, 100), duration_ref=100, units='s')
ascent.add_state('r', fix_initial=True, fix_final=False, rate_source='r_dot', units='m')
ascent.add_state('h', fix_initial=True, fix_final=False, units='m', rate_source='h_dot')
ascent.add_state('gam', fix_initial=False, fix_final=True, units='rad', rate_source='gam_dot')
ascent.add_state('v', fix_initial=False, fix_final=False, units='m/s', rate_source='v_dot')

ascent.add_parameter('S', targets=['S'], units='m**2', static_target=True)
ascent.add_parameter('mass', targets=['m'], units='kg', static_target=True)

ascent.add_control('CD', targets=['CD'], opt=False, val=0.05)

# Limit the muzzle energy
ascent.add_boundary_constraint('ke', loc='initial',
upper=400000, lower=0, ref=100000)

# Second Phase (descent)
transcription = dm.GaussLobatto(num_segments=5, order=3, compressed=True)
descent = dm.Phase(ode_class=CannonballODEVectorCD, transcription=transcription,
ode_init_kwargs={'static_gravity': False})

traj.add_phase('descent', descent)

# All initial states and time are free (they will be linked to the final states of ascent.
# Final altitude is fixed (we will set it to zero so that the phase ends at ground impact)
descent.set_time_options(initial_bounds=(.5, 100), duration_bounds=(.5, 100),
duration_ref=100, units='s')
descent.add_state('r', units='m', rate_source='r_dot')
descent.add_state('h', units='m', rate_source='h_dot', fix_initial=False, fix_final=True)
descent.add_state('gam', units='rad', rate_source='gam_dot', fix_initial=False, fix_final=False)
descent.add_state('v', units='m/s', rate_source='v_dot', fix_initial=False, fix_final=False)

descent.add_parameter('S', targets=['S'], units='m**2', static_target=True)
descent.add_parameter('mass', targets=['m'], units='kg', static_target=True)
descent.add_parameter('CD', targets=['CD'], val=0.01)

descent.add_objective('r', loc='final', scaler=-1.0)

# Add externally-provided design parameters to the trajectory.
# In this case, we connect 'm' to pre-existing input parameters named 'mass' in each phase.
traj.add_parameter('m', units='kg', val=1.0,
targets={'ascent': 'mass', 'descent': 'mass'}, static_target=True)

# In this case, by omitting targets, we're connecting these parameters to parameters
# with the same name in each phase.
traj.add_parameter('S', units='m**2', val=0.005, static_target=True)

# Link Phases (link time and all state variables)
traj.link_phases(phases=['ascent', 'descent'], vars=['*'])

# Issue Connections
p.model.connect('external_params.radius', 'size_comp.radius')
p.model.connect('external_params.dens', 'size_comp.dens')

p.model.connect('size_comp.mass', 'traj.parameters:m')
p.model.connect('size_comp.S', 'traj.parameters:S')

traj.connect('ascent.timeseries.CD', 'descent.parameters:CD', src_indices=[-1], flat_src_indices=True)

# A linear solver at the top level can improve performance.
p.model.linear_solver = om.DirectSolver()

# Finish Problem Setup
p.setup()

# Set Initial Guesses
p.set_val('external_params.radius', 0.05, units='m')
p.set_val('external_params.dens', 7.87, units='g/cm**3')

p.set_val('traj.ascent.controls:CD', 0.5)

p.set_val('traj.ascent.t_initial', 0.0)
p.set_val('traj.ascent.t_duration', 10.0)

p.set_val('traj.ascent.states:r', ascent.interp('r', [0, 100]))
p.set_val('traj.ascent.states:h', ascent.interp('h', [0, 100]))
p.set_val('traj.ascent.states:v', ascent.interp('v', [200, 150]))
p.set_val('traj.ascent.states:gam', ascent.interp('gam', [25, 0]), units='deg')

p.set_val('traj.descent.t_initial', 10.0)
p.set_val('traj.descent.t_duration', 10.0)

p.set_val('traj.descent.states:r', descent.interp('r', [100, 200]))
p.set_val('traj.descent.states:h', descent.interp('h', [100, 0]))
p.set_val('traj.descent.states:v', descent.interp('v', [150, 200]))
p.set_val('traj.descent.states:gam', descent.interp('gam', [0, -45]), units='deg')

dm.run_problem(p, simulate=True, make_plots=True)

assert_near_equal(p.get_val('traj.descent.states:r')[-1], 3183.25, tolerance=1.0E-2)
assert_near_equal(p.get_val('traj.ascent.timeseries.CD')[-1],
p.get_val('traj.descent.parameter_vals:CD')[0])


if __name__ == '__main__': # pragma: no cover
unittest.main()
14 changes: 7 additions & 7 deletions dymos/grid_refinement/grid_refinement_ode_system.py
Original file line number Diff line number Diff line change
Expand Up @@ -3,7 +3,7 @@

from ..phase.options import TimeOptionsDictionary
from ..utils.misc import get_rate_units
from ..utils.introspection import get_targets
from ..utils.introspection import get_targets, _get_targets_metadata
from ..transcriptions.grid_data import GridData


Expand Down Expand Up @@ -189,16 +189,16 @@ def configure(self):

# Configure the parameters
for name, options in self.options['parameters'].items():
static_target = options['static_target']
static_targets = options['static_targets']
shape = options['shape']
prom_name = f'parameters:{name}'
targets = get_targets(self.ode, name, options['targets'])
for tgt in targets:
if not static_target:
targets = _get_targets_metadata(self.ode, name, options['targets'])
for tgt, meta in targets.items():
if tgt in static_targets:
self.promotes('ode', inputs=[(tgt, prom_name)])
else:
self.promotes('ode', inputs=[(tgt, prom_name)],
src_indices=om.slicer[np.zeros(num_nodes, dtype=int), ...])
else:
self.promotes('ode', inputs=[(tgt, prom_name)])
if targets:
self.set_input_defaults(name=prom_name,
src_shape=shape,
Expand Down
1 change: 0 additions & 1 deletion dymos/grid_refinement/hp_adaptive/hp_adaptive.py
Original file line number Diff line number Diff line change
Expand Up @@ -102,7 +102,6 @@ def refine_first_iter(self, refine_results):
dict
A dictionary of phase paths : phases which were refined.
"""

for phase_path, phase_refinement_results in refine_results.items():
phase = self.phases[phase_path]
tx = phase.options['transcription']
Expand Down
16 changes: 13 additions & 3 deletions dymos/phase/options.py
Original file line number Diff line number Diff line change
Expand Up @@ -263,15 +263,25 @@ def __init__(self, read_only=False):
"option 'dynamic' set to False should now use 'static_target' set "
"to True.")

self.declare(name='static_target', values=[True, False, _unspecified], default=_unspecified,
self.declare(name='static_target', default=_unspecified,
desc='True if the target of this parameter does NOT have a unique value at '
'each node in the ODE.'
'If _unspecified, attempt to determine through introspection.')
'If _unspecified, attempt to determine through introspection.',
deprecation='Use option `static_targets` to specify whether all targets\n'
'are static (static_targegts=True), none are static (static_targets=False),\n'
'static_targets are determined via introspection (static_targets=_unspecified),\n'
'or give an explicit sequence of the static targets.')

self.declare(name='static_targets', default=_unspecified,
desc='If a boolean, specifies whether all targets are static (True), or no\n'
'targets are static (False). Otherwise, provide a list of the static\n'
'targets within the ODE. If left unspecified, static targets will be\n'
'determined by finding inptus tagged with \'dymos.static_target\'.')

self.declare(name='targets', allow_none=True, default=_unspecified,
desc='Targets in the ODE to which the state is connected')

self.declare(name='val', types=(Iterable, np.ndarray, Number), default=np.zeros(1),
self.declare(name='val', types=(Iterable, np.ndarray, Number), default=0.0,
desc='The default value of the parameter in the phase.')

self.declare(name='shape', check_valid=check_valid_shape, default=_unspecified,
Expand Down
Loading

0 comments on commit 99ec21e

Please sign in to comment.