From 3e55e8a79023ebefdaa8476cba5ed5b54af92873 Mon Sep 17 00:00:00 2001 From: Rob Falck Date: Mon, 14 Aug 2023 16:06:22 -0400 Subject: [PATCH] Better error messages with specific tests. Allow corner case where parameter has no targets since it may be used as a rate source. --- .../test_phase_param_introspection_failure.py | 480 ++++++++++++++++++ .../test_traj_param_static_and_dynamic.py | 220 ++++++++ dymos/utils/introspection.py | 61 ++- 3 files changed, 736 insertions(+), 25 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_phase_param_introspection_failure.py b/dymos/examples/brachistochrone/test/test_phase_param_introspection_failure.py new file mode 100644 index 000000000..5e76eea16 --- /dev/null +++ b/dymos/examples/brachistochrone/test/test_phase_param_introspection_failure.py @@ -0,0 +1,480 @@ +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) + + 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'] + + 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}, + 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'] + + 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}, + 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'] + + 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}, + 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)) + + +if __name__ == '__main__': # pragma: no cover + unittest.main() 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/utils/introspection.py b/dymos/utils/introspection.py index 76486fdff..1dd35e284 100644 --- a/dymos/utils/introspection.py +++ b/dymos/utils/introspection.py @@ -411,29 +411,40 @@ def configure_parameters_introspection(parameter_options, ode): if options['units'] is _unspecified: options['units'] = _get_common_metadata(targets, metadata_key='units') + # Check that all targets have the same shape. + static_shapes = {} + dynamic_shapes = {} + # First find the shapes of the static targets + for tgt, meta in targets.items(): + if tgt in options['static_targets']: + static_shapes[tgt] = meta['shape'] + else: + if len(meta['shape']) == 1: + dynamic_shapes[tgt] = (1,) + else: + dynamic_shapes[tgt] = meta['shape'][1:] + all_shapes = {**dynamic_shapes, **static_shapes} + # Check that they're unique + if len(set(all_shapes.values())) > 1: + raise RuntimeError(f'Invalid targets for parameter `{name}`.\n' + f'Targets have multiple shapes.\n' + f'{all_shapes}') + elif len(set(all_shapes.values())) == 1: + introspected_shape = next(iter(set(all_shapes.values()))) + else: + introspected_shape = None + if options['shape'] in {_unspecified, None}: if isinstance(options['val'], Number): - static_shapes = {} - dynamic_shapes = {} - # First find the shapes of the static targets - for tgt, meta in targets.items(): - if tgt in options['static_targets']: - static_shapes[tgt] = meta['shape'] - else: - if len(meta['shape']) == 1: - dynamic_shapes[tgt] = (1,) - else: - dynamic_shapes[tgt] = meta['shape'][1:] - all_shapes = {**dynamic_shapes, **static_shapes} - # Check that they're unique - if len(set(all_shapes.values())) != 1: - raise RuntimeError(f'Unable to obtain shape of parameter {name} via introspection.\n' - f'Targets have multiple shapes.\n' - f'{all_shapes}') - else: - options['shape'] = next(iter(set(all_shapes.values()))) + 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'{all_shapes}') options['val'], options['shape'] = ensure_compatible(name, options['val'], options['shape']) @@ -1187,14 +1198,14 @@ def _get_common_metadata(targets, metadata_key): if len(meta_set) == 1: return next(iter(meta_set)) elif len(meta_set) == 0: - raise ValueError(f'Unable to automatically assign {metadata_key} based on targets. \n' - f'No targets were found.') + 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 ValueError(f'Unable to automatically assign {metadata_key} based on targets. \n' - f'Targets have multiple {metadata_key} assigned: {err_dict}. \n' - f'Either promote targets and use set_input_defaults to assign common ' - f'{metadata_key}, or explicitly provide {metadata_key} to the variable.') + 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):