Repository: OpenMDAO/dymos

Issue:
Parameter issues.
### Description

1. Parameters should allow for a mixture of static and dynamic targets.
2. Trajectory parameters do not currently get sufficient introspected data from phases when under MPI.

### Example

```python
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', types=dict, 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']

        foo_tags = ['dymos.static_target'] if 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 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 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()
```


### Dymos Version

1.8.1-dev

### Relevant environment information

_No response_

Discussion:
coveralls: 
[![Coverage Status](https://coveralls.io/builds/61984239/badge)](https://coveralls.io/builds/61984239)

coverage: 92.682% (+0.1%) from 92.569% when pulling **651f199f9d04968a4467e5e5551639b1213e9ea0 on robfalck:traj_params_update** into **1aca42a4fedcca6e493c5516fadc2d78e4709797 on OpenMDAO:master**.


Patch:
diff --git a/dymos/grid_refinement/grid_refinement_ode_system.py b/dymos/grid_refinement/grid_refinement_ode_system.py
index 3fb76b90..fd86ec72 100644
--- a/dymos/grid_refinement/grid_refinement_ode_system.py
+++ b/dymos/grid_refinement/grid_refinement_ode_system.py
@@ -3,7 +3,7 @@ import openmdao.api as om
 
 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 @@ class GridRefinementODESystem(om.Group):
 
         # 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 a8a0e2bf..b675fe96 100644
--- a/dymos/grid_refinement/hp_adaptive/hp_adaptive.py
+++ b/dymos/grid_refinement/hp_adaptive/hp_adaptive.py
@@ -102,7 +102,6 @@ class HPAdaptive:
         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 48fe7cea..6f343806 100644
--- a/dymos/phase/options.py
+++ b/dymos/phase/options.py
@@ -263,15 +263,25 @@ class ParameterOptionsDictionary(om.OptionsDictionary):
                                  "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 bb32cb82..744e54a6 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 @@ class Phase(om.Group):
                       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 @@ class Phase(om.Group):
             (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 @@ class Phase(om.Group):
             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 @@ class Phase(om.Group):
             (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 @@ class Phase(om.Group):
                 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 @@ class Phase(om.Group):
 
             # 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/trajectory/trajectory.py b/dymos/trajectory/trajectory.py
index 62de5b65..46ea5469 100644
--- a/dymos/trajectory/trajectory.py
+++ b/dymos/trajectory/trajectory.py
@@ -20,8 +20,8 @@ from .phase_linkage_comp import PhaseLinkageComp
 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 @@ class Trajectory(om.Group):
     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 @@ class Trajectory(om.Group):
         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 @@ class Trajectory(om.Group):
 
         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 @@ class Trajectory(om.Group):
         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 @@ class Trajectory(om.Group):
                                    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 @@ class Trajectory(om.Group):
                         # 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 @@ class Trajectory(om.Group):
         """
         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 @@ class Trajectory(om.Group):
             # 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 @@ class Trajectory(om.Group):
                 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 @@ class Trajectory(om.Group):
                     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 @@ class Trajectory(om.Group):
                     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 @@ class Trajectory(om.Group):
                     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 @@ class Trajectory(om.Group):
                                     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 @@ class Trajectory(om.Group):
 
                     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 cbd64912..9380da70 100644
--- a/dymos/transcriptions/analytic/analytic.py
+++ b/dymos/transcriptions/analytic/analytic.py
@@ -425,18 +425,17 @@ class Analytic(TranscriptionBase):
 
         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 203c4129..b78d74f3 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 @@ class ParameterComp(ExplicitComponent):
 
         _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 @@ class ParameterComp(ExplicitComponent):
         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 aa244234..df26f120 100644
--- a/dymos/transcriptions/explicit_shooting/explicit_shooting.py
+++ b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
@@ -13,7 +13,7 @@ from .ode_integration_comp import ODEIntegrationComp
 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 @@ class ExplicitShooting(TranscriptionBase):
                 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 @@ class ExplicitShooting(TranscriptionBase):
                               [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 @@ class ExplicitShooting(TranscriptionBase):
                               [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 @@ class ExplicitShooting(TranscriptionBase):
 
         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 5b80b73b..186b6ca8 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 .vandermonde_control_interp_comp import VandermondeControlInterpComp
 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 @@ class ODEEvaluationGroup(om.Group):
                  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 @@ class ODEEvaluationGroup(om.Group):
 
     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 58424efc..93e635ea 100644
--- a/dymos/transcriptions/pseudospectral/gauss_lobatto.py
+++ b/dymos/transcriptions/pseudospectral/gauss_lobatto.py
@@ -150,39 +150,33 @@ class GaussLobatto(PseudospectralBase):
             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 @@ class GaussLobatto(PseudospectralBase):
             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 @@ class GaussLobatto(PseudospectralBase):
 
         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 d6f672e3..a074c54a 100644
--- a/dymos/transcriptions/pseudospectral/radau_pseudospectral.py
+++ b/dymos/transcriptions/pseudospectral/radau_pseudospectral.py
@@ -314,11 +314,7 @@ class Radau(PseudospectralBase):
             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 @@ class Radau(PseudospectralBase):
 
         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 b8f73364..a8d47839 100644
--- a/dymos/transcriptions/solve_ivp/solve_ivp.py
+++ b/dymos/transcriptions/solve_ivp/solve_ivp.py
@@ -10,7 +10,7 @@ from .components import SegmentSimulationComp, SegmentStateMuxComp, \
     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 bda17265..43246654 100644
--- a/dymos/transcriptions/transcription_base.py
+++ b/dymos/transcriptions/transcription_base.py
@@ -7,8 +7,8 @@ import openmdao.api as om
 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 @@ class TranscriptionBase(object):
         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 @@ class TranscriptionBase(object):
 
             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 @@ class TranscriptionBase(object):
                                          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 4e8b351f..cd3ed766 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=<tuple>` 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=<tuple>` 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 30b3f9f7..1576205c 100644
--- a/dymos/utils/misc.py
+++ b/dymos/utils/misc.py
@@ -9,6 +9,7 @@ from openmdao.core.constants import _ReprClass
 
 # 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):


Test patch:
diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_static_gravity.py b/dymos/examples/brachistochrone/test/test_brachistochrone_static_gravity.py
index 7ab79e96..e02913d2 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 @@ class TestBrachistochroneStaticGravity(unittest.TestCase):
 
         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 00000000..970fa801
--- /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 340c3185..f447bb77 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 @@ class CannonballODEVectorCD(om.ExplicitComponent):
 @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 00000000..d717576b
--- /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/phase/test/test_input_parameter_connections.py b/dymos/phase/test/test_input_parameter_connections.py
index 99dbbc16..34cd53b6 100644
--- a/dymos/phase/test/test_input_parameter_connections.py
+++ b/dymos/phase/test/test_input_parameter_connections.py
@@ -105,11 +105,12 @@ class TestStaticParameters(unittest.TestCase):
             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 @@ class TestStaticParameters(unittest.TestCase):
             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 4efe25c4..d25bf9ec 100644
--- a/dymos/phase/test/test_simulate.py
+++ b/dymos/phase/test/test_simulate.py
@@ -4,6 +4,7 @@ import numpy as np
 
 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 @@ class TestODE(om.Group):
     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 @@ class TestSimulateShapedParams(unittest.TestCase):
 
         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 @@ class TestSimulateShapedParams(unittest.TestCase):
 
         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 458cc313..0661627a 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 95fd0a95..f7b1b9e3 100755
--- a/dymos/test/test_run_problem.py
+++ b/dymos/test/test_run_problem.py
@@ -336,7 +336,6 @@ class TestRunProblem(unittest.TestCase):
                          '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/visualization/linkage/test/test_linkage_report.py b/dymos/visualization/linkage/test/test_linkage_report.py
index 0f01b0bc..467af31c 100644
--- a/dymos/visualization/linkage/test/test_linkage_report.py
+++ b/dymos/visualization/linkage/test/test_linkage_report.py
@@ -101,72 +101,72 @@ class TestLinkageReport(unittest.TestCase):
                            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']})