diff --git a/dymos/examples/brachistochrone/test/test_phase_param_introspection_failure.py b/dymos/examples/brachistochrone/test/test_phase_param_introspection_failure.py
new file mode 100644
index 000000000..5e76eea16
--- /dev/null
+++ b/dymos/examples/brachistochrone/test/test_phase_param_introspection_failure.py
@@ -0,0 +1,480 @@
+import unittest
+import numpy as np
+import openmdao.api as om
+import dymos as dm
+
+from openmdao.utils.assert_utils import assert_near_equal
+from openmdao.utils.testing_utils import use_tempdirs
+from openmdao.utils.units import convert_units
+
+
+class _TestEOM(om.Group):
+
+    def initialize(self):
+        self.options.declare('num_nodes', types=(int,))
+        self.options.declare('foo_units', allow_none=True, default=None)
+        self.options.declare('foo_shapes', allow_none=True, default=None)
+
+    def setup(self):
+        num_nodes = self.options['num_nodes']
+
+        foo_shape = (num_nodes,) if self.options['foo_shapes'] is None \
+            else (num_nodes,) + self.options['foo_shapes']['vdot_comp']
+
+        foo_unit = 'kg' if self.options['foo_units'] is None else self.options['foo_units']['vdot_comp']
+
+        vdot_comp = om.ExecComp(['vdot = g * cos(theta)',
+                                 'bar = foo'],
+                                vdot={'shape': (num_nodes,), 'units': 'm/s**2'},
+                                g={'val': 9.80665, 'units': 'm/s**2'},
+                                theta={'shape': (num_nodes,), 'units': 'rad'},
+                                foo={'shape': foo_shape, 'units': foo_unit},
+                                bar={'shape': foo_shape, 'units': foo_unit})
+
+        foo_shape = (num_nodes,) if self.options['foo_shapes'] is None \
+            else (num_nodes,) + self.options['foo_shapes']['xdot_comp']
+
+        foo_unit = 'kg' if self.options['foo_units'] is None else self.options['foo_units']['xdot_comp']
+
+        xdot_comp = om.ExecComp(['xdot = v * sin(theta)',
+                                 'bar = foo'],
+                                xdot={'shape': (num_nodes,), 'units': 'm/s'},
+                                v={'shape': (num_nodes,), 'units': 'm/s'},
+                                theta={'shape': (num_nodes,), 'units': 'rad'},
+                                foo={'shape': foo_shape, 'units': foo_unit},
+                                bar={'shape': foo_shape, 'units': foo_unit})
+
+        foo_shape = (num_nodes,) if self.options['foo_shapes'] is None \
+            else (num_nodes,) + self.options['foo_shapes']['ydot_comp']
+
+        foo_unit = 'kg' if self.options['foo_units'] is None else self.options['foo_units']['ydot_comp']
+
+        ydot_comp = om.ExecComp(['ydot = -v * cos(theta)',
+                                 'bar = foo'],
+                                ydot={'shape': (num_nodes,), 'units': 'm/s'},
+                                v={'shape': (num_nodes,), 'units': 'm/s'},
+                                theta={'shape': (num_nodes,), 'units': 'rad'},
+                                foo={'shape': foo_shape, 'units': foo_unit},
+                                bar={'shape': foo_shape, 'units': foo_unit})
+
+        self.add_subsystem('vdot_comp', vdot_comp)
+        self.add_subsystem('xdot_comp', xdot_comp)
+        self.add_subsystem('ydot_comp', ydot_comp)
+
+
+@use_tempdirs
+class TestParameterShapes(unittest.TestCase):
+
+    def test_valid_parameters(self):
+        import numpy as np
+        import openmdao.api as om
+        import dymos as dm
+
+        #
+        # Define the OpenMDAO problem
+        #
+        p = om.Problem(model=om.Group())
+
+        #
+        # Define a Trajectory object
+        #
+        traj = dm.Trajectory()
+
+        p.model.add_subsystem('traj', subsys=traj)
+
+        #
+        # Define a Dymos Phase object with GaussLobatto Transcription
+        #
+        phase = dm.Phase(ode_class=_TestEOM,
+                         transcription=dm.GaussLobatto(num_segments=10, order=3))
+
+        traj.add_phase(name='phase0', phase=phase)
+
+        #
+        # Set the time options
+        # Time has no targets in our ODE.
+        # We fix the initial time so that it is not a design variable in the optimization.
+        # The duration of the phase is allowed to be optimized, but is bounded on [0.5, 10].
+        #
+        phase.set_time_options(fix_initial=True, duration_bounds=(0.5, 10.0), units='s')
+
+        #
+        # Set the time options
+        # Initial values of positions and velocity are all fixed.
+        # The final value of position are fixed, but the final velocity is a free variable.
+        # The equations of motion are not functions of position, so 'x' and 'y' have no targets.
+        # The rate source points to the output in the ODE which provides the time derivative of the
+        # given state.
+        phase.add_state('x', fix_initial=True, fix_final=True, units='m',
+                        rate_source='xdot_comp.xdot')
+        phase.add_state('y', fix_initial=True, fix_final=True, units='m',
+                        rate_source='ydot_comp.ydot')
+        phase.add_state('v', fix_initial=True, fix_final=False, units='m/s',
+                        rate_source='vdot_comp.vdot', targets=['xdot_comp.v', 'ydot_comp.v'])
+
+        # Define theta as a control.
+        phase.add_control(name='theta', units='rad', lower=0, upper=np.pi,
+                          targets=['xdot_comp.theta', 'ydot_comp.theta', 'vdot_comp.theta'])
+
+        phase.add_parameter('foo',
+                            opt=False,
+                            targets=['xdot_comp.foo', 'ydot_comp.foo', 'vdot_comp.foo'])
+
+        # Minimize final time.
+        phase.add_objective('time', loc='final')
+
+        # Set the driver.
+        p.driver = om.ScipyOptimizeDriver()
+
+        # Allow OpenMDAO to automatically determine our sparsity pattern.
+        # Doing so can significant speed up the execution of Dymos.
+        p.driver.declare_coloring()
+
+        # Setup the problem
+        p.setup(check=True)
+
+        # Now that the OpenMDAO problem is setup, we can set the values of the states.
+        p.set_val('traj.phase0.t_initial', 0.0, units='s')
+        p.set_val('traj.phase0.t_duration', 2.0, units='s')
+        p.set_val('traj.phase0.states:x', phase.interp('x', [0, 10]), units='m')
+        p.set_val('traj.phase0.states:y', phase.interp('y', [10, 5]), units='m')
+        p.set_val('traj.phase0.states:v', phase.interp('v', [0, 5]), units='m/s')
+        p.set_val('traj.phase0.controls:theta', phase.interp('theta', [90, 90]), units='deg')
+        p.set_val('traj.phase0.parameters:foo', 5.0)
+
+        # Run the driver to solve the problem
+        p.run_driver()
+
+        self.assertEqual((1,), phase.parameter_options['foo']['shape'])
+        assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-5)
+        assert_near_equal(p.get_val('traj.phase0.parameter_vals:foo')[-1], 5.0, tolerance=1.0E-5)
+        assert_near_equal(p.get_val('traj.phase0.rhs_col.xdot_comp.foo'), 5.0*np.ones(10,), tolerance=1.0E-5)
+        assert_near_equal(p.get_val('traj.phase0.rhs_col.ydot_comp.foo'), 5.0*np.ones(10,), tolerance=1.0E-5)
+        assert_near_equal(p.get_val('traj.phase0.rhs_col.vdot_comp.foo'), 5.0*np.ones(10,), tolerance=1.0E-5)
+        assert_near_equal(p.get_val('traj.phase0.rhs_disc.xdot_comp.foo'), 5.0*np.ones(20,), tolerance=1.0E-5)
+        assert_near_equal(p.get_val('traj.phase0.rhs_disc.ydot_comp.foo'), 5.0*np.ones(20,), tolerance=1.0E-5)
+        assert_near_equal(p.get_val('traj.phase0.rhs_disc.vdot_comp.foo'), 5.0*np.ones(20,), tolerance=1.0E-5)
+
+    def test_invalid_params_different_target_shapes(self):
+        #
+        # Define the OpenMDAO problem
+        #
+        p = om.Problem(model=om.Group())
+
+        #
+        # Define a Trajectory object
+        #
+        traj = dm.Trajectory()
+
+        p.model.add_subsystem('traj', subsys=traj)
+
+        #
+        # Define a Dymos Phase object with GaussLobatto Transcription
+        #
+        phase = dm.Phase(ode_class=_TestEOM,
+                         transcription=dm.GaussLobatto(num_segments=10, order=3),
+                         ode_init_kwargs={'foo_shapes': {'xdot_comp': (2,),
+                                                         'ydot_comp': (2,),
+                                                         'vdot_comp': (2,)}})
+
+        traj.add_phase(name='phase0', phase=phase)
+
+        #
+        # Set the time options
+        # Time has no targets in our ODE.
+        # We fix the initial time so that it is not a design variable in the optimization.
+        # The duration of the phase is allowed to be optimized, but is bounded on [0.5, 10].
+        #
+        phase.set_time_options(fix_initial=True, duration_bounds=(0.5, 10.0), units='s')
+
+        #
+        # Set the time options
+        # Initial values of positions and velocity are all fixed.
+        # The final value of position are fixed, but the final velocity is a free variable.
+        # The equations of motion are not functions of position, so 'x' and 'y' have no targets.
+        # The rate source points to the output in the ODE which provides the time derivative of the
+        # given state.
+        phase.add_state('x', fix_initial=True, fix_final=True, units='m',
+                        rate_source='xdot_comp.xdot')
+        phase.add_state('y', fix_initial=True, fix_final=True, units='m',
+                        rate_source='ydot_comp.ydot')
+        phase.add_state('v', fix_initial=True, fix_final=False, units='m/s',
+                        rate_source='vdot_comp.vdot', targets=['xdot_comp.v', 'ydot_comp.v'])
+
+        # Define theta as a control.
+        phase.add_control(name='theta', units='rad', lower=0, upper=np.pi,
+                          targets=['xdot_comp.theta', 'ydot_comp.theta', 'vdot_comp.theta'])
+
+        phase.add_parameter('foo', shape=(1,),
+                            opt=False,
+                            targets=['xdot_comp.foo', 'ydot_comp.foo', 'vdot_comp.foo'])
+
+        # Minimize final time.
+        phase.add_objective('time', loc='final')
+
+        # Set the driver.
+        p.driver = om.ScipyOptimizeDriver()
+
+        # Allow OpenMDAO to automatically determine our sparsity pattern.
+        # Doing so can significant speed up the execution of Dymos.
+        p.driver.declare_coloring()
+
+        # Setup the problem
+        with self.assertRaises(RuntimeError) as e:
+            p.setup(check=True)
+
+        expected = ("Shape provided to parameter `foo` differs from its targets.\n"
+                    "Given shape: (1,)\n"
+                    "Target shapes:\n"
+                    "{'xdot_comp.foo': (2,), 'ydot_comp.foo': (2,), 'vdot_comp.foo': (2,)}")
+        self.assertEqual(expected, str(e.exception))
+
+    def test_invalid_params_different_target_shapes_introspection_failure(self):
+        #
+        # Define the OpenMDAO problem
+        #
+        p = om.Problem(model=om.Group())
+
+        #
+        # Define a Trajectory object
+        #
+        traj = dm.Trajectory()
+
+        p.model.add_subsystem('traj', subsys=traj)
+
+        #
+        # Define a Dymos Phase object with GaussLobatto Transcription
+        #
+        phase = dm.Phase(ode_class=_TestEOM,
+                         transcription=dm.GaussLobatto(num_segments=10, order=3),
+                         ode_init_kwargs={'foo_shapes': {'xdot_comp': (1,),
+                                                         'ydot_comp': (2,),
+                                                         'vdot_comp': (3,)}})
+
+        traj.add_phase(name='phase0', phase=phase)
+
+        #
+        # Set the time options
+        # Time has no targets in our ODE.
+        # We fix the initial time so that it is not a design variable in the optimization.
+        # The duration of the phase is allowed to be optimized, but is bounded on [0.5, 10].
+        #
+        phase.set_time_options(fix_initial=True, duration_bounds=(0.5, 10.0), units='s')
+
+        #
+        # Set the time options
+        # Initial values of positions and velocity are all fixed.
+        # The final value of position are fixed, but the final velocity is a free variable.
+        # The equations of motion are not functions of position, so 'x' and 'y' have no targets.
+        # The rate source points to the output in the ODE which provides the time derivative of the
+        # given state.
+        phase.add_state('x', fix_initial=True, fix_final=True, units='m',
+                        rate_source='xdot_comp.xdot')
+        phase.add_state('y', fix_initial=True, fix_final=True, units='m',
+                        rate_source='ydot_comp.ydot')
+        phase.add_state('v', fix_initial=True, fix_final=False, units='m/s',
+                        rate_source='vdot_comp.vdot', targets=['xdot_comp.v', 'ydot_comp.v'])
+
+        # Define theta as a control.
+        phase.add_control(name='theta', units='rad', lower=0, upper=np.pi,
+                          targets=['xdot_comp.theta', 'ydot_comp.theta', 'vdot_comp.theta'])
+
+        phase.add_parameter('foo',
+                            opt=False,
+                            targets=['xdot_comp.foo', 'ydot_comp.foo', 'vdot_comp.foo'])
+
+        # Minimize final time.
+        phase.add_objective('time', loc='final')
+
+        # Set the driver.
+        p.driver = om.ScipyOptimizeDriver()
+
+        # Allow OpenMDAO to automatically determine our sparsity pattern.
+        # Doing so can significant speed up the execution of Dymos.
+        p.driver.declare_coloring()
+
+        # Setup the problem
+        with self.assertRaises(RuntimeError) as e:
+            p.setup(check=True)
+
+        expected = ('Invalid targets for parameter `foo`.\n'
+                    'Targets have multiple shapes.\n'
+                    "{'xdot_comp.foo': (1,), 'ydot_comp.foo': (2,), 'vdot_comp.foo': (3,)}")
+        self.assertEqual(expected, str(e.exception))
+
+
+@use_tempdirs
+class TestParameterUnits(unittest.TestCase):
+
+    def test_valid_parameters(self):
+        import numpy as np
+        import openmdao.api as om
+        import dymos as dm
+
+        #
+        # Define the OpenMDAO problem
+        #
+        p = om.Problem(model=om.Group())
+
+        #
+        # Define a Trajectory object
+        #
+        traj = dm.Trajectory()
+
+        p.model.add_subsystem('traj', subsys=traj)
+
+        #
+        # Define a Dymos Phase object with GaussLobatto Transcription
+        #
+        phase = dm.Phase(ode_class=_TestEOM,
+                         transcription=dm.GaussLobatto(num_segments=10, order=3))
+
+        traj.add_phase(name='phase0', phase=phase)
+
+        #
+        # Set the time options
+        # Time has no targets in our ODE.
+        # We fix the initial time so that it is not a design variable in the optimization.
+        # The duration of the phase is allowed to be optimized, but is bounded on [0.5, 10].
+        #
+        phase.set_time_options(fix_initial=True, duration_bounds=(0.5, 10.0), units='s')
+
+        #
+        # Set the time options
+        # Initial values of positions and velocity are all fixed.
+        # The final value of position are fixed, but the final velocity is a free variable.
+        # The equations of motion are not functions of position, so 'x' and 'y' have no targets.
+        # The rate source points to the output in the ODE which provides the time derivative of the
+        # given state.
+        phase.add_state('x', fix_initial=True, fix_final=True, units='m',
+                        rate_source='xdot_comp.xdot')
+        phase.add_state('y', fix_initial=True, fix_final=True, units='m',
+                        rate_source='ydot_comp.ydot')
+        phase.add_state('v', fix_initial=True, fix_final=False, units='m/s',
+                        rate_source='vdot_comp.vdot', targets=['xdot_comp.v', 'ydot_comp.v'])
+
+        # Define theta as a control.
+        phase.add_control(name='theta', units='rad', lower=0, upper=np.pi,
+                          targets=['xdot_comp.theta', 'ydot_comp.theta', 'vdot_comp.theta'])
+
+        phase.add_parameter('foo',
+                            units='lbm',
+                            opt=False,
+                            targets=['xdot_comp.foo', 'ydot_comp.foo', 'vdot_comp.foo'])
+
+        # Minimize final time.
+        phase.add_objective('time', loc='final')
+
+        # Set the driver.
+        p.driver = om.ScipyOptimizeDriver()
+
+        # Allow OpenMDAO to automatically determine our sparsity pattern.
+        # Doing so can significant speed up the execution of Dymos.
+        p.driver.declare_coloring()
+
+        # Setup the problem
+        p.setup(check=True)
+
+        # Now that the OpenMDAO problem is setup, we can set the values of the states.
+        p.set_val('traj.phase0.t_initial', 0.0, units='s')
+        p.set_val('traj.phase0.t_duration', 2.0, units='s')
+        p.set_val('traj.phase0.states:x', phase.interp('x', [0, 10]), units='m')
+        p.set_val('traj.phase0.states:y', phase.interp('y', [10, 5]), units='m')
+        p.set_val('traj.phase0.states:v', phase.interp('v', [0, 5]), units='m/s')
+        p.set_val('traj.phase0.controls:theta', phase.interp('theta', [90, 90]), units='deg')
+        p.set_val('traj.phase0.parameters:foo', 5.0)
+
+        # Run the driver to solve the problem
+        p.run_driver()
+
+        expected = convert_units(5.0, 'lbm', 'kg')
+
+        self.assertEqual((1,), phase.parameter_options['foo']['shape'])
+        assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-5)
+        assert_near_equal(p.get_val('traj.phase0.parameter_vals:foo')[-1], 5.0, tolerance=1.0E-5)
+        assert_near_equal(p.get_val('traj.phase0.rhs_col.xdot_comp.foo'), expected*np.ones(10,), tolerance=1.0E-5)
+        assert_near_equal(p.get_val('traj.phase0.rhs_col.ydot_comp.foo'), expected*np.ones(10,), tolerance=1.0E-5)
+        assert_near_equal(p.get_val('traj.phase0.rhs_col.vdot_comp.foo'), expected*np.ones(10,), tolerance=1.0E-5)
+        assert_near_equal(p.get_val('traj.phase0.rhs_disc.xdot_comp.foo'), expected*np.ones(20,), tolerance=1.0E-5)
+        assert_near_equal(p.get_val('traj.phase0.rhs_disc.ydot_comp.foo'), expected*np.ones(20,), tolerance=1.0E-5)
+        assert_near_equal(p.get_val('traj.phase0.rhs_disc.vdot_comp.foo'), expected*np.ones(20,), tolerance=1.0E-5)
+
+    def test_invalid_params_different_target_units_introspection_failure(self):
+        #
+        # Define the OpenMDAO problem
+        #
+        p = om.Problem(model=om.Group())
+
+        #
+        # Define a Trajectory object
+        #
+        traj = dm.Trajectory()
+
+        p.model.add_subsystem('traj', subsys=traj)
+
+        #
+        # Define a Dymos Phase object with GaussLobatto Transcription
+        #
+        phase = dm.Phase(ode_class=_TestEOM,
+                         transcription=dm.GaussLobatto(num_segments=10, order=3),
+                         ode_init_kwargs={'foo_units': {'xdot_comp': 'kg',
+                                                        'ydot_comp': 'lbm',
+                                                        'vdot_comp': 'slug'}})
+
+        traj.add_phase(name='phase0', phase=phase)
+
+        #
+        # Set the time options
+        # Time has no targets in our ODE.
+        # We fix the initial time so that it is not a design variable in the optimization.
+        # The duration of the phase is allowed to be optimized, but is bounded on [0.5, 10].
+        #
+        phase.set_time_options(fix_initial=True, duration_bounds=(0.5, 10.0), units='s')
+
+        #
+        # Set the time options
+        # Initial values of positions and velocity are all fixed.
+        # The final value of position are fixed, but the final velocity is a free variable.
+        # The equations of motion are not functions of position, so 'x' and 'y' have no targets.
+        # The rate source points to the output in the ODE which provides the time derivative of the
+        # given state.
+        phase.add_state('x', fix_initial=True, fix_final=True, units='m',
+                        rate_source='xdot_comp.xdot')
+        phase.add_state('y', fix_initial=True, fix_final=True, units='m',
+                        rate_source='ydot_comp.ydot')
+        phase.add_state('v', fix_initial=True, fix_final=False, units='m/s',
+                        rate_source='vdot_comp.vdot', targets=['xdot_comp.v', 'ydot_comp.v'])
+
+        # Define theta as a control.
+        phase.add_control(name='theta', units='rad', lower=0, upper=np.pi,
+                          targets=['xdot_comp.theta', 'ydot_comp.theta', 'vdot_comp.theta'])
+
+        phase.add_parameter('foo',
+                            opt=False,
+                            targets=['xdot_comp.foo', 'ydot_comp.foo', 'vdot_comp.foo'])
+
+        # Minimize final time.
+        phase.add_objective('time', loc='final')
+
+        # Set the driver.
+        p.driver = om.ScipyOptimizeDriver()
+
+        # Allow OpenMDAO to automatically determine our sparsity pattern.
+        # Doing so can significant speed up the execution of Dymos.
+        p.driver.declare_coloring()
+
+        # Setup the problem
+        with self.assertRaises(RuntimeError) as e:
+            p.setup(check=True)
+
+        expected = ("Unable to automatically assign units based on targets.\n"
+                    "Targets have multiple units assigned:\n"
+                    "{'xdot_comp.foo': 'kg', 'ydot_comp.foo': 'lbm', 'vdot_comp.foo': 'slug'}.\n"
+                    "Either promote targets and use set_input_defaults to assign common\n"
+                    "units, or explicitly provide units to the variable.")
+
+        self.assertEqual(expected, str(e.exception))
+
+
+if __name__ == '__main__':  # pragma: no cover
+    unittest.main()
diff --git a/dymos/examples/cannonball/test/test_traj_param_static_and_dynamic.py b/dymos/examples/cannonball/test/test_traj_param_static_and_dynamic.py
new file mode 100644
index 000000000..d717576b8
--- /dev/null
+++ b/dymos/examples/cannonball/test/test_traj_param_static_and_dynamic.py
@@ -0,0 +1,220 @@
+import unittest
+
+import numpy as np
+
+import openmdao
+import openmdao.api as om
+
+from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
+
+from dymos.examples.cannonball.cannonball_ode import rho_interp
+
+om_dev_version = openmdao.__version__.endswith('dev')
+om_version = tuple(int(s) for s in openmdao.__version__.split('-')[0].split('.'))
+
+
+GRAVITY = 9.80665  # m/s**2
+
+
+class CannonballODEVectorCD(om.ExplicitComponent):
+    """
+    Cannonball ODE assuming flat earth and accounting for air resistance
+    """
+
+    def initialize(self):
+        self.options.declare('num_nodes', types=int)
+        self.options.declare('static_gravity', types=bool)
+
+    def setup(self):
+        nn = self.options['num_nodes']
+
+        # static parameters
+        self.add_input('m', units='kg')
+        self.add_input('S', units='m**2')
+
+        if self.options['static_gravity']:
+            self.add_input('g', units='m/s**2', val=GRAVITY)
+        else:
+            self.add_input('g', units='m/s**2', val=GRAVITY * np.ones(nn,))
+
+        # This will be used as both a control and a parameter
+        self.add_input('CD', 0.5, shape=nn)
+
+        # time varying inputs
+        self.add_input('h', units='m', shape=nn)
+        self.add_input('v', units='m/s', shape=nn)
+        self.add_input('gam', units='rad', shape=nn)
+
+        # state rates
+        self.add_output('v_dot', shape=nn, units='m/s**2')
+        self.add_output('gam_dot', shape=nn, units='rad/s')
+        self.add_output('h_dot', shape=nn, units='m/s')
+        self.add_output('r_dot', shape=nn, units='m/s')
+        self.add_output('ke', shape=nn, units='J')
+
+        # Ask OpenMDAO to compute the partial derivatives using complex-step
+        # with a partial coloring algorithm for improved performance
+        self.declare_partials('*', '*', method='cs')
+        self.declare_coloring(wrt='*', method='cs')
+
+    def compute(self, inputs, outputs):
+
+        gam = inputs['gam']
+        v = inputs['v']
+        h = inputs['h']
+        m = inputs['m']
+        S = inputs['S']
+        CD = inputs['CD']
+
+        # handle complex-step gracefully from the interpolant
+        if np.iscomplexobj(h):
+            rho = rho_interp(inputs['h'])
+        else:
+            rho = rho_interp(inputs['h']).real
+
+        q = 0.5*rho*inputs['v']**2
+        qS = q * S
+        D = qS * CD
+        cgam = np.cos(gam)
+        sgam = np.sin(gam)
+        outputs['v_dot'] = - D/m-GRAVITY*sgam
+        outputs['gam_dot'] = -(GRAVITY/v)*cgam
+        outputs['h_dot'] = v*sgam
+        outputs['r_dot'] = v*cgam
+        outputs['ke'] = 0.5*m*v**2
+
+
+@use_tempdirs
+class TestTrajParamStaticAndDynamic(unittest.TestCase):
+
+    @require_pyoptsparse(optimizer='SLSQP')
+    def test_traj_param_static_and_dynamic(self):
+        """ Test that the final value of a control in one phase can be connected as the value
+        of a parameter in a subsequent phase. """
+        import openmdao.api as om
+        from openmdao.utils.assert_utils import assert_near_equal
+
+        import dymos as dm
+        from dymos.examples.cannonball.size_comp import CannonballSizeComp
+
+        p = om.Problem(model=om.Group())
+
+        p.driver = om.pyOptSparseDriver()
+        p.driver.options['optimizer'] = 'SLSQP'
+        p.driver.declare_coloring()
+
+        external_params = p.model.add_subsystem('external_params', om.IndepVarComp())
+
+        external_params.add_output('radius', val=0.10, units='m')
+        external_params.add_output('dens', val=7.87, units='g/cm**3')
+
+        external_params.add_design_var('radius', lower=0.01, upper=0.10, ref0=0.01, ref=0.10)
+
+        p.model.add_subsystem('size_comp', CannonballSizeComp())
+
+        traj = p.model.add_subsystem('traj', dm.Trajectory())
+
+        transcription = dm.Radau(num_segments=5, order=3, compressed=True)
+        ascent = dm.Phase(ode_class=CannonballODEVectorCD, transcription=transcription,
+                          ode_init_kwargs={'static_gravity': True})
+
+        ascent = traj.add_phase('ascent', ascent)
+
+        # All initial states except flight path angle are fixed
+        # Final flight path angle is fixed (we will set it to zero so that the phase ends at apogee)
+
+        ascent.set_time_options(fix_initial=True, duration_bounds=(1, 100), duration_ref=100, units='s')
+        ascent.add_state('r', fix_initial=True, fix_final=False, rate_source='r_dot', units='m')
+        ascent.add_state('h', fix_initial=True, fix_final=False, units='m', rate_source='h_dot')
+        ascent.add_state('gam', fix_initial=False, fix_final=True, units='rad', rate_source='gam_dot')
+        ascent.add_state('v', fix_initial=False, fix_final=False, units='m/s', rate_source='v_dot')
+
+        ascent.add_parameter('S', targets=['S'], units='m**2', static_target=True)
+        ascent.add_parameter('mass', targets=['m'], units='kg', static_target=True)
+
+        ascent.add_control('CD', targets=['CD'], opt=False, val=0.05)
+
+        # Limit the muzzle energy
+        ascent.add_boundary_constraint('ke', loc='initial',
+                                       upper=400000, lower=0, ref=100000)
+
+        # Second Phase (descent)
+        transcription = dm.GaussLobatto(num_segments=5, order=3, compressed=True)
+        descent = dm.Phase(ode_class=CannonballODEVectorCD, transcription=transcription,
+                           ode_init_kwargs={'static_gravity': False})
+
+        traj.add_phase('descent', descent)
+
+        # All initial states and time are free (they will be linked to the final states of ascent.
+        # Final altitude is fixed (we will set it to zero so that the phase ends at ground impact)
+        descent.set_time_options(initial_bounds=(.5, 100), duration_bounds=(.5, 100),
+                                 duration_ref=100, units='s')
+        descent.add_state('r', units='m', rate_source='r_dot')
+        descent.add_state('h', units='m', rate_source='h_dot', fix_initial=False, fix_final=True)
+        descent.add_state('gam', units='rad', rate_source='gam_dot', fix_initial=False, fix_final=False)
+        descent.add_state('v', units='m/s', rate_source='v_dot', fix_initial=False, fix_final=False)
+
+        descent.add_parameter('S', targets=['S'], units='m**2', static_target=True)
+        descent.add_parameter('mass', targets=['m'], units='kg', static_target=True)
+        descent.add_parameter('CD', targets=['CD'], val=0.01)
+
+        descent.add_objective('r', loc='final', scaler=-1.0)
+
+        # Add externally-provided design parameters to the trajectory.
+        # In this case, we connect 'm' to pre-existing input parameters named 'mass' in each phase.
+        traj.add_parameter('m', units='kg', val=1.0,
+                           targets={'ascent': 'mass', 'descent': 'mass'}, static_target=True)
+
+        # In this case, by omitting targets, we're connecting these parameters to parameters
+        # with the same name in each phase.
+        traj.add_parameter('S', units='m**2', val=0.005, static_target=True)
+
+        # Link Phases (link time and all state variables)
+        traj.link_phases(phases=['ascent', 'descent'], vars=['*'])
+
+        # Issue Connections
+        p.model.connect('external_params.radius', 'size_comp.radius')
+        p.model.connect('external_params.dens', 'size_comp.dens')
+
+        p.model.connect('size_comp.mass', 'traj.parameters:m')
+        p.model.connect('size_comp.S', 'traj.parameters:S')
+
+        traj.connect('ascent.timeseries.CD', 'descent.parameters:CD', src_indices=[-1], flat_src_indices=True)
+
+        # A linear solver at the top level can improve performance.
+        p.model.linear_solver = om.DirectSolver()
+
+        # Finish Problem Setup
+        p.setup()
+
+        # Set Initial Guesses
+        p.set_val('external_params.radius', 0.05, units='m')
+        p.set_val('external_params.dens', 7.87, units='g/cm**3')
+
+        p.set_val('traj.ascent.controls:CD', 0.5)
+
+        p.set_val('traj.ascent.t_initial', 0.0)
+        p.set_val('traj.ascent.t_duration', 10.0)
+
+        p.set_val('traj.ascent.states:r', ascent.interp('r', [0, 100]))
+        p.set_val('traj.ascent.states:h', ascent.interp('h', [0, 100]))
+        p.set_val('traj.ascent.states:v', ascent.interp('v', [200, 150]))
+        p.set_val('traj.ascent.states:gam', ascent.interp('gam', [25, 0]), units='deg')
+
+        p.set_val('traj.descent.t_initial', 10.0)
+        p.set_val('traj.descent.t_duration', 10.0)
+
+        p.set_val('traj.descent.states:r', descent.interp('r', [100, 200]))
+        p.set_val('traj.descent.states:h', descent.interp('h', [100, 0]))
+        p.set_val('traj.descent.states:v', descent.interp('v', [150, 200]))
+        p.set_val('traj.descent.states:gam', descent.interp('gam', [0, -45]), units='deg')
+
+        dm.run_problem(p, simulate=True, make_plots=True)
+
+        assert_near_equal(p.get_val('traj.descent.states:r')[-1], 3183.25, tolerance=1.0E-2)
+        assert_near_equal(p.get_val('traj.ascent.timeseries.CD')[-1],
+                          p.get_val('traj.descent.parameter_vals:CD')[0])
+
+
+if __name__ == '__main__':  # pragma: no cover
+    unittest.main()
diff --git a/dymos/utils/introspection.py b/dymos/utils/introspection.py
index 76486fdff..1dd35e284 100644
--- a/dymos/utils/introspection.py
+++ b/dymos/utils/introspection.py
@@ -411,29 +411,40 @@ def configure_parameters_introspection(parameter_options, ode):
         if options['units'] is _unspecified:
             options['units'] = _get_common_metadata(targets, metadata_key='units')
 
+        # Check that all targets have the same shape.
+        static_shapes = {}
+        dynamic_shapes = {}
+        # First find the shapes of the static targets
+        for tgt, meta in targets.items():
+            if tgt in options['static_targets']:
+                static_shapes[tgt] = meta['shape']
+            else:
+                if len(meta['shape']) == 1:
+                    dynamic_shapes[tgt] = (1,)
+                else:
+                    dynamic_shapes[tgt] = meta['shape'][1:]
+        all_shapes = {**dynamic_shapes, **static_shapes}
+        # Check that they're unique
+        if len(set(all_shapes.values())) > 1:
+            raise RuntimeError(f'Invalid targets for parameter `{name}`.\n'
+                               f'Targets have multiple shapes.\n'
+                               f'{all_shapes}')
+        elif len(set(all_shapes.values())) == 1:
+            introspected_shape = next(iter(set(all_shapes.values())))
+        else:
+            introspected_shape = None
+
         if options['shape'] in {_unspecified, None}:
             if isinstance(options['val'], Number):
-                static_shapes = {}
-                dynamic_shapes = {}
-                # First find the shapes of the static targets
-                for tgt, meta in targets.items():
-                    if tgt in options['static_targets']:
-                        static_shapes[tgt] = meta['shape']
-                    else:
-                        if len(meta['shape']) == 1:
-                            dynamic_shapes[tgt] = (1,)
-                        else:
-                            dynamic_shapes[tgt] = meta['shape'][1:]
-                all_shapes = {**dynamic_shapes, **static_shapes}
-                # Check that they're unique
-                if len(set(all_shapes.values())) != 1:
-                    raise RuntimeError(f'Unable to obtain shape of parameter {name} via introspection.\n'
-                                       f'Targets have multiple shapes.\n'
-                                       f'{all_shapes}')
-                else:
-                    options['shape'] = next(iter(set(all_shapes.values())))
+                options['shape'] = introspected_shape
             else:
                 options['shape'] = np.asarray(options['val']).shape
+        else:
+            if introspected_shape is not None and options['shape'] != introspected_shape:
+                raise RuntimeError(f'Shape provided to parameter `{name}` differs from its targets.\n'
+                                   f'Given shape: {options["shape"]}\n'
+                                   f'Target shapes:\n'
+                                   f'{all_shapes}')
 
         options['val'], options['shape'] = ensure_compatible(name, options['val'], options['shape'])
 
@@ -1187,14 +1198,14 @@ def _get_common_metadata(targets, metadata_key):
     if len(meta_set) == 1:
         return next(iter(meta_set))
     elif len(meta_set) == 0:
-        raise ValueError(f'Unable to automatically assign {metadata_key} based on targets. \n'
-                         f'No targets were found.')
+        raise RuntimeError(f'Unable to automatically assign {metadata_key} based on targets. \n'
+                           f'No targets were found.')
     else:
         err_dict = {tgt: meta[metadata_key] for tgt, meta in targets.items()}
-        raise ValueError(f'Unable to automatically assign {metadata_key} based on targets. \n'
-                         f'Targets have multiple {metadata_key} assigned: {err_dict}. \n'
-                         f'Either promote targets and use set_input_defaults to assign common '
-                         f'{metadata_key}, or explicitly provide {metadata_key} to the variable.')
+        raise RuntimeError(f'Unable to automatically assign {metadata_key} based on targets.\n'
+                           f'Targets have multiple {metadata_key} assigned:\n{err_dict}.\n'
+                           f'Either promote targets and use set_input_defaults to assign common\n'
+                           f'{metadata_key}, or explicitly provide {metadata_key} to the variable.')
 
 
 def get_source_metadata(ode, src, user_units=_unspecified, user_shape=_unspecified):