From 99ec21ee168071fcc1bfdfab9916ae037c83509f Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Wed, 23 Aug 2023 18:20:06 -0400 Subject: [PATCH] Updated parameter introspection and traj parameters (#964) * 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 --- .../test_brachistochrone_static_gravity.py | 6 +- .../test_phase_param_introspection_failure.py | 689 ++++++++++++++++++ .../test/test_connect_control_to_parameter.py | 2 - .../test_traj_param_static_and_dynamic.py | 220 ++++++ .../grid_refinement_ode_system.py | 14 +- .../hp_adaptive/hp_adaptive.py | 1 - dymos/phase/options.py | 16 +- dymos/phase/phase.py | 51 +- .../test/test_input_parameter_connections.py | 18 +- dymos/phase/test/test_simulate.py | 64 +- dymos/test/test_load_case.py | 1 - dymos/test/test_run_problem.py | 1 - dymos/trajectory/trajectory.py | 113 ++- dymos/transcriptions/analytic/analytic.py | 23 +- dymos/transcriptions/common/parameter_comp.py | 9 +- .../explicit_shooting/explicit_shooting.py | 44 +- .../explicit_shooting/ode_evaluation_group.py | 61 +- .../pseudospectral/gauss_lobatto.py | 91 +-- .../pseudospectral/radau_pseudospectral.py | 28 +- dymos/transcriptions/solve_ivp/solve_ivp.py | 2 +- dymos/transcriptions/transcription_base.py | 14 +- dymos/utils/introspection.py | 487 +++++-------- dymos/utils/misc.py | 1 + .../linkage/test/test_linkage_report.py | 30 +- 24 files changed, 1455 insertions(+), 531 deletions(-) create mode 100644 dymos/examples/brachistochrone/test/test_phase_param_introspection_failure.py create mode 100644 dymos/examples/cannonball/test/test_traj_param_static_and_dynamic.py diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_static_gravity.py b/dymos/examples/brachistochrone/test/test_brachistochrone_static_gravity.py index 7ab79e967..e02913d22 100644 --- a/dymos/examples/brachistochrone/test/test_brachistochrone_static_gravity.py +++ b/dymos/examples/brachistochrone/test/test_brachistochrone_static_gravity.py @@ -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'. """ diff --git a/dymos/examples/brachistochrone/test/test_phase_param_introspection_failure.py b/dymos/examples/brachistochrone/test/test_phase_param_introspection_failure.py new file mode 100644 index 000000000..970fa8012 --- /dev/null +++ b/dymos/examples/brachistochrone/test/test_phase_param_introspection_failure.py @@ -0,0 +1,689 @@ +import unittest +import numpy as np +import openmdao.api as om +import dymos as dm + +from openmdao.utils.assert_utils import assert_near_equal +from openmdao.utils.testing_utils import use_tempdirs +from openmdao.utils.units import convert_units + + +class _TestEOM(om.Group): + + def initialize(self): + self.options.declare('num_nodes', types=(int,)) + self.options.declare('foo_units', allow_none=True, default=None) + self.options.declare('foo_shapes', allow_none=True, default=None) + self.options.declare('foo_static', default=[]) + + def setup(self): + num_nodes = self.options['num_nodes'] + + foo_shape = (num_nodes,) if self.options['foo_shapes'] is None \ + else (num_nodes,) + self.options['foo_shapes']['vdot_comp'] + + foo_unit = 'kg' if self.options['foo_units'] is None else self.options['foo_units']['vdot_comp'] + + foo_tags = ['dymos.static_target']\ + if 'vdot_comp' in self.options['foo_static'] and self.options['foo_static']['vdot_comp'] else [] + + foo_shape = (1,) if 'dymos.static_target' in foo_tags else foo_shape + + vdot_comp = om.ExecComp(['vdot = g * cos(theta)', + 'bar = foo'], + vdot={'shape': (num_nodes,), 'units': 'm/s**2'}, + g={'val': 9.80665, 'units': 'm/s**2'}, + theta={'shape': (num_nodes,), 'units': 'rad'}, + foo={'shape': foo_shape, 'units': foo_unit, 'tags': foo_tags}, + bar={'shape': foo_shape, 'units': foo_unit}) + + foo_shape = (num_nodes,) if self.options['foo_shapes'] is None \ + else (num_nodes,) + self.options['foo_shapes']['xdot_comp'] + + foo_unit = 'kg' if self.options['foo_units'] is None else self.options['foo_units']['xdot_comp'] + + foo_tags = ['dymos.static_target']\ + if 'xdot_comp' in self.options['foo_static'] and self.options['foo_static']['xdot_comp'] else [] + + foo_shape = (1,) if 'dymos.static_target' in foo_tags else foo_shape + + xdot_comp = om.ExecComp(['xdot = v * sin(theta)', + 'bar = foo'], + xdot={'shape': (num_nodes,), 'units': 'm/s'}, + v={'shape': (num_nodes,), 'units': 'm/s'}, + theta={'shape': (num_nodes,), 'units': 'rad'}, + foo={'shape': foo_shape, 'units': foo_unit, 'tags': foo_tags}, + bar={'shape': foo_shape, 'units': foo_unit}) + + foo_shape = (num_nodes,) if self.options['foo_shapes'] is None \ + else (num_nodes,) + self.options['foo_shapes']['ydot_comp'] + + foo_unit = 'kg' if self.options['foo_units'] is None else self.options['foo_units']['ydot_comp'] + + foo_tags = ['dymos.static_target']\ + if 'ydot_comp' in self.options['foo_static'] and self.options['foo_static']['ydot_comp'] else [] + + foo_shape = (1,) if 'dymos.static_target' in foo_tags else foo_shape + + ydot_comp = om.ExecComp(['ydot = -v * cos(theta)', + 'bar = foo'], + ydot={'shape': (num_nodes,), 'units': 'm/s'}, + v={'shape': (num_nodes,), 'units': 'm/s'}, + theta={'shape': (num_nodes,), 'units': 'rad'}, + foo={'shape': foo_shape, 'units': foo_unit, 'tags': foo_tags}, + bar={'shape': foo_shape, 'units': foo_unit}) + + self.add_subsystem('vdot_comp', vdot_comp) + self.add_subsystem('xdot_comp', xdot_comp) + self.add_subsystem('ydot_comp', ydot_comp) + + +@use_tempdirs +class TestParameterShapes(unittest.TestCase): + + def test_valid_parameters(self): + import numpy as np + import openmdao.api as om + import dymos as dm + + # + # Define the OpenMDAO problem + # + p = om.Problem(model=om.Group()) + + # + # Define a Trajectory object + # + traj = dm.Trajectory() + + p.model.add_subsystem('traj', subsys=traj) + + # + # Define a Dymos Phase object with GaussLobatto Transcription + # + phase = dm.Phase(ode_class=_TestEOM, + transcription=dm.GaussLobatto(num_segments=10, order=3)) + + traj.add_phase(name='phase0', phase=phase) + + # + # Set the time options + # Time has no targets in our ODE. + # We fix the initial time so that it is not a design variable in the optimization. + # The duration of the phase is allowed to be optimized, but is bounded on [0.5, 10]. + # + phase.set_time_options(fix_initial=True, duration_bounds=(0.5, 10.0), units='s') + + # + # Set the time options + # Initial values of positions and velocity are all fixed. + # The final value of position are fixed, but the final velocity is a free variable. + # The equations of motion are not functions of position, so 'x' and 'y' have no targets. + # The rate source points to the output in the ODE which provides the time derivative of the + # given state. + phase.add_state('x', fix_initial=True, fix_final=True, units='m', + rate_source='xdot_comp.xdot') + phase.add_state('y', fix_initial=True, fix_final=True, units='m', + rate_source='ydot_comp.ydot') + phase.add_state('v', fix_initial=True, fix_final=False, units='m/s', + rate_source='vdot_comp.vdot', targets=['xdot_comp.v', 'ydot_comp.v']) + + # Define theta as a control. + phase.add_control(name='theta', units='rad', lower=0, upper=np.pi, + targets=['xdot_comp.theta', 'ydot_comp.theta', 'vdot_comp.theta']) + + phase.add_parameter('foo', + opt=False, + targets=['xdot_comp.foo', 'ydot_comp.foo', 'vdot_comp.foo']) + + # Minimize final time. + phase.add_objective('time', loc='final') + + # Set the driver. + p.driver = om.ScipyOptimizeDriver() + + # Allow OpenMDAO to automatically determine our sparsity pattern. + # Doing so can significant speed up the execution of Dymos. + p.driver.declare_coloring() + + # Setup the problem + p.setup(check=True) + + # Now that the OpenMDAO problem is setup, we can set the values of the states. + p.set_val('traj.phase0.t_initial', 0.0, units='s') + p.set_val('traj.phase0.t_duration', 2.0, units='s') + p.set_val('traj.phase0.states:x', phase.interp('x', [0, 10]), units='m') + p.set_val('traj.phase0.states:y', phase.interp('y', [10, 5]), units='m') + p.set_val('traj.phase0.states:v', phase.interp('v', [0, 5]), units='m/s') + p.set_val('traj.phase0.controls:theta', phase.interp('theta', [90, 90]), units='deg') + p.set_val('traj.phase0.parameters:foo', 5.0) + + # Run the driver to solve the problem + p.run_driver() + + self.assertEqual((1,), phase.parameter_options['foo']['shape']) + assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.parameter_vals:foo')[-1], 5.0, tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_col.xdot_comp.foo'), 5.0*np.ones(10,), tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_col.ydot_comp.foo'), 5.0*np.ones(10,), tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_col.vdot_comp.foo'), 5.0*np.ones(10,), tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_disc.xdot_comp.foo'), 5.0*np.ones(20,), tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_disc.ydot_comp.foo'), 5.0*np.ones(20,), tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_disc.vdot_comp.foo'), 5.0*np.ones(20,), tolerance=1.0E-5) + + def test_invalid_params_different_target_shapes(self): + # + # Define the OpenMDAO problem + # + p = om.Problem(model=om.Group()) + + # + # Define a Trajectory object + # + traj = dm.Trajectory() + + p.model.add_subsystem('traj', subsys=traj) + + # + # Define a Dymos Phase object with GaussLobatto Transcription + # + phase = dm.Phase(ode_class=_TestEOM, + transcription=dm.GaussLobatto(num_segments=10, order=3), + ode_init_kwargs={'foo_shapes': {'xdot_comp': (2,), + 'ydot_comp': (2,), + 'vdot_comp': (2,)}}) + + traj.add_phase(name='phase0', phase=phase) + + # + # Set the time options + # Time has no targets in our ODE. + # We fix the initial time so that it is not a design variable in the optimization. + # The duration of the phase is allowed to be optimized, but is bounded on [0.5, 10]. + # + phase.set_time_options(fix_initial=True, duration_bounds=(0.5, 10.0), units='s') + + # + # Set the time options + # Initial values of positions and velocity are all fixed. + # The final value of position are fixed, but the final velocity is a free variable. + # The equations of motion are not functions of position, so 'x' and 'y' have no targets. + # The rate source points to the output in the ODE which provides the time derivative of the + # given state. + phase.add_state('x', fix_initial=True, fix_final=True, units='m', + rate_source='xdot_comp.xdot') + phase.add_state('y', fix_initial=True, fix_final=True, units='m', + rate_source='ydot_comp.ydot') + phase.add_state('v', fix_initial=True, fix_final=False, units='m/s', + rate_source='vdot_comp.vdot', targets=['xdot_comp.v', 'ydot_comp.v']) + + # Define theta as a control. + phase.add_control(name='theta', units='rad', lower=0, upper=np.pi, + targets=['xdot_comp.theta', 'ydot_comp.theta', 'vdot_comp.theta']) + + phase.add_parameter('foo', shape=(1,), + opt=False, + targets=['xdot_comp.foo', 'ydot_comp.foo', 'vdot_comp.foo']) + + # Minimize final time. + phase.add_objective('time', loc='final') + + # Set the driver. + p.driver = om.ScipyOptimizeDriver() + + # Allow OpenMDAO to automatically determine our sparsity pattern. + # Doing so can significant speed up the execution of Dymos. + p.driver.declare_coloring() + + # Setup the problem + with self.assertRaises(RuntimeError) as e: + p.setup(check=True) + + expected = ("Shape provided to parameter `foo` differs from its targets.\n" + "Given shape: (1,)\n" + "Target shapes:\n" + "{'xdot_comp.foo': (2,), 'ydot_comp.foo': (2,), 'vdot_comp.foo': (2,)}") + self.assertEqual(expected, str(e.exception)) + + def test_invalid_params_different_target_shapes_introspection_failure(self): + # + # Define the OpenMDAO problem + # + p = om.Problem(model=om.Group()) + + # + # Define a Trajectory object + # + traj = dm.Trajectory() + + p.model.add_subsystem('traj', subsys=traj) + + # + # Define a Dymos Phase object with GaussLobatto Transcription + # + phase = dm.Phase(ode_class=_TestEOM, + transcription=dm.GaussLobatto(num_segments=10, order=3), + ode_init_kwargs={'foo_shapes': {'xdot_comp': (1,), + 'ydot_comp': (2,), + 'vdot_comp': (3,)}}) + + traj.add_phase(name='phase0', phase=phase) + + # + # Set the time options + # Time has no targets in our ODE. + # We fix the initial time so that it is not a design variable in the optimization. + # The duration of the phase is allowed to be optimized, but is bounded on [0.5, 10]. + # + phase.set_time_options(fix_initial=True, duration_bounds=(0.5, 10.0), units='s') + + # + # Set the time options + # Initial values of positions and velocity are all fixed. + # The final value of position are fixed, but the final velocity is a free variable. + # The equations of motion are not functions of position, so 'x' and 'y' have no targets. + # The rate source points to the output in the ODE which provides the time derivative of the + # given state. + phase.add_state('x', fix_initial=True, fix_final=True, units='m', + rate_source='xdot_comp.xdot') + phase.add_state('y', fix_initial=True, fix_final=True, units='m', + rate_source='ydot_comp.ydot') + phase.add_state('v', fix_initial=True, fix_final=False, units='m/s', + rate_source='vdot_comp.vdot', targets=['xdot_comp.v', 'ydot_comp.v']) + + # Define theta as a control. + phase.add_control(name='theta', units='rad', lower=0, upper=np.pi, + targets=['xdot_comp.theta', 'ydot_comp.theta', 'vdot_comp.theta']) + + phase.add_parameter('foo', + opt=False, + targets=['xdot_comp.foo', 'ydot_comp.foo', 'vdot_comp.foo']) + + # Minimize final time. + phase.add_objective('time', loc='final') + + # Set the driver. + p.driver = om.ScipyOptimizeDriver() + + # Allow OpenMDAO to automatically determine our sparsity pattern. + # Doing so can significant speed up the execution of Dymos. + p.driver.declare_coloring() + + # Setup the problem + with self.assertRaises(RuntimeError) as e: + p.setup(check=True) + + expected = ('Invalid targets for parameter `foo`.\n' + 'Targets have multiple shapes.\n' + "{'xdot_comp.foo': (1,), 'ydot_comp.foo': (2,), 'vdot_comp.foo': (3,)}") + self.assertEqual(expected, str(e.exception)) + + +@use_tempdirs +class TestParameterUnits(unittest.TestCase): + + def test_valid_parameters(self): + import numpy as np + import openmdao.api as om + import dymos as dm + + # + # Define the OpenMDAO problem + # + p = om.Problem(model=om.Group()) + + # + # Define a Trajectory object + # + traj = dm.Trajectory() + + p.model.add_subsystem('traj', subsys=traj) + + # + # Define a Dymos Phase object with GaussLobatto Transcription + # + phase = dm.Phase(ode_class=_TestEOM, + transcription=dm.GaussLobatto(num_segments=10, order=3)) + + traj.add_phase(name='phase0', phase=phase) + + # + # Set the time options + # Time has no targets in our ODE. + # We fix the initial time so that it is not a design variable in the optimization. + # The duration of the phase is allowed to be optimized, but is bounded on [0.5, 10]. + # + phase.set_time_options(fix_initial=True, duration_bounds=(0.5, 10.0), units='s') + + # + # Set the time options + # Initial values of positions and velocity are all fixed. + # The final value of position are fixed, but the final velocity is a free variable. + # The equations of motion are not functions of position, so 'x' and 'y' have no targets. + # The rate source points to the output in the ODE which provides the time derivative of the + # given state. + phase.add_state('x', fix_initial=True, fix_final=True, units='m', + rate_source='xdot_comp.xdot') + phase.add_state('y', fix_initial=True, fix_final=True, units='m', + rate_source='ydot_comp.ydot') + phase.add_state('v', fix_initial=True, fix_final=False, units='m/s', + rate_source='vdot_comp.vdot', targets=['xdot_comp.v', 'ydot_comp.v']) + + # Define theta as a control. + phase.add_control(name='theta', units='rad', lower=0, upper=np.pi, + targets=['xdot_comp.theta', 'ydot_comp.theta', 'vdot_comp.theta']) + + phase.add_parameter('foo', + units='lbm', + opt=False, + targets=['xdot_comp.foo', 'ydot_comp.foo', 'vdot_comp.foo']) + + # Minimize final time. + phase.add_objective('time', loc='final') + + # Set the driver. + p.driver = om.ScipyOptimizeDriver() + + # Allow OpenMDAO to automatically determine our sparsity pattern. + # Doing so can significant speed up the execution of Dymos. + p.driver.declare_coloring() + + # Setup the problem + p.setup(check=True) + + # Now that the OpenMDAO problem is setup, we can set the values of the states. + p.set_val('traj.phase0.t_initial', 0.0, units='s') + p.set_val('traj.phase0.t_duration', 2.0, units='s') + p.set_val('traj.phase0.states:x', phase.interp('x', [0, 10]), units='m') + p.set_val('traj.phase0.states:y', phase.interp('y', [10, 5]), units='m') + p.set_val('traj.phase0.states:v', phase.interp('v', [0, 5]), units='m/s') + p.set_val('traj.phase0.controls:theta', phase.interp('theta', [90, 90]), units='deg') + p.set_val('traj.phase0.parameters:foo', 5.0) + + # Run the driver to solve the problem + p.run_driver() + + expected = convert_units(5.0, 'lbm', 'kg') + + self.assertEqual((1,), phase.parameter_options['foo']['shape']) + assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.parameter_vals:foo')[-1], 5.0, tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_col.xdot_comp.foo'), expected*np.ones(10,), tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_col.ydot_comp.foo'), expected*np.ones(10,), tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_col.vdot_comp.foo'), expected*np.ones(10,), tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_disc.xdot_comp.foo'), expected*np.ones(20,), tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_disc.ydot_comp.foo'), expected*np.ones(20,), tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_disc.vdot_comp.foo'), expected*np.ones(20,), tolerance=1.0E-5) + + def test_invalid_params_different_target_units_introspection_failure(self): + # + # Define the OpenMDAO problem + # + p = om.Problem(model=om.Group()) + + # + # Define a Trajectory object + # + traj = dm.Trajectory() + + p.model.add_subsystem('traj', subsys=traj) + + # + # Define a Dymos Phase object with GaussLobatto Transcription + # + phase = dm.Phase(ode_class=_TestEOM, + transcription=dm.GaussLobatto(num_segments=10, order=3), + ode_init_kwargs={'foo_units': {'xdot_comp': 'kg', + 'ydot_comp': 'lbm', + 'vdot_comp': 'slug'}}) + + traj.add_phase(name='phase0', phase=phase) + + # + # Set the time options + # Time has no targets in our ODE. + # We fix the initial time so that it is not a design variable in the optimization. + # The duration of the phase is allowed to be optimized, but is bounded on [0.5, 10]. + # + phase.set_time_options(fix_initial=True, duration_bounds=(0.5, 10.0), units='s') + + # + # Set the time options + # Initial values of positions and velocity are all fixed. + # The final value of position are fixed, but the final velocity is a free variable. + # The equations of motion are not functions of position, so 'x' and 'y' have no targets. + # The rate source points to the output in the ODE which provides the time derivative of the + # given state. + phase.add_state('x', fix_initial=True, fix_final=True, units='m', + rate_source='xdot_comp.xdot') + phase.add_state('y', fix_initial=True, fix_final=True, units='m', + rate_source='ydot_comp.ydot') + phase.add_state('v', fix_initial=True, fix_final=False, units='m/s', + rate_source='vdot_comp.vdot', targets=['xdot_comp.v', 'ydot_comp.v']) + + # Define theta as a control. + phase.add_control(name='theta', units='rad', lower=0, upper=np.pi, + targets=['xdot_comp.theta', 'ydot_comp.theta', 'vdot_comp.theta']) + + phase.add_parameter('foo', + opt=False, + targets=['xdot_comp.foo', 'ydot_comp.foo', 'vdot_comp.foo']) + + # Minimize final time. + phase.add_objective('time', loc='final') + + # Set the driver. + p.driver = om.ScipyOptimizeDriver() + + # Allow OpenMDAO to automatically determine our sparsity pattern. + # Doing so can significant speed up the execution of Dymos. + p.driver.declare_coloring() + + # Setup the problem + with self.assertRaises(RuntimeError) as e: + p.setup(check=True) + + expected = ("Unable to automatically assign units based on targets.\n" + "Targets have multiple units assigned:\n" + "{'xdot_comp.foo': 'kg', 'ydot_comp.foo': 'lbm', 'vdot_comp.foo': 'slug'}.\n" + "Either promote targets and use set_input_defaults to assign common\n" + "units, or explicitly provide units to the variable.") + + self.assertEqual(expected, str(e.exception)) + + +@use_tempdirs +class TestMixedStaticDynamicParameterTargets(unittest.TestCase): + + def test_all_static(self): + import numpy as np + import openmdao.api as om + import dymos as dm + + # + # Define the OpenMDAO problem + # + p = om.Problem(model=om.Group()) + + # + # Define a Trajectory object + # + traj = dm.Trajectory() + + p.model.add_subsystem('traj', subsys=traj) + + # + # Define a Dymos Phase object with GaussLobatto Transcription + # + phase = dm.Phase(ode_class=_TestEOM, + transcription=dm.GaussLobatto(num_segments=10, order=3), + ode_init_kwargs={'foo_static': {'xdot_comp': True, 'ydot_comp': True, 'vdot_comp': True}}) + + traj.add_phase(name='phase0', phase=phase) + + # + # Set the time options + # Time has no targets in our ODE. + # We fix the initial time so that it is not a design variable in the optimization. + # The duration of the phase is allowed to be optimized, but is bounded on [0.5, 10]. + # + phase.set_time_options(fix_initial=True, duration_bounds=(0.5, 10.0), units='s') + + # + # Set the time options + # Initial values of positions and velocity are all fixed. + # The final value of position are fixed, but the final velocity is a free variable. + # The equations of motion are not functions of position, so 'x' and 'y' have no targets. + # The rate source points to the output in the ODE which provides the time derivative of the + # given state. + phase.add_state('x', fix_initial=True, fix_final=True, units='m', + rate_source='xdot_comp.xdot') + phase.add_state('y', fix_initial=True, fix_final=True, units='m', + rate_source='ydot_comp.ydot') + phase.add_state('v', fix_initial=True, fix_final=False, units='m/s', + rate_source='vdot_comp.vdot', targets=['xdot_comp.v', 'ydot_comp.v']) + + # Define theta as a control. + phase.add_control(name='theta', units='rad', lower=0, upper=np.pi, + targets=['xdot_comp.theta', 'ydot_comp.theta', 'vdot_comp.theta']) + + phase.add_parameter('foo', + units='lbm', + opt=False, + targets=['xdot_comp.foo', 'ydot_comp.foo', 'vdot_comp.foo']) + + # Minimize final time. + phase.add_objective('time', loc='final') + + # Set the driver. + p.driver = om.ScipyOptimizeDriver() + + # Allow OpenMDAO to automatically determine our sparsity pattern. + # Doing so can significant speed up the execution of Dymos. + p.driver.declare_coloring() + + # Setup the problem + p.setup(check=True) + + # Now that the OpenMDAO problem is setup, we can set the values of the states. + p.set_val('traj.phase0.t_initial', 0.0, units='s') + p.set_val('traj.phase0.t_duration', 2.0, units='s') + p.set_val('traj.phase0.states:x', phase.interp('x', [0, 10]), units='m') + p.set_val('traj.phase0.states:y', phase.interp('y', [10, 5]), units='m') + p.set_val('traj.phase0.states:v', phase.interp('v', [0, 5]), units='m/s') + p.set_val('traj.phase0.controls:theta', phase.interp('theta', [90, 90]), units='deg') + p.set_val('traj.phase0.parameters:foo', 5.0) + + # Run the driver to solve the problem + p.run_driver() + + expected = convert_units(5.0, 'lbm', 'kg') + + self.assertEqual((1,), phase.parameter_options['foo']['shape']) + assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.parameter_vals:foo')[-1], 5.0, tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_col.xdot_comp.foo'), expected, tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_col.ydot_comp.foo'), expected, tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_col.vdot_comp.foo'), expected, tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_disc.xdot_comp.foo'), expected, tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_disc.ydot_comp.foo'), expected, tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_disc.vdot_comp.foo'), expected, tolerance=1.0E-5) + + def test_mixed_static(self): + import numpy as np + import openmdao.api as om + import dymos as dm + + # + # Define the OpenMDAO problem + # + p = om.Problem(model=om.Group()) + + # + # Define a Trajectory object + # + traj = dm.Trajectory() + + p.model.add_subsystem('traj', subsys=traj) + + # + # Define a Dymos Phase object with GaussLobatto Transcription + # + phase = dm.Phase(ode_class=_TestEOM, + transcription=dm.GaussLobatto(num_segments=10, order=3), + ode_init_kwargs={'foo_static': {'xdot_comp': True, 'ydot_comp': True, 'vdot_comp': False}}) + + traj.add_phase(name='phase0', phase=phase) + + # + # Set the time options + # Time has no targets in our ODE. + # We fix the initial time so that it is not a design variable in the optimization. + # The duration of the phase is allowed to be optimized, but is bounded on [0.5, 10]. + # + phase.set_time_options(fix_initial=True, duration_bounds=(0.5, 10.0), units='s') + + # + # Set the time options + # Initial values of positions and velocity are all fixed. + # The final value of position are fixed, but the final velocity is a free variable. + # The equations of motion are not functions of position, so 'x' and 'y' have no targets. + # The rate source points to the output in the ODE which provides the time derivative of the + # given state. + phase.add_state('x', fix_initial=True, fix_final=True, units='m', + rate_source='xdot_comp.xdot') + phase.add_state('y', fix_initial=True, fix_final=True, units='m', + rate_source='ydot_comp.ydot') + phase.add_state('v', fix_initial=True, fix_final=False, units='m/s', + rate_source='vdot_comp.vdot', targets=['xdot_comp.v', 'ydot_comp.v']) + + # Define theta as a control. + phase.add_control(name='theta', units='rad', lower=0, upper=np.pi, + targets=['xdot_comp.theta', 'ydot_comp.theta', 'vdot_comp.theta']) + + phase.add_parameter('foo', + units='lbm', + opt=False, + targets=['xdot_comp.foo', 'ydot_comp.foo', 'vdot_comp.foo']) + + # Minimize final time. + phase.add_objective('time', loc='final') + + # Set the driver. + p.driver = om.ScipyOptimizeDriver() + + # Allow OpenMDAO to automatically determine our sparsity pattern. + # Doing so can significant speed up the execution of Dymos. + p.driver.declare_coloring() + + # Setup the problem + p.setup(check=True) + + # Now that the OpenMDAO problem is setup, we can set the values of the states. + p.set_val('traj.phase0.t_initial', 0.0, units='s') + p.set_val('traj.phase0.t_duration', 2.0, units='s') + p.set_val('traj.phase0.states:x', phase.interp('x', [0, 10]), units='m') + p.set_val('traj.phase0.states:y', phase.interp('y', [10, 5]), units='m') + p.set_val('traj.phase0.states:v', phase.interp('v', [0, 5]), units='m/s') + p.set_val('traj.phase0.controls:theta', phase.interp('theta', [90, 90]), units='deg') + p.set_val('traj.phase0.parameters:foo', 5.0) + + # Run the driver to solve the problem + p.run_driver() + + expected = convert_units(5.0, 'lbm', 'kg') + + self.assertEqual((1,), phase.parameter_options['foo']['shape']) + self.assertEqual({'xdot_comp.foo', 'ydot_comp.foo'}, phase.parameter_options['foo']['static_targets']) + assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.parameter_vals:foo')[-1], 5.0, tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_col.xdot_comp.foo'), expected, tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_col.ydot_comp.foo'), expected, tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_col.vdot_comp.foo'), expected*np.ones(10), tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_disc.xdot_comp.foo'), expected, tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_disc.ydot_comp.foo'), expected, tolerance=1.0E-5) + assert_near_equal(p.get_val('traj.phase0.rhs_disc.vdot_comp.foo'), expected*np.ones(20), tolerance=1.0E-5) + + +if __name__ == '__main__': # pragma: no cover + unittest.main() diff --git a/dymos/examples/cannonball/test/test_connect_control_to_parameter.py b/dymos/examples/cannonball/test/test_connect_control_to_parameter.py index 340c31855..f447bb770 100644 --- a/dymos/examples/cannonball/test/test_connect_control_to_parameter.py +++ b/dymos/examples/cannonball/test/test_connect_control_to_parameter.py @@ -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 diff --git a/dymos/examples/cannonball/test/test_traj_param_static_and_dynamic.py b/dymos/examples/cannonball/test/test_traj_param_static_and_dynamic.py new file mode 100644 index 000000000..d717576b8 --- /dev/null +++ b/dymos/examples/cannonball/test/test_traj_param_static_and_dynamic.py @@ -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() diff --git a/dymos/grid_refinement/grid_refinement_ode_system.py b/dymos/grid_refinement/grid_refinement_ode_system.py index 3fb76b900..fd86ec729 100644 --- a/dymos/grid_refinement/grid_refinement_ode_system.py +++ b/dymos/grid_refinement/grid_refinement_ode_system.py @@ -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 @@ -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, diff --git a/dymos/grid_refinement/hp_adaptive/hp_adaptive.py b/dymos/grid_refinement/hp_adaptive/hp_adaptive.py index a8a0e2bf2..b675fe96c 100644 --- a/dymos/grid_refinement/hp_adaptive/hp_adaptive.py +++ b/dymos/grid_refinement/hp_adaptive/hp_adaptive.py @@ -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'] diff --git a/dymos/phase/options.py b/dymos/phase/options.py index 48fe7cead..6f343806c 100644 --- a/dymos/phase/options.py +++ b/dymos/phase/options.py @@ -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, diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py index bb32cb829..744e54a67 100644 --- a/dymos/phase/phase.py +++ b/dymos/phase/phase.py @@ -57,7 +57,6 @@ class Phase(om.Group): **kwargs : dict Dictionary of optional phase arguments. """ - def __init__(self, from_phase=None, **kwargs): _kwargs = kwargs.copy() @@ -902,7 +901,7 @@ def add_parameter(self, name, val=_unspecified, units=_unspecified, opt=False, desc=_unspecified, lower=_unspecified, upper=_unspecified, scaler=_unspecified, adder=_unspecified, ref0=_unspecified, ref=_unspecified, targets=_unspecified, shape=_unspecified, dynamic=_unspecified, static_target=_unspecified, - include_timeseries=_unspecified): + include_timeseries=_unspecified, static_targets=_unspecified): """ Add a parameter (static control variable) to the phase. @@ -948,6 +947,12 @@ def add_parameter(self, name, val=_unspecified, units=_unspecified, opt=False, (meaning they cannot have a unique value at each node). Otherwise False. include_timeseries : bool True if the static parameters should be included in output timeseries, else False. + static_targets : bool or Sequence or _unspecified + True if ALL targets in the ODE are not shaped with num_nodes as the first dimension + (meaning they cannot have a unique value at each node). If False, ALL targets are + expected to be shaped with the first dimension as the number of nodes. If given + as a Sequence, it provides those targets not shaped with num_nodes. If left _unspecified, + static targets will be determined automatically. """ self.check_parameter(name) @@ -955,16 +960,18 @@ def add_parameter(self, name, val=_unspecified, units=_unspecified, opt=False, self.parameter_options[name] = ParameterOptionsDictionary() self.parameter_options[name]['name'] = name - self.set_parameter_options(name, val, units, opt, desc, lower, upper, - scaler, adder, ref0, ref, targets, shape, dynamic, - static_target, include_timeseries) + self.set_parameter_options(name, val=val, units=units, opt=opt, desc=desc, + lower=lower, upper=upper, scaler=scaler, adder=adder, + ref0=ref0, ref=ref, targets=targets, shape=shape, dynamic=dynamic, + static_target=static_target, static_targets=static_targets, + include_timeseries=include_timeseries) def set_parameter_options(self, name, val=_unspecified, units=_unspecified, opt=False, desc=_unspecified, lower=_unspecified, upper=_unspecified, scaler=_unspecified, adder=_unspecified, ref0=_unspecified, ref=_unspecified, targets=_unspecified, shape=_unspecified, dynamic=_unspecified, static_target=_unspecified, - include_timeseries=_unspecified): + include_timeseries=_unspecified, static_targets=_unspecified): """ Set options for an existing parameter (static control variable) in the phase. @@ -1010,6 +1017,12 @@ def set_parameter_options(self, name, val=_unspecified, units=_unspecified, opt= (meaning they cannot have a unique value at each node). Otherwise False. include_timeseries : bool True if the static parameters should be included in output timeseries, else False. + static_targets : bool or Sequence or _unspecified + True if ALL targets in the ODE are not shaped with num_nodes as the first dimension + (meaning they cannot have a unique value at each node). If False, ALL targets are + expected to be shaped with the first dimension as the number of nodes. If given + as a Sequence, it provides those targets not shaped with num_nodes. If left _unspecified, + static targets will be determined automatically. """ if units is not _unspecified: self.parameter_options[name]['units'] = units @@ -1035,22 +1048,29 @@ def set_parameter_options(self, name, val=_unspecified, units=_unspecified, opt= self.parameter_options[name]['shape'] = tuple(shape) else: self.parameter_options[name]['shape'] = shape - elif val is not _unspecified: - if isinstance(val, float) or isinstance(val, int) or isinstance(val, complex): - self.parameter_options[name]['shape'] = (1,) - else: - self.parameter_options[name]['shape'] = tuple(np.asarray(val).shape) if dynamic is not _unspecified: self.parameter_options[name]['static_target'] = not dynamic + self.parameter_options[name]['static_targets'] = not dynamic if static_target is not _unspecified: self.parameter_options[name]['static_target'] = static_target + self.parameter_options[name]['static_targets'] = static_target + + if static_targets is not _unspecified: + self.parameter_options[name]['static_target'] = static_targets + self.parameter_options[name]['static_targets'] = static_targets if dynamic is not _unspecified and static_target is not _unspecified: - raise ValueError("Both the deprecated 'dynamic' option and option 'static_target' were " + raise ValueError("Both the deprecated 'dynamic' option and option 'static_target' were\n" + f"specified for parameter '{name}'. Going forward, please use only\n" + "option static_targets. Options 'dynamic' and 'static_target'\n" + "will be removed in Dymos 2.0.0.") + + if dynamic is not _unspecified and static_targets is not _unspecified: + raise ValueError("Both the deprecated 'dynamic' option and option 'static_targets' were " f"specified for parameter '{name}'. Going forward, please use only " - "option static_target. Option 'dynamic' will be removed in " + "option static_targets. Option 'dynamic' will be removed in " "Dymos 2.0.0.") if lower is not _unspecified: @@ -2292,10 +2312,7 @@ def initialize_values_from_phase(self, prob, from_phase, phase_path='', skip_par # We use this private function to grab the correctly sized variable from the # auto_ivc source. - if om_version < (3, 4, 1): - val = phs.get_val(f'parameters:{name}', units=units)[0, ...] - else: - val = phs.get_val(f'parameters:{name}', units=units) + val = phs.get_val(f'parameters:{name}', units=units) if phase_path: prob_path = f'{phase_path}.{self.name}.parameters:{name}' diff --git a/dymos/phase/test/test_input_parameter_connections.py b/dymos/phase/test/test_input_parameter_connections.py index 99dbbc168..34cd53b63 100644 --- a/dymos/phase/test/test_input_parameter_connections.py +++ b/dymos/phase/test/test_input_parameter_connections.py @@ -105,11 +105,12 @@ def test_static_and_dynamic_error(self): phase.add_parameter('alpha', val=np.ones((n_traj, 2)), units='m', targets='comp.alpha', dynamic=False, static_target=True) - expected_msg = "Both the deprecated 'dynamic' option and option 'static_target' were " \ - "specified for parameter 'alpha'. Going forward, please use only option " \ - "static_target. Option 'dynamic' will be removed in Dymos 2.0.0." + expected_msg = "Both the deprecated 'dynamic' option and option 'static_target' were\n" \ + "specified for parameter 'alpha'. Going forward, please use only\n" \ + "option static_targets. Options 'dynamic' and 'static_target'\n" \ + "will be removed in Dymos 2.0.0." - self.assertEqual(str(e.exception), expected_msg) + self.assertEqual(expected_msg, str(e.exception)) def test_static_and_dynamic_error_in_traj(self): @@ -119,11 +120,12 @@ def test_static_and_dynamic_error_in_traj(self): t.add_parameter('alpha', val=np.ones((n_traj, 2)), units='m', targets={'p': ['comp.alpha']}, dynamic=False, static_target=True) - expected_msg = "Both the deprecated 'dynamic' option and option 'static_target' were " \ - "specified for parameter 'alpha'. Going forward, please use only option " \ - "static_target. Option 'dynamic' will be removed in Dymos 2.0.0." + expected_msg = "Both the deprecated 'dynamic' option and option 'static_target' were\n" \ + "specified for parameter 'alpha'. Going forward, please use only\n" \ + "option static_targets. Options 'dynamic' and 'static_target'\n" \ + "will be removed in Dymos 2.0.0." - self.assertEqual(str(e.exception), expected_msg) + self.assertEqual(expected_msg, str(e.exception)) if __name__ == '__main__': # pragma: no cover diff --git a/dymos/phase/test/test_simulate.py b/dymos/phase/test/test_simulate.py index 4efe25c49..d25bf9ec2 100644 --- a/dymos/phase/test/test_simulate.py +++ b/dymos/phase/test/test_simulate.py @@ -4,6 +4,7 @@ import openmdao.api as om from openmdao.utils.testing_utils import use_tempdirs +from openmdao.utils.assert_utils import assert_near_equal import dymos as dm @@ -26,8 +27,10 @@ def initialize(self): def setup(self): nn = self.options['num_nodes'] self.add_subsystem('Sink', - om.ExecComp('sink = chord[0]', sink={'val': 0.0, 'units': None}, - chord={'val': np.zeros(4), 'units': 'm'}), + om.ExecComp('sink = chord[0] + chord[1] + chord[2] + chord[3]', + sink={'val': 0.0, 'units': None}, + chord={'val': np.zeros(4), 'units': 'm', + 'tags': ['dymos.static_target']}), promotes_inputs=['chord']) self.add_subsystem('calc', om.ExecComp('Out = Thrust * 2', @@ -58,8 +61,7 @@ def test_shaped_params(self): main_phase.set_time_options(fix_initial=True, fix_duration=True, units='s') - main_phase.add_parameter('chord', targets='chord', shape=(4,), units='inch', - static_target=True) + main_phase.add_parameter('chord', units='inch') p.model.connect('chord', 'hop0.main_phase.parameters:chord') main_phase.add_state('impulse', fix_initial=True, fix_final=False, units='N*s', @@ -75,13 +77,59 @@ def test_shaped_params(self): p.setup(mode='auto', check=['unconnected_inputs'], force_alloc_complex=True) - p['hop0.main_phase.t_initial'] = 0.0 - p['hop0.main_phase.t_duration'] = 10 - p['hop0.main_phase.polynomial_controls:Thrust'][:, 0] = -3400 - p['hop0.main_phase.states:impulse'] = main_phase.interp('impulse', [0, 0]) + p.set_val('hop0.main_phase.t_initial', 0.0) + p.set_val('hop0.main_phase.t_duration', 10) + p.set_val('hop0.main_phase.polynomial_controls:Thrust', val=-3400, indices=om.slicer[:, 0]) + p.set_val('hop0.main_phase.states:impulse', main_phase.interp('impulse', [0, 0])) p.run_driver() + assert_near_equal(p.get_val('hop0.main_phase.timeseries.impulse')[-1, 0], -7836.66666, tolerance=1.0E-4) + + try: + hop0.simulate() + except: + self.fail('Simulate did not correctly complete.') + + def test_shaped_traj_params(self): + + main_tx = dm.Radau(num_segments=1, order=3, compressed=False) + + p = om.Problem(model=om.Group()) + + p.driver = om.ScipyOptimizeDriver() + + hop0 = dm.Trajectory() + p.model.add_subsystem('hop0', hop0) + main_phase = hop0.add_phase(name='main_phase', + phase=MainPhase(transcription=main_tx)) + + main_phase.set_time_options(fix_initial=True, fix_duration=True, units='s') + + hop0.add_parameter('chord', units='inch', targets={'main_phase': ['chord']}) + + main_phase.add_state('impulse', fix_initial=True, fix_final=False, units='N*s', + rate_source='Out', + solve_segments=False) + + main_phase.add_polynomial_control('Thrust', units='N', + targets='Thrust', + lower=-3450, upper=-500, + order=5, opt=True) + + main_phase.add_objective('impulse', loc='final', ref=-1) + + p.setup(mode='auto', check=['unconnected_inputs'], force_alloc_complex=True) + + p.set_val('hop0.main_phase.t_initial', 0.0) + p.set_val('hop0.main_phase.t_duration', 10) + p.set_val('hop0.main_phase.polynomial_controls:Thrust', val=-3400, indices=om.slicer[:, 0]) + p.set_val('hop0.main_phase.states:impulse', main_phase.interp('impulse', [0, 0])) + + p.run_driver() + + assert_near_equal(p.get_val('hop0.main_phase.timeseries.impulse')[-1, 0], -7836.66666, tolerance=1.0E-4) + try: hop0.simulate() except: diff --git a/dymos/test/test_load_case.py b/dymos/test/test_load_case.py index 458cc3138..0661627ae 100644 --- a/dymos/test/test_load_case.py +++ b/dymos/test/test_load_case.py @@ -56,7 +56,6 @@ def setup_problem(trans=dm.GaussLobatto(num_segments=10), polynomial_control=Fal return p -@unittest.skipIf(om_version <= (2, 9, 0), 'load_case requires an OpenMDAO version later than 2.9.0') @use_tempdirs class TestLoadCase(unittest.TestCase): diff --git a/dymos/test/test_run_problem.py b/dymos/test/test_run_problem.py index 95fd0a950..f7b1b9e33 100755 --- a/dymos/test/test_run_problem.py +++ b/dymos/test/test_run_problem.py @@ -336,7 +336,6 @@ def test_illegal_simulate_kwargs(self): 'Key "case_prefix" was found in simulate_kwargs but should instead by provided by ' 'the argument "case_prefix", not part of the simulate_kwargs dictionary.') - @unittest.skipIf(om_version < (3, 18, 0), 'test requires OpenMDAO >= 3.18.01') @require_pyoptsparse(optimizer='SLSQP') def test_run_brachistochrone_problem_refine_case_driver_case_prefix(self): p = om.Problem(model=om.Group()) diff --git a/dymos/trajectory/trajectory.py b/dymos/trajectory/trajectory.py index 62de5b655..46ea5469c 100644 --- a/dymos/trajectory/trajectory.py +++ b/dymos/trajectory/trajectory.py @@ -20,8 +20,8 @@ from ..phase.analytic_phase import AnalyticPhase from ..phase.options import TrajParameterOptionsDictionary from ..transcriptions.common import ParameterComp -from ..utils.misc import get_rate_units, _unspecified -from ..utils.introspection import get_promoted_vars, get_source_metadata +from ..utils.misc import get_rate_units, _unspecified, _none_or_unspecified +from ..utils.introspection import get_promoted_vars, get_source_metadata, _get_common_metadata from .._options import options as dymos_options @@ -94,7 +94,8 @@ def add_phase(self, name, phase, **kwargs): def set_parameter_options(self, name, units=_unspecified, val=_unspecified, desc=_unspecified, opt=False, targets=_unspecified, lower=_unspecified, upper=_unspecified, scaler=_unspecified, adder=_unspecified, ref0=_unspecified, ref=_unspecified, - shape=_unspecified, dynamic=_unspecified, static_target=_unspecified): + shape=_unspecified, dynamic=_unspecified, static_target=_unspecified, + static_targets=_unspecified): """ Set the options of an existing a parameter in the trajectory. @@ -138,6 +139,12 @@ def set_parameter_options(self, name, units=_unspecified, val=_unspecified, desc static_target : bool or _unspecified True if the targets in the ODE are not shaped with num_nodes as the first dimension (meaning they cannot have a unique value at each node). Otherwise False. + static_targets : bool or Sequence or _unspecified + True if ALL targets in the ODE are not shaped with num_nodes as the first dimension + (meaning they cannot have a unique value at each node). If False, ALL targets are + expected to be shaped with the first dimension as the number of nodes. If given + as a Sequence, it provides those targets not shaped with num_nodes. If left _unspecified, + static targets will be determined automatically. """ if name not in self.parameter_options: self.parameter_options[name] = TrajParameterOptionsDictionary() @@ -186,17 +193,28 @@ def set_parameter_options(self, name, units=_unspecified, val=_unspecified, desc if static_target is not _unspecified: self.parameter_options[name]['static_target'] = static_target + self.parameter_options[name]['static_targets'] = static_target + + if static_targets is not _unspecified: + self.parameter_options[name]['static_targets'] = static_target if dynamic is not _unspecified and static_target is not _unspecified: - raise ValueError("Both the deprecated 'dynamic' option and option 'static_target' were " + raise ValueError("Both the deprecated 'dynamic' option and option 'static_target' were\n" + f"specified for parameter '{name}'. Going forward, please use only\n" + "option static_targets. Options 'dynamic' and 'static_target'\n" + "will be removed in Dymos 2.0.0.") + + if dynamic is not _unspecified and static_targets is not _unspecified: + raise ValueError("Both the deprecated 'dynamic' option and option 'static_targets' were " f"specified for parameter '{name}'. " - f"Going forward, please use only option static_target. Option " + f"Going forward, please use only option static_targets. Option " f"'dynamic' will be removed in Dymos 2.0.0.") def add_parameter(self, name, units=_unspecified, val=_unspecified, desc=_unspecified, opt=False, targets=_unspecified, lower=_unspecified, upper=_unspecified, scaler=_unspecified, adder=_unspecified, ref0=_unspecified, ref=_unspecified, - shape=_unspecified, dynamic=_unspecified, static_target=_unspecified): + shape=_unspecified, dynamic=_unspecified, static_target=_unspecified, + static_targets=_unspecified): """ Add a parameter (static control) to the trajectory. @@ -240,6 +258,12 @@ def add_parameter(self, name, units=_unspecified, val=_unspecified, desc=_unspec static_target : bool or _unspecified True if the targets in the ODE are not shaped with num_nodes as the first dimension (meaning they cannot have a unique value at each node). Otherwise False. + static_targets : bool or Sequence or _unspecified + True if ALL targets in the ODE are not shaped with num_nodes as the first dimension + (meaning they cannot have a unique value at each node). If False, ALL targets are + expected to be shaped with the first dimension as the number of nodes. If given + as a Sequence, it provides those targets not shaped with num_nodes. If left _unspecified, + static targets will be determined automatically. """ if name not in self.parameter_options: self.parameter_options[name] = TrajParameterOptionsDictionary() @@ -252,6 +276,29 @@ def add_parameter(self, name, units=_unspecified, val=_unspecified, desc=_unspec upper=upper, scaler=scaler, adder=adder, ref0=ref0, ref=ref, shape=shape, dynamic=dynamic, static_target=static_target) + def _get_phase_parameters(self): + """ + Retrieve a dict of parameter options for each phase within the trajectory. + + Returns + ------- + dict + A dictionary keyed by phase name. Each associated value is a dictionary + keyed by parameter name and the associated values are parameter options + for each parameter. + + """ + phase_param_options = {} + for phs in self.phases._subsystems_myproc: + phase_param_options[phs.name] = phs.parameter_options + + if self.comm.size > 1: + data = self.comm.allgather(phase_param_options) + if data: + for d in data: + phase_param_options.update(d) + return phase_param_options + def _setup_parameters(self): """ Adds an IndepVarComp if necessary and issues appropriate connections based @@ -286,7 +333,7 @@ def _setup_parameters(self): # We need to add an input parameter to this phase. # The default target in the phase is name unless otherwise specified. - kwargs = {'static_target': options['static_target'], + kwargs = {'static_targets': options['static_targets'], 'units': options['units'], 'val': options['val'], 'shape': options['shape'], @@ -371,7 +418,6 @@ def _configure_parameters(self): """ parameter_options = self.parameter_options promoted_inputs = [] - for name, options in parameter_options.items(): promoted_inputs.append(f'parameters:{name}') targets = options['targets'] @@ -390,19 +436,16 @@ def _configure_parameters(self): # For each phase, use introspection to get the units and shape. # If units do not match across all phases, require user to set them. # If shapes do not match across all phases, this is an error. - tgts = [] - tgt_units = {} - tgt_shapes = {} + targets_per_phase = {} for phase_name, phs in self._phases.items(): + target_param = None if targets is None or phase_name not in targets: # Attempt to connect to an input parameter of the same name in the phase, if # it exists. if name in phs.parameter_options: - tgt = f'{phase_name}.parameters:{name}' - tgt_shapes[phs.name] = phs.parameter_options[name]['shape'] - tgt_units[phs.name] = phs.parameter_options[name]['units'] + target_param = name else: continue elif targets[phase_name] is None: @@ -411,9 +454,7 @@ def _configure_parameters(self): elif isinstance(targets[phase_name], str): if targets[phase_name] in phs.parameter_options: # Connect to an input parameter with a different name in this phase - tgt = f'{phase_name}.parameters:{targets[phase_name]}' - tgt_shapes[phs.name] = phs.parameter_options[targets[phase_name]]['shape'] - tgt_units[phs.name] = phs.parameter_options[targets[phase_name]]['units'] + target_param = targets[phase_name] else: msg = f'Invalid target for trajectory `{self.pathname}` parameter `{name}` in phase ' \ f"`{phase_name}`.\nTarget for phase `{phase_name}` is '{targets[phase_name]}' but " \ @@ -423,9 +464,7 @@ def _configure_parameters(self): if name in phs.parameter_options: # User gave a list of ODE targets which were passed to the creation of a # new input parameter in setup, just connect to that new input parameter - tgt = f'{phase_name}.parameters:{name}' - tgt_shapes[phs.name] = phs.parameter_options[name]['shape'] - tgt_units[phs.name] = phs.parameter_options[name]['units'] + target_param = name else: msg = f'Invalid target for trajectory `{self.pathname}` parameter `{name}` in phase ' \ f"`{phase_name}`.\nThe phase did not add the parameter as expected. Please file an " \ @@ -435,9 +474,11 @@ def _configure_parameters(self): raise ValueError(f'Unhandled target(s) ({targets[phase_name]}) for parameter {name} in ' f'phase {phase_name}. If connecting to ODE inputs in the phase, ' f'format the targets as a sequence of strings.') - tgts.append(tgt) - if not tgts: + if target_param is not None: + targets_per_phase[phase_name] = target_param + + if not targets_per_phase: # Find the reason if targets is None: reason = f'Option `targets=None` but no phase in the trajectory has a parameter named `{name}`.' @@ -447,22 +488,19 @@ def _configure_parameters(self): reason = '' raise ValueError(f'No target was found for trajectory parameter `{name}` in any phase.\n{reason}') - if options['shape'] is _unspecified: - if len(set(tgt_shapes.values())) == 1: - options['shape'] = next(iter(tgt_shapes.values())) - else: - raise ValueError(f'Parameter {name} in Trajectory {self.pathname} is connected to ' - f'targets in multiple phases that have different shapes.') + # If metadata is unspecified, use introspection to find + # it based on common values among the targets. + params_by_phase = self._get_phase_parameters() + + targets = {phase_name: phs_params[targets_per_phase[phase_name]] + for phase_name, phs_params in params_by_phase.items() + if phase_name in targets_per_phase and targets_per_phase[phase_name] in phs_params} if options['units'] is _unspecified: - tgt_units_set = set(tgt_units.values()) - if len(tgt_units_set) == 1: - options['units'] = tgt_units_set.pop() - else: - ValueError(f'Parameter {name} in Trajectory {self.pathname} is connected to ' - f'targets in multiple phases that have different units. You must ' - f'explicitly provide units for the parameter since they cannot be ' - f'inferred.') + options['units'] = _get_common_metadata(targets, metadata_key='units') + + if options['shape'] in _none_or_unspecified: + options['shape'] = _get_common_metadata(targets, metadata_key='shape') param_comp = self._get_subsystem('param_comp') param_comp.add_parameter(name, val=options['val'], shape=options['shape'], units=options['units']) @@ -478,6 +516,7 @@ def _configure_parameters(self): ref0=options['ref0'], ref=options['ref']) + tgts = [f'{phase_name}.parameters:{param_name}' for phase_name, param_name in targets_per_phase.items()] self.connect(f'parameter_vals:{name}', tgts) return promoted_inputs @@ -496,7 +535,7 @@ def _configure_phase_options_dicts(self): all_ranks = self.comm.allgather(options['shape']) for item in all_ranks: - if item not in {None, _unspecified}: + if item not in _none_or_unspecified: options['shape'] = item break else: diff --git a/dymos/transcriptions/analytic/analytic.py b/dymos/transcriptions/analytic/analytic.py index cbd649120..9380da70c 100644 --- a/dymos/transcriptions/analytic/analytic.py +++ b/dymos/transcriptions/analytic/analytic.py @@ -425,18 +425,17 @@ def get_parameter_connections(self, name, phase): if name in phase.parameter_options: options = phase.parameter_options[name] - if not options['static_target']: - src_idxs_raw = np.zeros(self.grid_data.subset_num_nodes['all'], dtype=int) - src_idxs = get_src_indices_by_row(src_idxs_raw, options['shape']) - if options['shape'] == (1,): - src_idxs = src_idxs.ravel() - else: - src_idxs_raw = np.zeros(1, dtype=int) - src_idxs = get_src_indices_by_row(src_idxs_raw, options['shape']) - src_idxs = np.squeeze(src_idxs, axis=0) - - rhs_tgts = [f'rhs.{t}' for t in options['targets']] - connection_info.append((rhs_tgts, (src_idxs,))) + for tgt in options['targets']: + if tgt in options['static_targets']: + src_idxs_raw = np.zeros(1, dtype=int) + src_idxs = get_src_indices_by_row(src_idxs_raw, options['shape']) + src_idxs = np.squeeze(src_idxs, axis=0) + else: + src_idxs_raw = np.zeros(self.grid_data.subset_num_nodes['all'], dtype=int) + src_idxs = get_src_indices_by_row(src_idxs_raw, options['shape']) + if options['shape'] == (1,): + src_idxs = src_idxs.ravel() + connection_info.append((f'rhs.{tgt}', (src_idxs,))) return connection_info diff --git a/dymos/transcriptions/common/parameter_comp.py b/dymos/transcriptions/common/parameter_comp.py index 203c4129b..b78d74f34 100644 --- a/dymos/transcriptions/common/parameter_comp.py +++ b/dymos/transcriptions/common/parameter_comp.py @@ -4,7 +4,7 @@ import numpy as np from openmdao.core.explicitcomponent import ExplicitComponent -from ...utils.misc import _unspecified +from ...utils.misc import _none_or_unspecified from ..._options import options as dymos_options @@ -133,7 +133,7 @@ def add_parameter(self, name, val=1.0, shape=None, output_name=None, _val = np.asarray(val) - if shape in {None, _unspecified}: + if shape in _none_or_unspecified: _shape = (1,) size = _val.size else: @@ -157,10 +157,7 @@ def add_parameter(self, name, val=1.0, shape=None, output_name=None, elif isinstance(output_tags, str): output_tags = [output_tags] - if np.ndim(val) == 0 or _val.shape == (1,): - in_val = np.full(_shape, _val) - else: - in_val = _val + in_val = _val out_val = np.expand_dims(in_val, axis=0) i_meta = self.add_input(name=f'parameters:{name}', val=in_val, shape=_shape, units=units, desc=desc, diff --git a/dymos/transcriptions/explicit_shooting/explicit_shooting.py b/dymos/transcriptions/explicit_shooting/explicit_shooting.py index aa244234b..df26f1209 100644 --- a/dymos/transcriptions/explicit_shooting/explicit_shooting.py +++ b/dymos/transcriptions/explicit_shooting/explicit_shooting.py @@ -13,7 +13,7 @@ from ..._options import options as dymos_options from ...utils.misc import get_rate_units, CoerceDesvar from ...utils.indexing import get_src_indices_by_row -from ...utils.introspection import get_promoted_vars, get_source_metadata, get_targets, get_target_metadata +from ...utils.introspection import get_promoted_vars, get_source_metadata, get_targets, _get_targets_metadata from ...utils.constants import INF_BOUND from ..common import TimeComp, TimeseriesOutputGroup, ControlGroup, PolynomialControlGroup, \ ParameterComp @@ -217,12 +217,13 @@ def configure_time(self, phase): phase.connect(f'integrator.{name}', [f'ode.{t}' for t in targets], src_indices=src_idxs, flat_src_indices=True if dynamic else None) - for name, targets in [('t_initial', time_options['t_initial_targets']), - ('t_duration', time_options['t_duration_targets'])]: - for t in targets: - tgt_shape, _, static_tgt = get_target_metadata(ode, name=name, - user_targets=t, - user_units=time_options['units']) + for name, tgts in [('t_initial', time_options['t_initial_targets']), + ('t_duration', time_options['t_duration_targets'])]: + + targets = _get_targets_metadata(ode, name, user_targets=tgts) + for t, meta in targets.items(): + tgt_shape = meta['shape'] + if tgt_shape == (1,): src_idxs = None flat_src_idxs = None @@ -424,14 +425,14 @@ def configure_controls(self, phase): [f'ode.{t}' for t in targets]) # Rate targets - rate_targets = get_targets(ode_inputs, control_name, options['rate_targets'], control_rates=1) + rate_targets = get_targets(ode_inputs, control_name, options['rate_targets']) if rate_targets: phase.connect(f'control_rates:{control_name}_rate', [f'ode.{t}' for t in rate_targets]) # Second time derivative targets must be specified explicitly - rate2_targets = get_targets(ode_inputs, control_name, options['rate2_targets'], control_rates=2) + rate2_targets = get_targets(ode_inputs, control_name, options['rate2_targets']) if rate2_targets: phase.connect(f'control_rates:{control_name}_rate2', @@ -518,14 +519,14 @@ def configure_polynomial_controls(self, phase): [f'ode.{t}' for t in targets]) # Rate targets - rate_targets = get_targets(ode_inputs, control_name, options['rate_targets'], control_rates=1) + rate_targets = get_targets(ode_inputs, control_name, options['rate_targets']) if rate_targets: phase.connect(f'polynomial_control_rates:{control_name}_rate', [f'ode.{t}' for t in rate_targets]) # Second time derivative targets must be specified explicitly - rate2_targets = get_targets(ode_inputs, control_name, options['rate2_targets'], control_rates=2) + rate2_targets = get_targets(ode_inputs, control_name, options['rate2_targets']) if rate2_targets: phase.connect(f'polynomial_control_rates:{control_name}_rate2', @@ -707,18 +708,19 @@ def get_parameter_connections(self, name, phase): if name in phase.parameter_options: options = phase.parameter_options[name] - if not options['static_target']: - src_idxs_raw = np.zeros(self._output_grid_data.subset_num_nodes['all'], dtype=int) - src_idxs = get_src_indices_by_row(src_idxs_raw, options['shape']) - if options['shape'] == (1,): - src_idxs = src_idxs.ravel() - else: - src_idxs_raw = np.zeros(1, dtype=int) - src_idxs = get_src_indices_by_row(src_idxs_raw, options['shape']) - src_idxs = np.squeeze(src_idxs, axis=0) + for tgt in options['targets']: + if tgt in options['static_targets']: + src_idxs_raw = np.zeros(1, dtype=int) + src_idxs = get_src_indices_by_row(src_idxs_raw, options['shape']) + src_idxs = np.squeeze(src_idxs, axis=0) + else: + src_idxs_raw = np.zeros(self._output_grid_data.subset_num_nodes['all'], dtype=int) + src_idxs = get_src_indices_by_row(src_idxs_raw, options['shape']) + if options['shape'] == (1,): + src_idxs = src_idxs.ravel() + connection_info.append((f'ode.{tgt}', (src_idxs,))) connection_info.append(([f'integrator.parameters:{name}'], None)) - connection_info.append(([f'ode.{tgt}' for tgt in options['targets']], (src_idxs,))) return connection_info diff --git a/dymos/transcriptions/explicit_shooting/ode_evaluation_group.py b/dymos/transcriptions/explicit_shooting/ode_evaluation_group.py index 5b80b73b3..186b6ca88 100644 --- a/dymos/transcriptions/explicit_shooting/ode_evaluation_group.py +++ b/dymos/transcriptions/explicit_shooting/ode_evaluation_group.py @@ -1,3 +1,4 @@ +from copy import deepcopy import numpy as np import openmdao.api as om @@ -5,10 +6,11 @@ from .state_rate_collector_comp import StateRateCollectorComp from .tau_comp import TauComp -from ...utils.introspection import get_targets, configure_controls_introspection, \ +from ...utils.introspection import configure_controls_introspection, \ configure_time_introspection, configure_parameters_introspection, \ - configure_states_discovery, configure_states_introspection, get_target_metadata -from ...utils.misc import get_rate_units + configure_states_discovery, configure_states_introspection, _get_targets_metadata, \ + _get_common_metadata, get_promoted_vars +from ...utils.misc import get_rate_units, _unspecified, _none_or_unspecified class ODEEvaluationGroup(om.Group): @@ -43,11 +45,16 @@ def __init__(self, ode_class, input_grid_data, time_options, state_options, para polynomial_control_options, ode_init_kwargs=None, vec_size=1, **kwargs): super().__init__(**kwargs) - self._state_options = state_options - self._parameter_options = parameter_options - self._time_options = time_options - self._control_options = control_options - self._polynomial_control_options = polynomial_control_options + # This component creates copies of the variable options from the phase. + # It needs to perform its own introspection with respect to its ODE instance, + # and this would override unspecified variables for parameter introspection + # at the phase level. + self._state_options = deepcopy(state_options) + self._parameter_options = deepcopy(parameter_options) + self._time_options = deepcopy(time_options) + self._control_options = deepcopy(control_options) + self._polynomial_control_options = deepcopy(polynomial_control_options) + self._control_interpolants = {} self._polynomial_control_interpolants = {} self._ode_class = ode_class @@ -196,38 +203,42 @@ def _configure_states(self): def _configure_params(self): vec_size = self._vec_size + ode_inputs = get_promoted_vars(self.ode, iotypes='input', metadata_keys=['shape', 'units', 'val', 'tags']) for name, options in self._parameter_options.items(): var_name = f'parameters:{name}' - targets = get_targets(ode=self.ode, name=name, user_targets=options['targets']) + targets = _get_targets_metadata(ode_inputs, name=name, user_targets=options['targets']) + + if options['units'] is _unspecified: + units = _get_common_metadata(targets, 'units') + else: + units = options['units'] - shape, units, static = get_target_metadata(self.ode, name=name, - user_targets=targets, - user_shape=options['shape'], - user_units=options['units'], - user_static_target=options['static_target']) - options['units'] = units - options['shape'] = shape - options['static_target'] = static + if options['shape'] in _none_or_unspecified: + shape = _get_common_metadata(targets, 'shape') + else: + shape = options['shape'] + + static_target = [tgt for tgt, meta in targets.items() if 'dymos.static_target' in meta['tags']] self._ivc.add_output(var_name, shape=shape, units=units) self.add_design_var(var_name) - if options['static_target']: - src_idxs = None - shape = None - else: - src_rows = np.zeros(vec_size, dtype=int) - src_idxs = om.slicer[src_rows, ...] - # Promote targets from the ODE for tgt in targets: + if tgt in options['static_targets']: + src_idxs = None + shape = None + else: + src_rows = np.zeros(vec_size, dtype=int) + src_idxs = om.slicer[src_rows, ...] + self.promotes('ode', inputs=[(tgt, var_name)], src_indices=src_idxs, src_shape=shape) if targets: self.set_input_defaults(name=var_name, - val=np.ones(shape), + val=1.0, units=options['units']) def _configure_controls(self): diff --git a/dymos/transcriptions/pseudospectral/gauss_lobatto.py b/dymos/transcriptions/pseudospectral/gauss_lobatto.py index 58424efc6..93e635eaf 100644 --- a/dymos/transcriptions/pseudospectral/gauss_lobatto.py +++ b/dymos/transcriptions/pseudospectral/gauss_lobatto.py @@ -150,39 +150,33 @@ def configure_controls(self, phase): disc_src_idxs = (disc_src_idxs,) col_src_idxs = (col_src_idxs,) - # Control targets are detected automatically - targets = get_targets(ode_inputs, name, options['targets']) - - if targets: + if options['targets']: phase.connect(f'control_values:{name}', - [f'rhs_disc.{t}' for t in targets], + [f'rhs_disc.{t}' for t in options['targets']], src_indices=disc_src_idxs, flat_src_indices=True) phase.connect(f'control_values:{name}', - [f'rhs_col.{t}' for t in targets], + [f'rhs_col.{t}' for t in options['targets']], src_indices=col_src_idxs, flat_src_indices=True) # Rate targets - targets = get_targets(ode_inputs, name, options['rate_targets'], control_rates=1) - - if targets: + if options['rate_targets']: phase.connect(f'control_rates:{name}_rate', - [f'rhs_disc.{t}' for t in targets], + [f'rhs_disc.{t}' for t in options['rate_targets']], src_indices=disc_src_idxs, flat_src_indices=True) phase.connect(f'control_rates:{name}_rate', - [f'rhs_col.{t}' for t in targets], + [f'rhs_col.{t}' for t in options['rate_targets']], src_indices=col_src_idxs, flat_src_indices=True) # Second time derivative targets must be specified explicitly - targets = get_targets(ode_inputs, name, options['rate2_targets'], control_rates=2) - if targets: + if options['rate2_targets']: phase.connect(f'control_rates:{name}_rate2', - [f'rhs_disc.{t}' for t in targets], + [f'rhs_disc.{t}' for t in options['rate2_targets']], src_indices=disc_src_idxs, flat_src_indices=True) phase.connect(f'control_rates:{name}_rate2', - [f'rhs_col.{t}' for t in targets], + [f'rhs_col.{t}' for t in options['rate2_targets']], src_indices=col_src_idxs, flat_src_indices=True) def configure_polynomial_controls(self, phase): @@ -213,33 +207,30 @@ def configure_polynomial_controls(self, phase): disc_src_idxs = (disc_src_idxs,) col_src_idxs = (col_src_idxs,) - targets = get_targets(ode=ode_inputs, name=name, user_targets=options['targets']) - if targets: + if options['targets']: phase.connect(f'polynomial_control_values:{name}', - [f'rhs_disc.{t}' for t in targets], + [f'rhs_disc.{t}' for t in options['targets']], src_indices=disc_src_idxs, flat_src_indices=True) phase.connect(f'polynomial_control_values:{name}', - [f'rhs_col.{t}' for t in targets], + [f'rhs_col.{t}' for t in options['targets']], src_indices=col_src_idxs, flat_src_indices=True) - targets = get_targets(ode=ode_inputs, name=name, user_targets=options['rate_targets']) - if targets: + if options['rate_targets']: phase.connect(f'polynomial_control_rates:{name}_rate', - [f'rhs_disc.{t}' for t in targets], + [f'rhs_disc.{t}' for t in options['rate_targets']], src_indices=disc_src_idxs, flat_src_indices=True) phase.connect(f'polynomial_control_rates:{name}_rate', - [f'rhs_col.{t}' for t in targets], + [f'rhs_col.{t}' for t in options['rate_targets']], src_indices=col_src_idxs, flat_src_indices=True) - targets = get_targets(ode=ode_inputs, name=name, user_targets=options['rate2_targets']) - if targets: + if options['rate2_targets']: phase.connect(f'polynomial_control_rates:{name}_rate2', - [f'rhs_disc.{t}' for t in targets], + [f'rhs_disc.{t}' for t in options['rate2_targets']], src_indices=disc_src_idxs, flat_src_indices=True) phase.connect(f'polynomial_control_rates:{name}_rate2', - [f'rhs_col.{t}' for t in targets], + [f'rhs_col.{t}' for t in options['rate2_targets']], src_indices=col_src_idxs, flat_src_indices=True) def setup_ode(self, phase): @@ -639,32 +630,30 @@ def get_parameter_connections(self, name, phase): if name in phase.parameter_options: options = phase.parameter_options[name] - targets = options['targets'] - static = options['static_target'] + # targets = options['targets'] + # static = options['static_targets'] shape = options['shape'] - if not static: - disc_rows = np.zeros(self.grid_data.subset_num_nodes['state_disc'], dtype=int) - col_rows = np.zeros(self.grid_data.subset_num_nodes['col'], dtype=int) - disc_src_idxs = get_src_indices_by_row(disc_rows, shape) - col_src_idxs = get_src_indices_by_row(col_rows, shape) - if shape == (1,): - disc_src_idxs = disc_src_idxs.ravel() - col_src_idxs = col_src_idxs.ravel() - else: - inds = np.squeeze(get_src_indices_by_row([0], shape), axis=0) - disc_src_idxs = inds - col_src_idxs = inds - - # enclose indices in tuple to ensure shaping of indices works - disc_src_idxs = (disc_src_idxs,) - col_src_idxs = (col_src_idxs,) - - rhs_disc_tgts = [f'rhs_disc.{t}' for t in targets] - connection_info.append((rhs_disc_tgts, disc_src_idxs)) - - rhs_col_tgts = [f'rhs_col.{t}' for t in targets] - connection_info.append((rhs_col_tgts, col_src_idxs)) + for tgt in options['targets']: + if tgt in options['static_targets']: + inds = np.squeeze(get_src_indices_by_row([0], shape), axis=0) + disc_src_idxs = inds + col_src_idxs = inds + else: + disc_rows = np.zeros(self.grid_data.subset_num_nodes['state_disc'], dtype=int) + col_rows = np.zeros(self.grid_data.subset_num_nodes['col'], dtype=int) + disc_src_idxs = get_src_indices_by_row(disc_rows, shape) + col_src_idxs = get_src_indices_by_row(col_rows, shape) + if shape == (1,): + disc_src_idxs = disc_src_idxs.ravel() + col_src_idxs = col_src_idxs.ravel() + + # enclose indices in tuple to ensure shaping of indices works + disc_src_idxs = (disc_src_idxs,) + col_src_idxs = (col_src_idxs,) + + connection_info.append((f'rhs_disc.{tgt}', disc_src_idxs)) + connection_info.append((f'rhs_col.{tgt}', col_src_idxs)) return connection_info diff --git a/dymos/transcriptions/pseudospectral/radau_pseudospectral.py b/dymos/transcriptions/pseudospectral/radau_pseudospectral.py index d6f672e32..a074c54ad 100644 --- a/dymos/transcriptions/pseudospectral/radau_pseudospectral.py +++ b/dymos/transcriptions/pseudospectral/radau_pseudospectral.py @@ -314,11 +314,7 @@ def _get_rate_source_path(self, state_name, nodes, phase): node_idxs = gd.subset_node_indices[nodes] elif var_type == 'parameter': rate_path = f'parameter_vals:{var}' - dynamic = not phase.parameter_options[var]['static_target'] - if dynamic: - node_idxs = np.zeros(gd.subset_num_nodes[nodes], dtype=int) - else: - node_idxs = np.zeros(1, dtype=int) + node_idxs = np.zeros(gd.subset_num_nodes[nodes], dtype=int) else: # Failed to find variable, assume it is in the ODE rate_path = f'rhs_all.{var}' @@ -469,18 +465,16 @@ def get_parameter_connections(self, name, phase): if name in phase.parameter_options: options = phase.parameter_options[name] - if not options['static_target']: - src_idxs_raw = np.zeros(self.grid_data.subset_num_nodes['all'], dtype=int) - src_idxs = get_src_indices_by_row(src_idxs_raw, options['shape']) - if options['shape'] == (1,): - src_idxs = src_idxs.ravel() - else: - src_idxs_raw = np.zeros(1, dtype=int) - src_idxs = get_src_indices_by_row(src_idxs_raw, options['shape']) - src_idxs = np.squeeze(src_idxs, axis=0) - - rhs_all_tgts = [f'rhs_all.{t}' for t in options['targets']] - connection_info.append((rhs_all_tgts, (src_idxs,))) + for tgt in options['targets']: + if tgt in options['static_targets']: + src_idxs = np.squeeze(get_src_indices_by_row([0], options['shape']), axis=0) + else: + src_idxs_raw = np.zeros(self.grid_data.subset_num_nodes['all'], dtype=int) + src_idxs = get_src_indices_by_row(src_idxs_raw, options['shape']) + if options['shape'] == (1,): + src_idxs = src_idxs.ravel() + + connection_info.append((f'rhs_all.{tgt}', (src_idxs,))) return connection_info diff --git a/dymos/transcriptions/solve_ivp/solve_ivp.py b/dymos/transcriptions/solve_ivp/solve_ivp.py index b8f733649..a8d478399 100644 --- a/dymos/transcriptions/solve_ivp/solve_ivp.py +++ b/dymos/transcriptions/solve_ivp/solve_ivp.py @@ -10,7 +10,7 @@ SolveIVPControlGroup, SolveIVPPolynomialControlGroup, SolveIVPTimeseriesOutputComp from ..common import TimeComp, TimeseriesOutputGroup from ...utils.misc import get_rate_units -from ...utils.introspection import get_promoted_vars, get_targets, get_source_metadata, get_target_metadata +from ...utils.introspection import get_promoted_vars, get_targets, get_source_metadata from ...utils.indexing import get_src_indices_by_row diff --git a/dymos/transcriptions/transcription_base.py b/dymos/transcriptions/transcription_base.py index bda172659..432466541 100644 --- a/dymos/transcriptions/transcription_base.py +++ b/dymos/transcriptions/transcription_base.py @@ -7,8 +7,8 @@ from .common import ControlGroup, PolynomialControlGroup, ParameterComp from ..utils.constants import INF_BOUND from ..utils.indexing import get_constraint_flat_idxs -from ..utils.misc import _unspecified -from ..utils.introspection import configure_states_introspection, get_promoted_vars, get_target_metadata, \ +from ..utils.misc import _none_or_unspecified +from ..utils.introspection import configure_states_introspection, get_promoted_vars, \ configure_states_discovery @@ -105,7 +105,7 @@ def configure_time(self, phase): time_options = phase.time_options # Determine the time unit. - if time_options['units'] in {None, _unspecified}: + if time_options['units'] in _none_or_unspecified: if time_options['targets']: ode = phase._get_subsystem(self._rhs_source) @@ -273,7 +273,6 @@ def configure_parameters(self, phase): for name, options in phase.parameter_options.items(): param_comp.add_parameter(name, val=options['val'], shape=options['shape'], units=options['units']) - if options['opt']: lb = -INF_BOUND if options['lower'] is None else options['lower'] ub = INF_BOUND if options['upper'] is None else options['upper'] @@ -286,11 +285,8 @@ def configure_parameters(self, phase): ref=options['ref']) for tgts, src_idxs in self.get_parameter_connections(name, phase): - if not options['static_target']: - phase.connect(f'parameter_vals:{name}', tgts, src_indices=src_idxs, - flat_src_indices=True) - else: - phase.connect(f'parameter_vals:{name}', tgts) + phase.connect(f'parameter_vals:{name}', tgts, src_indices=src_idxs, + flat_src_indices=True) def setup_states(self, phase): """ diff --git a/dymos/utils/introspection.py b/dymos/utils/introspection.py index 4e8b351f0..cd3ed766d 100644 --- a/dymos/utils/introspection.py +++ b/dymos/utils/introspection.py @@ -1,11 +1,12 @@ -from collections.abc import Iterable +from collections.abc import Iterable, Sequence import fnmatch +from numbers import Number import re import openmdao.api as om import numpy as np -from openmdao.utils.array_utils import shape_to_len -from dymos.utils.misc import _unspecified +from openmdao.utils.general_utils import ensure_compatible +from dymos.utils.misc import _unspecified, _none_or_unspecified from .._options import options as dymos_options from ..phase.options import StateOptionsDictionary, TimeseriesOutputOptionsDictionary from .misc import get_rate_units @@ -115,7 +116,7 @@ def get_promoted_vars(ode, iotypes, metadata_keys=None, get_remote=True): metadata_keys=metadata_keys).values()} -def get_targets(ode, name, user_targets, control_rates=False): +def get_targets(ode, name, user_targets): """ Return the targets of a variable in a given ODE system. @@ -135,10 +136,6 @@ def get_targets(ode, name, user_targets, control_rates=False): The name of the variable whose targets are desired. user_targets : str or None or Sequence or _unspecified Targets for the variable as given by the user. - control_rates : bool or int - If True, search for the target of the variable with the given name. If 1, search for - the first rate of the variable '{control_name}_rate', and if 2, search for the second - derivative of the variable '{control_name}_rate2'. Returns ------- @@ -156,18 +153,13 @@ def get_targets(ode, name, user_targets, control_rates=False): else: ode_inputs = get_promoted_vars(ode, iotypes=('input',)) if user_targets is _unspecified: - if name in ode_inputs and control_rates not in {1, 2}: + if name in ode_inputs: return [name] - elif control_rates == 1 and f'{name}_rate' in ode_inputs: - return [f'{name}_rate'] - elif control_rates == 2 and f'{name}_rate2' in ode_inputs: - return [f'{name}_rate2'] elif user_targets: if isinstance(user_targets, str): return [user_targets] else: return user_targets - return [] @@ -314,53 +306,68 @@ def configure_controls_introspection(control_options, ode, time_units='s'): If the control or one of its rates are connected to a variable that is tagged as static within the ODE. """ - ode_inputs = get_promoted_vars(ode, iotypes='input') + ode_inputs = get_promoted_vars(ode, iotypes='input', metadata_keys=['shape', 'units', 'val', 'tags']) for name, options in control_options.items(): + targets = _get_targets_metadata(ode_inputs, name=name, user_targets=options['targets']) - targets, shape, units, static_target = _get_targets_metadata(ode_inputs, - name=name, - user_targets=options['targets'], - user_units=options['units'], - user_shape=options['shape'], - control_rate=False) - options['targets'] = targets + options['targets'] = list(targets.keys()) if targets: - options['units'] = units - options['shape'] = shape + if options['units'] is _unspecified: + options['units'] = _get_common_metadata(targets, metadata_key='units') - if static_target: - raise ValueError(f"Control '{name}' cannot be connected to its targets because one " - f"or more targets are tagged with 'dymos.static_target'.") + if options['shape'] in _none_or_unspecified: + shape = _get_common_metadata(targets, metadata_key='shape') + if len(shape) == 1: + options['shape'] = (1,) + else: + options['shape'] = shape[1:] + + if any(['dymos.static_target' in meta['tags'] for meta in targets.values()]): + raise ValueError(f"Control '{name}' cannot be connected to its targets because one " + f"or more targets are tagged with 'dymos.static_target'.") # Now check rate targets - rate_targets, shape, units, static_target = _get_targets_metadata(ode_inputs, - name=name, - user_targets=options['rate_targets'], - user_units=options['units'], - user_shape=options['shape'], - control_rate=1) - options['rate_targets'] = rate_targets - if options['units'] is _unspecified: - options['units'] = time_units if units is None else f'{units}*{time_units}' - options['shape'] = shape - if static_target: - raise ValueError(f"Control rate of '{name}' cannot be connected to its targets " - f"because one or more targets are tagged with 'dymos.static_target'.") + rate_targets = _get_targets_metadata(ode_inputs, name=f'{name}_rate', + user_targets=options['rate_targets']) + + options['rate_targets'] = list(rate_targets.keys()) + if rate_targets: + if options['units'] is _unspecified: + rate_target_units = _get_common_metadata(rate_targets, metadata_key='units') + options['units'] = time_units if rate_target_units is None else f'{rate_target_units}*{time_units}' + + if options['shape'] in _none_or_unspecified: + shape = _get_common_metadata(rate_targets, metadata_key='shape') + if len(shape) == 1: + options['shape'] = (1,) + else: + options['shape'] = shape[1:] + + if any(['dymos.static_target' in meta['tags'] for meta in rate_targets.values()]): + raise ValueError(f"Control rate of '{name}' cannot be connected to its targets because one " + f"or more targets are tagged with 'dymos.static_target'.") # Now check rate2 targets - rate2_targets, shape, units, static_target = _get_targets_metadata(ode_inputs, - name=name, - user_targets=options['rate2_targets'], - user_units=options['units'], - user_shape=options['shape'], - control_rate=2) - options['rate2_targets'] = rate2_targets - if options['units'] is _unspecified: - options['units'] = time_units if units is None else f'{units}*{time_units}**2' - options['shape'] = shape - if static_target: - raise ValueError(f"Control rate2 of '{name}' cannot be connected to its targets " - f"because one or more targets are tagged with 'dymos.static_target'.") + rate2_targets = _get_targets_metadata(ode_inputs, name=f'{name}_rate2', + user_targets=options['rate2_targets']) + + options['rate2_targets'] = list(rate2_targets.keys()) + if rate2_targets: + if options['units'] is _unspecified: + rate2_target_units = _get_common_metadata(rate_targets, metadata_key='units') + options['units'] = f'{time_units**2}' if rate2_target_units is None \ + else f'{rate2_target_units}*{time_units}**2' + + if options['shape'] in _none_or_unspecified: + shape = _get_common_metadata(rate2_targets, metadata_key='shape') + if len(shape) == 1: + options['shape'] = (1,) + else: + options['shape'] = shape[1:] + + if any(['dymos.static_target' in meta['tags'] for meta in rate2_targets.values()]): + raise ValueError(f"Control rate2 of '{name}' cannot be connected to its targets because one " + f"or more targets are tagged with 'dymos.static_target'.") def configure_parameters_introspection(parameter_options, ode): @@ -375,21 +382,68 @@ def configure_parameters_introspection(parameter_options, ode): ode : om.System An instantiated System that serves as the ODE to which the parameters should be applied. """ - ode_inputs = get_promoted_vars(ode, iotypes='input') - for name, options in parameter_options.items(): try: - targets, shape, units, static_target = _get_targets_metadata(ode_inputs, name=name, - user_targets=options['targets'], - user_units=options['units'], - user_shape=options['shape'], - user_static_target=options['static_target']) + targets = _get_targets_metadata(ode, name=name, user_targets=options['targets']) except ValueError as e: raise ValueError(f'Parameter `{name}` has invalid target(s).\n{str(e)}') from e - options['targets'] = targets - options['units'] = units - options['shape'] = shape - options['static_target'] = static_target + + options['targets'] = list(targets.keys()) + + static_tagged_targets = {tgt for tgt, meta in targets.items() if 'dymos.static_target' in meta['tags']} + + # This is a bit of a hack. Any target with a shape of (1,) is unambiguously static. + # We may want to consider forcing users to tag these as static for dymos 2.0.0 + shape_1_targets = {tgt for tgt, meta in targets.items() if meta['shape'] == (1,)} + if options['static_targets'] is _unspecified: + options['static_targets'] = static_tagged_targets.union(shape_1_targets) + elif options['static_targets']: + options['static_targets'] = options['targets'].copy() + else: + options['static_targets'] = [] + + if static_tagged_targets and not options['static_targets']: + raise ValueError(f"Parameter `{name}` has invalid target(s).\n" + f"User has specified 'static_target = False' for parameter `{name}`,\nbut one or more " + f"targets is tagged with 'dymos.static_target':\n{static_tagged_targets}") + + if options['units'] is _unspecified: + options['units'] = _get_common_metadata(targets, metadata_key='units') + + # Check that all targets have the same shape. + tgt_shapes = {} + # First find the shapes of the static targets + for tgt, meta in targets.items(): + if tgt in options['static_targets']: + tgt_shapes[tgt] = meta['shape'] + else: + if len(meta['shape']) == 1: + tgt_shapes[tgt] = (1,) + else: + tgt_shapes[tgt] = meta['shape'][1:] + # Check that they're unique + if len(set(tgt_shapes.values())) > 1: + raise RuntimeError(f'Invalid targets for parameter `{name}`.\n' + f'Targets have multiple shapes.\n' + f'{tgt_shapes}') + elif len(set(tgt_shapes.values())) == 1: + introspected_shape = next(iter(set(tgt_shapes.values()))) + else: + introspected_shape = None + + if options['shape'] in _none_or_unspecified: + if isinstance(options['val'], Number): + options['shape'] = introspected_shape + else: + options['shape'] = np.asarray(options['val']).shape + else: + if introspected_shape is not None and options['shape'] != introspected_shape: + raise RuntimeError(f'Shape provided to parameter `{name}` differs from its targets.\n' + f'Given shape: {options["shape"]}\n' + f'Target shapes:\n' + f'{tgt_shapes}') + + options['val'], options['shape'] = ensure_compatible(name, options['val'], options['shape']) def configure_time_introspection(time_options, ode): @@ -410,33 +464,27 @@ def configure_time_introspection(time_options, ode): If time or time_phase are connected to a variable that is tagged as static within the ODE. """ - ode_inputs = get_promoted_vars(ode, 'input') + ode_inputs = get_promoted_vars(ode, 'input', metadata_keys=['shape', 'val', 'units', 'tags']) time_name = time_options['name'] t_phase_name = f'{time_name}_phase' - targets, shape, units, static_target = _get_targets_metadata(ode_inputs, - name=time_name, - user_targets=time_options['targets'], - user_units=time_options['units'], - user_shape=(1,)) + targets = _get_targets_metadata(ode_inputs, name=time_name, user_targets=time_options['targets']) time_options['targets'] = targets - time_options['units'] = units - if static_target: + if time_options['units'] is _unspecified: + time_options['units'] = _get_common_metadata(targets, 'units') + + if any(['dymos.static_target' in meta['tags'] for meta in targets.values()]): raise ValueError(f"The integration variable {time_name} cannot be connected to its targets because one " f"or more targets are tagged with 'dymos.static_target'.") # t_phase - targets, shape, units, static_target = _get_targets_metadata(ode_inputs, - name=t_phase_name, - user_targets=time_options['time_phase_targets'], - user_units=time_options['units'], - user_shape=(1,)) + targets = _get_targets_metadata(ode_inputs, name=t_phase_name, user_targets=time_options['time_phase_targets']) time_options['time_phase_targets'] = targets - if static_target: + if any(['dymos.static_target' in meta['tags'] for meta in targets.values()]): raise ValueError(f"'t_phase' cannot be connected to its targets because one " f"or more targets are tagged with 'dymos.static_target'.") @@ -470,21 +518,29 @@ def configure_states_introspection(state_options, time_options, control_options, The OpenMDAO system which provides the state rates as outputs. """ time_units = time_options['units'] - ode_inputs = get_promoted_vars(ode, 'input') - ode_outputs = get_promoted_vars(ode, 'output') + ode_inputs = get_promoted_vars(ode, 'input', metadata_keys=['units', 'shape', 'val', 'tags']) + ode_outputs = get_promoted_vars(ode, 'output', metadata_keys=['units', 'shape', 'val', 'tags']) for state_name, options in state_options.items(): # Automatically determine targets of state if left _unspecified - targets, tgt_shape, tgt_units, static_target = _get_targets_metadata(ode_inputs, - state_name, - user_targets=options['targets'], - user_units=options['units'], - user_shape=options['shape']) + targets = _get_targets_metadata(ode_inputs, state_name, + user_targets=options['targets']) - options['targets'] = targets + options['targets'] = list(targets.keys()) if targets: - options['shape'] = tgt_shape - options['units'] = tgt_units + if options['units'] is _unspecified: + options['units'] = _get_common_metadata(targets, metadata_key='units') + + if options['shape'] in _none_or_unspecified: + shape = _get_common_metadata(targets, metadata_key='shape') + if len(shape) == 1: + options['shape'] = (1,) + else: + options['shape'] = shape[1:] + + if any(['dymos.static_target' in meta['tags'] for meta in targets.values()]): + raise ValueError(f"State '{name}' cannot be connected to its targets because one " + f"or more targets are tagged with 'dymos.static_target'.") # 3. Attempt rate-source introspection rate_src = options['rate_source'] @@ -534,7 +590,7 @@ def configure_states_introspection(state_options, time_options, control_options, rate_src_shape = (1,) rate_src_units = None - if options['shape'] in {None, _unspecified}: + if options['shape'] in _none_or_unspecified: options['shape'] = rate_src_shape if options['units'] is _unspecified: @@ -584,7 +640,7 @@ def configure_analytic_states_introspection(state_options, ode): raise RuntimeError(f'ODE output {source} is tagged with `dymos.static_output` and cannot be used as a ' f'state variable in an AnalyticPhase.') - if options['shape'] in {None, _unspecified}: + if options['shape'] in _none_or_unspecified: options['shape'] = src_shape if options['units'] is _unspecified: @@ -820,7 +876,7 @@ def configure_timeseries_expr_introspection(phase): expr_reduced = expr units = output_options['units'] if output_options['units'] is not _unspecified else None - shape = output_options['shape'] if output_options['shape'] not in {_unspecified, None} else (1,) + shape = output_options['shape'] if output_options['shape'] not in _none_or_unspecified else (1,) abs_names = [x.strip() for x in re.findall(var_names_regex, expr) if not x.endswith('(') and not x.endswith(':')] @@ -942,125 +998,6 @@ def filter_outputs(patterns, sys): return results -def get_target_metadata(ode, name, user_targets=_unspecified, user_units=_unspecified, - user_shape=_unspecified, control_rate=False, user_static_target=_unspecified): - """ - Return the targets of a state variable in a given ODE system. - - If the targets of the state is _unspecified, and the state name is a top level input name - in the ODE, then the state values are automatically connected to that top-level input. - If _unspecified and not a top-level input of the ODE, no connection is made. - If targets is explicitly None, then no connection is made. - Otherwise, if the user specified some other string or sequence of strings as targets, then - those are returned. - - Parameters - ---------- - ode : om.System or dict - The OpenMDAO system which serves as the ODE for dymos, or a dictionary of inputs as returned by - utils.introspection.get_promoted_vars. If a system, it should already have had its setup and configure - methods called. - name : str - The name of the variable whose targets are desired. - user_targets : str or None or Sequence or _unspecified - Targets for the variable as given by the user. - user_units : str or None or _unspecified - Units for the variable as given by the user. - user_shape : None or Sequence or _unspecified - Shape for the variable as given by the user. - control_rate : bool - When True, check for the control rate if the name is not in the ODE. - user_static_target : bool or None or _unspecified - When False, assume the shape of the target in the ODE includes the number of nodes as the - first dimension. If True, the connecting parameter does not need to be "fanned out" to - connect to each node. If _unspecified, attempt to resolve by the presence of a tag - `dymos.static_target` on the target variable, which is the same as `static_target=True`. - - Returns - ------- - shape : tuple - The shape of the variable. If not specified, shape is taken from the ODE targets. - units : str - The units of the variable. If not specified, units are taken from the ODE targets. - static_target : bool - True if the target is static, otherwise False. - - Notes - ----- - This method requires that the ODE has run its setup and configure methods. Thus, - this method should be called from configure of some parent Group, and the ODE should - be a system within that Group. - """ - ode_inputs = ode if isinstance(ode, dict) else get_promoted_vars(ode, iotypes='input') - - if user_targets is _unspecified: - if name in ode_inputs: - targets = [name] - elif control_rate and f'{name}_rate' in ode_inputs: - targets = [f'{name}_rate'] - elif control_rate and f'{name}_rate2' in ode_inputs: - targets = [f'{name}_rate2'] - else: - targets = [] - elif user_targets: - if isinstance(user_targets, str): - targets = [user_targets] - else: - targets = user_targets - else: - targets = [] - - if user_units is _unspecified: - target_units_set = {ode_inputs[tgt]['units'] for tgt in targets} - if len(target_units_set) == 1: - units = next(iter(target_units_set)) - else: - raise ValueError(f'Unable to automatically assign units to {name}. ' - f'Targets have multiple units: {target_units_set}. ' - f'Either promote targets and use set_input_defaults to assign common ' - f'units, or explicitly provide them to {name}.') - else: - units = user_units - - # Resolve whether the targets is static or dynamic - static_target_tags = [tgt for tgt in targets if 'dymos.static_target' in ode_inputs[tgt]['tags']] - if static_target_tags: - static_target = True - if not user_static_target: - raise ValueError(f"User has specified 'static_target = False' for parameter {name}," - f"but one or more targets is tagged with " - f"'dymos.static_target': {' '.join(static_target_tags)}") - else: - if user_static_target is _unspecified: - static_target = False - else: - static_target = user_static_target - - if user_shape in {None, _unspecified}: - # Resolve target shape - target_shape_set = {ode_inputs[tgt]['shape'] for tgt in targets} - if len(target_shape_set) == 1: - shape = next(iter(target_shape_set)) - if not static_target: - if len(shape) == 1: - shape = (1,) - else: - shape = shape[1:] - elif len(target_shape_set) == 0: - raise ValueError(f'Unable to automatically assign a shape to {name}.\n' - 'Targets for this variable either do not exist or have no shape set.\n' - 'The shape for this variable must be set explicitly via the ' - '`shape=` argument.') - else: - raise ValueError(f'Unable to automatically assign a shape to {name} based on targets. ' - f'Targets have multiple shapes assigned: {target_shape_set}. ' - f'Change targets such that all have common shapes.') - else: - shape = user_shape - - return shape, units, static_target - - def configure_duration_balance_introspection(phase): """ Modify duration balance options in-place using introspection of the phase and its ODE. @@ -1177,8 +1114,7 @@ def configure_duration_balance_introspection(phase): f' has no index specified. The balance may only have shape (1,) or a single index') -def _get_targets_metadata(ode, name, user_targets=_unspecified, user_units=_unspecified, - user_shape=_unspecified, control_rate=False, user_static_target=_unspecified): +def _get_targets_metadata(ode, name, user_targets=_unspecified): """ Return the targets of a variable in a given ODE system and their metadata. @@ -1199,28 +1135,11 @@ def _get_targets_metadata(ode, name, user_targets=_unspecified, user_units=_unsp The name of the variable whose targets are desired. user_targets : str or None or Sequence or _unspecified Targets for the variable as given by the user. - user_units : str or None or _unspecified - Units for the variable as given by the user. - user_shape : None or Sequence or _unspecified - Shape for the variable as given by the user. - control_rate : bool - When True, check for the control rate if the name is not in the ODE. - user_static_target : bool or None or _unspecified - When False, assume the shape of the target in the ODE includes the number of nodes as the - first dimension. If True, the connecting parameter does not need to be "fanned out" to - connect to each node. If _unspecified, attempt to resolve by the presence of a tag - `dymos.static_target` on the target variable, which is the same as `static_target=True`. Returns ------- - targets : list - The target inputs of the variable in the ODE, as a list. - shape : tuple - The shape of the variable. If not specified, shape is taken from the ODE targets. - units : str - The units of the variable. If not specified, units are taken from the ODE targets. - static_target : bool - True if the target is static, otherwise False. + targets : dict + A dictionary of target inputs in the ODE and metadata associated with each target. Notes ----- @@ -1228,66 +1147,62 @@ def _get_targets_metadata(ode, name, user_targets=_unspecified, user_units=_unsp this method should be called from configure of some parent Group, and the ODE should be a system within that Group. """ - ode_inputs = ode if isinstance(ode, dict) else get_promoted_vars(ode, iotypes='input') - - targets = get_targets(ode_inputs, name, user_targets=user_targets, control_rates=control_rate) + ode_inputs = ode if isinstance(ode, dict) else get_promoted_vars(ode, + iotypes='input', + metadata_keys=['shape', 'units', + 'val', 'tags']) - if not targets: - return targets, user_shape, user_units, False + targets = {t: {} for t in get_targets(ode_inputs, name, user_targets=user_targets)} for tgt in targets: if tgt not in ode_inputs: raise ValueError(f"No such ODE input: '{tgt}'.") - if user_units is _unspecified: - target_units_set = {ode_inputs[tgt]['units'] for tgt in targets} - if len(target_units_set) == 1: - units = next(iter(target_units_set)) - else: - raise ValueError(f'Unable to automatically assign units to {name}. ' - f'Targets have multiple units: {target_units_set}. ' - f'Either promote targets and use set_input_defaults to assign common ' - f'units, or explicitly provide them to {name}.') - else: - units = user_units - - # Resolve whether the targets is static or dynamic - static_target_tags = [tgt for tgt in targets if 'dymos.static_target' in ode_inputs[tgt]['tags']] - if static_target_tags: - static_target = True - if not user_static_target: - raise ValueError(f"User has specified 'static_target = False' for parameter {name}," - f"but one or more targets is tagged with " - f"'dymos.static_target': {' '.join(static_target_tags)}") - else: - if user_static_target is _unspecified: - static_target = False - else: - static_target = user_static_target - - if user_shape in {None, _unspecified}: - # Resolve target shape - target_shape_set = {ode_inputs[tgt]['shape'] for tgt in targets} - if len(target_shape_set) == 1: - shape = next(iter(target_shape_set)) - if not static_target: - if len(shape) == 1: - shape = (1,) - else: - shape = shape[1:] - elif len(target_shape_set) == 0: - raise ValueError(f'Unable to automatically assign a shape to {name}.\n' - 'Targets for this variable either do not exist or have no shape set.\n' - 'The shape for this variable must be set explicitly via the ' - '`shape=` argument.') - else: - raise ValueError(f'Unable to automatically assign a shape to {name} based on targets. ' - f'Targets have multiple shapes assigned: {target_shape_set}. ' - f'Change targets such that all have common shapes.') - else: - shape = user_shape + for tgt, options in targets.items(): + options['units'] = ode_inputs[tgt]['units'] + options['shape'] = ode_inputs[tgt]['shape'] + options['val'] = ode_inputs[tgt]['val'] + options['tags'] = ode_inputs[tgt]['tags'] + + return targets + - return targets, shape, units, static_target +def _get_common_metadata(targets, metadata_key): + """ + Given a dictionary containing targets and their metadata, return the value associated + with the given metadata key if that value is common to all targets, otherwise raise an + Exception. + + Parameters + ---------- + targets : dict + A dictionary of targets and their metadata which must include the desired metadata key. + metadata_key : str + The metadata key desired. + + Returns + ------- + object + The common metadata value shared by all targets. + + Raises + ------ + ValueError + ValueError is raised if the targets do not all have the same metadata value. + """ + meta_set = {meta[metadata_key] for meta in targets.values()} + + if len(meta_set) == 1: + return next(iter(meta_set)) + elif len(meta_set) == 0: + raise RuntimeError(f'Unable to automatically assign {metadata_key} based on targets. \n' + f'No targets were found.') + else: + err_dict = {tgt: meta[metadata_key] for tgt, meta in targets.items()} + raise RuntimeError(f'Unable to automatically assign {metadata_key} based on targets.\n' + f'Targets have multiple {metadata_key} assigned:\n{err_dict}.\n' + f'Either promote targets and use set_input_defaults to assign common\n' + f'{metadata_key}, or explicitly provide {metadata_key} to the variable.') def get_source_metadata(ode, src, user_units=_unspecified, user_shape=_unspecified): @@ -1331,12 +1246,12 @@ def get_source_metadata(ode, src, user_units=_unspecified, user_shape=_unspecifi if src not in ode_outputs: raise ValueError(f"Unable to find the source '{src}' in the ODE.") - if user_units in {None, _unspecified}: + if user_units in _none_or_unspecified: meta['units'] = ode_outputs[src]['units'] else: meta['units'] = user_units - if user_shape in {None, _unspecified}: + if user_shape in _none_or_unspecified: ode_shape = ode_outputs[src]['shape'] meta['shape'] = (1,) if len(ode_shape) == 1 else ode_shape[1:] else: diff --git a/dymos/utils/misc.py b/dymos/utils/misc.py index 30b3f9f78..1576205c2 100644 --- a/dymos/utils/misc.py +++ b/dymos/utils/misc.py @@ -9,6 +9,7 @@ # unique object to check if default is given (when None is an allowed value) _unspecified = _ReprClass("unspecified") +_none_or_unspecified = {None, _unspecified} def get_rate_units(units, time_units, deriv=1): diff --git a/dymos/visualization/linkage/test/test_linkage_report.py b/dymos/visualization/linkage/test/test_linkage_report.py index 0f01b0bc7..467af31cf 100644 --- a/dymos/visualization/linkage/test/test_linkage_report.py +++ b/dymos/visualization/linkage/test/test_linkage_report.py @@ -101,72 +101,72 @@ def test_model_data(self): targets={'br_to_v1': ['m'], 'v1_to_vr': ['m'], 'rto': ['m'], 'rotate': ['m'], 'climb': ['m']}) - traj.add_parameter('T_nominal', val=27000 * 2, opt=False, units='lbf', static_target=True, + traj.add_parameter('T_nominal', val=27000 * 2, opt=False, units='lbf', static_targets=True, desc='nominal aircraft thrust', targets={'br_to_v1': ['T']}) - traj.add_parameter('T_engine_out', val=27000, opt=False, units='lbf', static_target=True, + traj.add_parameter('T_engine_out', val=27000, opt=False, units='lbf', static_targets=True, desc='thrust under a single engine', targets={'v1_to_vr': ['T'], 'rotate': ['T'], 'climb': ['T']}) - traj.add_parameter('T_shutdown', val=0.0, opt=False, units='lbf', static_target=True, + traj.add_parameter('T_shutdown', val=0.0, opt=False, units='lbf', static_targets=True, desc='thrust when engines are shut down for rejected takeoff', targets={'rto': ['T']}) - traj.add_parameter('mu_r_nominal', val=0.03, opt=False, units=None, static_target=True, + traj.add_parameter('mu_r_nominal', val=0.03, opt=False, units=None, static_targets=True, desc='nominal runway friction coeffcient', targets={'br_to_v1': ['mu_r'], 'v1_to_vr': ['mu_r'], 'rotate': ['mu_r']}) - traj.add_parameter('mu_r_braking', val=0.3, opt=False, units=None, static_target=True, + traj.add_parameter('mu_r_braking', val=0.3, opt=False, units=None, static_targets=True, desc='runway friction coefficient under braking', targets={'rto': ['mu_r']}) - traj.add_parameter('h_runway', val=0., opt=False, units='ft', static_target=False, + traj.add_parameter('h_runway', val=0., opt=False, units='ft', static_targets=False, desc='runway altitude', targets={'br_to_v1': ['h'], 'v1_to_vr': ['h'], 'rto': ['h'], 'rotate': ['h']}) - traj.add_parameter('rho', val=1.225, opt=False, units='kg/m**3', static_target=True, + traj.add_parameter('rho', val=1.225, opt=False, units='kg/m**3', static_targets=True, desc='atmospheric density', targets={'br_to_v1': ['rho'], 'v1_to_vr': ['rho'], 'rto': ['rho'], 'rotate': ['rho']}) - traj.add_parameter('S', val=124.7, opt=False, units='m**2', static_target=True, + traj.add_parameter('S', val=124.7, opt=False, units='m**2', static_targets=True, desc='aerodynamic reference area', targets={'br_to_v1': ['S'], 'v1_to_vr': ['S'], 'rto': ['S'], 'rotate': ['S'], 'climb': ['S']}) - traj.add_parameter('CD0', val=0.03, opt=False, units=None, static_target=True, + traj.add_parameter('CD0', val=0.03, opt=False, units=None, static_targets=True, desc='zero-lift drag coefficient', targets={f'{phase}': ['CD0'] for phase in ['br_to_v1', 'v1_to_vr', 'rto', 'rotate', 'climb']}) - traj.add_parameter('AR', val=9.45, opt=False, units=None, static_target=True, + traj.add_parameter('AR', val=9.45, opt=False, units=None, static_targets=True, desc='wing aspect ratio', targets={f'{phase}': ['AR'] for phase in ['br_to_v1', 'v1_to_vr', 'rto', 'rotate', 'climb']}) - traj.add_parameter('e', val=801, opt=False, units=None, static_target=True, + traj.add_parameter('e', val=801, opt=False, units=None, static_targets=True, desc='Oswald span efficiency factor', targets={f'{phase}': ['e'] for phase in ['br_to_v1', 'v1_to_vr', 'rto', 'rotate', 'climb']}) - traj.add_parameter('span', val=35.7, opt=False, units='m', static_target=True, + traj.add_parameter('span', val=35.7, opt=False, units='m', static_targets=True, desc='wingspan', targets={f'{phase}': ['span'] for phase in ['br_to_v1', 'v1_to_vr', 'rto', 'rotate', 'climb']}) - traj.add_parameter('h_w', val=1.0, opt=False, units='m', static_target=True, + traj.add_parameter('h_w', val=1.0, opt=False, units='m', static_targets=True, desc='height of wing above CG', targets={f'{phase}': ['h_w'] for phase in ['br_to_v1', 'v1_to_vr', 'rto', 'rotate', 'climb']}) - traj.add_parameter('CL0', val=0.5, opt=False, units=None, static_target=True, + traj.add_parameter('CL0', val=0.5, opt=False, units=None, static_targets=True, desc='zero-alpha lift coefficient', targets={f'{phase}': ['CL0'] for phase in ['br_to_v1', 'v1_to_vr', 'rto', 'rotate', 'climb']}) - traj.add_parameter('CL_max', val=2.0, opt=False, units=None, static_target=True, + traj.add_parameter('CL_max', val=2.0, opt=False, units=None, static_targets=True, desc='maximum lift coefficient for linear fit', targets={f'{phase}': ['CL_max'] for phase in ['br_to_v1', 'v1_to_vr', 'rto', 'rotate', 'climb']})