Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Updated parameter introspection and traj parameters #964

Merged
merged 11 commits into from
Aug 23, 2023
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