Repository: OpenMDAO/dymos

Issue:
Polynomial control rates are not correctly evaluated in explicit shooting.
### Description

When a polynomial control has `rate_targets` in the ODE, the resulting simulation obtains an incorrect value for that rate because it relies on the ration of phase duration to phase "tau" space. (`0.5 * t_duration`).

`t_duration` doesn't seem to be correctly connected to the integration component, resulting in incorrect evaluations of the control rate and therefore incorrect propagation.

### Example

The following code tests control rate targets in both the optimization and simulation for multiple transcriptions and both polynomial and "full" controls.

In this case, make_plots is True to demonstrate the issue with the polynomial control cases.  They will also fail the assertion that the simulation and solution be nearly the same.

```python3
import unittest

import numpy as np

try:
    import matplotlib
except ImportError:
    matplotlib = None

import openmdao.api as om
from openmdao.utils.testing_utils import use_tempdirs
from dymos.utils.misc import om_version
from dymos.utils.testing_utils import assert_timeseries_near_equal
from dymos.utils.misc import _unspecified

class BrachistochroneRateTargetODE(om.ExplicitComponent):
    #
    # The following dictionaries provide a way of 'tagging' the Brachistochrone ODE with
    # information about states and parameters that can be accessed from Phase.
    #
    # In this case these are class attributes, but the choice of whether or not to tie
    # this information to the ODE itself (and how to do so) is entirely up to the user.
    #
    # In a dynamic ODE model these might be instance attributes whose values vary depending on
    # the arguments to the instantiation.
    #
    states = {'x': {'rate_source': 'xdot',
                    'units': 'm'},
              'y': {'rate_source': 'ydot',
                    'units': 'm'},
              'v': {'rate_source': 'vdot',
                    'units': 'm/s'}}

    parameters = {'theta': {'units': 'rad'},
                  'g': {'units': 'm/s**2'}}

    def initialize(self):
        self.options.declare('num_nodes', types=int)
        self.options.declare('static_gravity', types=(bool,), default=False,
                             desc='If True, treat gravity as a static (scalar) input, rather than '
                                  'having different values at each node.')
        self.options.declare('control_name', values=('theta', 'int_theta_rate'), default='int_theta_rate')

    def setup(self):
        nn = self.options['num_nodes']
        control_name = self.options['control_name']
        g_default_val = 9.80665 if self.options['static_gravity'] else 9.80665 * np.ones(nn)

        # Inputs
        self.add_input('v', val=np.zeros(nn), desc='velocity', units='m/s')

        self.add_input('g', val=g_default_val, desc='grav. acceleration', units='m/s/s')

        # Note, kind of strange for this to be named theta_rate, but this is just a demonstration that
        # we can connect to a control rate source in dymos.
        self.add_input(control_name, val=np.ones(nn), desc='angle of wire', units='rad')

        self.add_output('xdot', val=np.zeros(nn), desc='velocity component in x', units='m/s')

        self.add_output('ydot', val=np.zeros(nn), desc='velocity component in y', units='m/s')

        self.add_output('vdot', val=np.zeros(nn), desc='acceleration magnitude', units='m/s**2')

        self.add_output('check', val=np.zeros(nn),
                        desc='check solution: v/sin(theta) = constant',
                        units='m/s')

        # Setup partials
        arange = np.arange(self.options['num_nodes'])
        self.declare_partials(of='vdot', wrt=control_name, rows=arange, cols=arange)

        self.declare_partials(of='xdot', wrt='v', rows=arange, cols=arange)
        self.declare_partials(of='xdot', wrt=control_name, rows=arange, cols=arange)

        self.declare_partials(of='ydot', wrt='v', rows=arange, cols=arange)
        self.declare_partials(of='ydot', wrt=control_name, rows=arange, cols=arange)

        self.declare_partials(of='check', wrt='v', rows=arange, cols=arange)
        self.declare_partials(of='check', wrt=control_name, rows=arange, cols=arange)

        if self.options['static_gravity']:
            c = np.zeros(self.options['num_nodes'])
            self.declare_partials(of='vdot', wrt='g', rows=arange, cols=c)
        else:
            self.declare_partials(of='vdot', wrt='g', rows=arange, cols=arange)

    def compute(self, inputs, outputs):
        u_name = self.options['control_name']
        theta = inputs[u_name]
        cos_theta = np.cos(theta)
        sin_theta = np.sin(theta)
        g = inputs['g']
        v = inputs['v']

        outputs['vdot'] = g * cos_theta
        outputs['xdot'] = v * sin_theta
        outputs['ydot'] = -v * cos_theta
        outputs['check'] = v / sin_theta

    def compute_partials(self, inputs, partials):
        u_name = self.options['control_name']
        theta = inputs[u_name]
        cos_theta = np.cos(theta)
        sin_theta = np.sin(theta)
        g = inputs['g']
        v = inputs['v']

        partials['vdot', 'g'] = cos_theta
        partials['vdot', u_name] = -g * sin_theta

        partials['xdot', 'v'] = sin_theta
        partials['xdot', u_name] = v * cos_theta

        partials['ydot', 'v'] = -cos_theta
        partials['ydot', u_name] = v * sin_theta

        partials['check', 'v'] = 1 / sin_theta
        partials['check', u_name] = -v * cos_theta / sin_theta ** 2


#@use_tempdirs
class TestBrachistochroneControlRateTargets(unittest.TestCase):

    def test_brachistochrone_control_rate_targets(self):

        import openmdao.api as om
        from openmdao.utils.assert_utils import assert_near_equal
        import dymos as dm

        transcriptions = {'gauss-lobatto': dm.GaussLobatto(num_segments=10),}
                         # 'radau': dm.Radau(num_segments=10),
                         # 'birkhoff': dm.Birkhoff(num_nodes=30)}

        for tx_name, tx in transcriptions.items():
            for control_target_method in ('implicit', 'explicit'):
                for control_type in ('full', 'polynomial'):

                    with self.subTest(f'{tx_name=} {control_target_method=} {control_type=}'):

                        p = om.Problem(model=om.Group())
                        p.driver = om.pyOptSparseDriver(optimizer='SLSQP', print_results=False)
                        p.driver.declare_coloring()

                        traj = dm.Trajectory()

                        ode_kwargs = {'control_name': 'int_theta_rate' if control_target_method == 'implicit' else 'theta'}

                        phase = dm.Phase(ode_class=BrachistochroneRateTargetODE,
                                        ode_init_kwargs=ode_kwargs,
                                        transcription=tx)

                        traj.add_phase('phase0', phase)

                        p.model.add_subsystem('traj0', traj)

                        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))

                        phase.add_state('x', rate_source='xdot',
                                        units='m',
                                        fix_initial=True, fix_final=True, solve_segments=False)

                        phase.add_state('y', rate_source='ydot',
                                        units='m',
                                        fix_initial=True, fix_final=True, solve_segments=False)

                        phase.add_state('v', rate_source='vdot',
                                        units='m/s',
                                        fix_initial=True, fix_final=False, solve_segments=False)

                        # our control is the time integral of theta.  We pass its rate to the ODE and it should implicitly find it since
                        # it is named `int_theta_rate` in the ODE.
                        phase.add_control('int_theta', lower=0.0, upper=None, fix_initial=True,
                                        rate_targets=_unspecified if control_target_method == 'implicit' else ['theta'],
                                        control_type=control_type, order=7)

                        phase.add_parameter('g', units='m/s**2', opt=False, val=9.80665)

                        # Minimize time at the end of the phase
                        phase.add_objective('time', loc='final', scaler=10)

                        p.model.linear_solver = om.DirectSolver()

                        p.setup()

                        phase.set_time_val(initial=0.0, duration=2.0)

                        phase.set_state_val('x', [0, 10])
                        phase.set_state_val('y', [10, 5])
                        phase.set_state_val('v', [0, 9.9])
                        phase.set_control_val('int_theta', [0, 100], units='deg*s')

                        p.run_model()

                        # Solve for the optimal trajectory
                        dm.run_problem(p, simulate=True, make_plots=True, simulate_kwargs={'reports': True})

                        sol_db = 'dymos_solution.db'
                        sim_db = 'dymos_simulation.db'
                        if om_version()[0] > (3, 34, 2):
                            sol_db = p.get_outputs_dir() / sol_db
                            sim_db = traj.sim_prob.get_outputs_dir() / sim_db

                        sol_case = om.CaseReader(sol_db).get_case('final')
                        sim_case = om.CaseReader(sim_db).get_case('final')

                        t_sol = sol_case.get_val('traj0.phase0.timeseries.time')
                        x_sol = sol_case.get_val('traj0.phase0.timeseries.x')
                        y_sol = sol_case.get_val('traj0.phase0.timeseries.y')
                        v_sol = sol_case.get_val('traj0.phase0.timeseries.v')
                        theta_sol = sol_case.get_val('traj0.phase0.timeseries.int_theta_rate')

                        t_sim = sim_case.get_val('traj0.phase0.timeseries.time')
                        x_sim = sim_case.get_val('traj0.phase0.timeseries.x')
                        y_sim = sim_case.get_val('traj0.phase0.timeseries.y')
                        v_sim = sim_case.get_val('traj0.phase0.timeseries.v')
                        theta_sim = sim_case.get_val('traj0.phase0.timeseries.int_theta_rate')

                        assert_timeseries_near_equal(t_sol, x_sol, t_sim, x_sim, rel_tolerance=5.0E-3, abs_tolerance=0.01)
                        assert_timeseries_near_equal(t_sol, y_sol, t_sim, y_sim, rel_tolerance=5.0E-3, abs_tolerance=0.01)
                        assert_timeseries_near_equal(t_sol, v_sol, t_sim, v_sim, rel_tolerance=5.0E-3, abs_tolerance=0.01)
                        assert_timeseries_near_equal(t_sol, theta_sol, t_sim, theta_sim, rel_tolerance=5.0E-3, abs_tolerance=0.01)
```

### Dymos Version

1.12.1-dev

### Relevant environment information

_No response_

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

coverage: 93.589% (-0.3%) from 93.88%
when pulling **010ec79bd378607151dd92a037724dac56e410e2 on robfalck:control_rate_bug**
into **42fea76c302d94fd31c58f67799a95b6a5debf35 on OpenMDAO:master**.


Patch:
diff --git a/dymos/transcriptions/explicit_shooting/explicit_shooting.py b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
index a194f9e5..d60f2912 100644
--- a/dymos/transcriptions/explicit_shooting/explicit_shooting.py
+++ b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
@@ -355,7 +355,8 @@ class ExplicitShooting(TranscriptionBase):
 
             phase.add_subsystem('control_group',
                                 subsys=control_group,
-                                promotes=['dt_dstau', '*controls:*', '*control_values:*', '*control_rates:*'])
+                                promotes=[('t_duration', 't_duration_val'), 'dt_dstau',
+                                          '*controls:*', '*control_values:*', '*control_rates:*'])
 
             control_prefix = 'controls:' if phase.timeseries_options['use_prefix'] else ''
             control_rate_prefix = 'control_rates:' if phase.timeseries_options['use_prefix'] else ''
diff --git a/dymos/visualization/timeseries/bokeh_timeseries_report.py b/dymos/visualization/timeseries/bokeh_timeseries_report.py
index 99b501ca..ca7ef5bb 100644
--- a/dymos/visualization/timeseries/bokeh_timeseries_report.py
+++ b/dymos/visualization/timeseries/bokeh_timeseries_report.py
@@ -286,7 +286,9 @@ def _load_data_sources(traj_and_phase_meta=None, solution_record_file=None, simu
     if sim_outputs is not None and sol_outputs is not None:
         if not set(sim_outputs.keys()).issubset(sol_outputs.keys()):
             om.issue_warning('Simulation file does not contain the same outputs as the solution '
-                             'file. Skipping plotting of simulation timeseries data.')
+                             'file. Skipping plotting of simulation timeseries data.\nThe following '
+                             'outputs are in the simulation results but not in the solution results:\n'
+                             f'{set(sim_outputs.keys()) - set(sol_outputs.keys())}')
             sim_case = None
 
     for traj_path, traj_data in traj_and_phase_meta.items():


Test patch:
diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py b/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py
index 13578dea..caa52d4d 100644
--- a/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py
+++ b/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py
@@ -9,6 +9,9 @@ except ImportError:
 
 import openmdao.api as om
 from openmdao.utils.testing_utils import use_tempdirs
+from dymos.utils.misc import om_version
+from dymos.utils.testing_utils import assert_timeseries_near_equal
+from dymos.utils.misc import _unspecified
 
 
 class BrachistochroneRateTargetODE(om.ExplicitComponent):
@@ -37,9 +40,11 @@ class BrachistochroneRateTargetODE(om.ExplicitComponent):
         self.options.declare('static_gravity', types=(bool,), default=False,
                              desc='If True, treat gravity as a static (scalar) input, rather than '
                                   'having different values at each node.')
+        self.options.declare('control_name', values=('theta', 'int_theta_rate'), default='int_theta_rate')
 
     def setup(self):
         nn = self.options['num_nodes']
+        control_name = self.options['control_name']
         g_default_val = 9.80665 if self.options['static_gravity'] else 9.80665 * np.ones(nn)
 
         # Inputs
@@ -49,7 +54,7 @@ class BrachistochroneRateTargetODE(om.ExplicitComponent):
 
         # Note, kind of strange for this to be named theta_rate, but this is just a demonstration that
         # we can connect to a control rate source in dymos.
-        self.add_input('theta_rate', val=np.ones(nn), desc='angle of wire', units='rad')
+        self.add_input(control_name, val=np.ones(nn), desc='angle of wire', units='rad')
 
         self.add_output('xdot', val=np.zeros(nn), desc='velocity component in x', units='m/s')
 
@@ -63,16 +68,16 @@ class BrachistochroneRateTargetODE(om.ExplicitComponent):
 
         # Setup partials
         arange = np.arange(self.options['num_nodes'])
-        self.declare_partials(of='vdot', wrt='theta_rate', rows=arange, cols=arange)
+        self.declare_partials(of='vdot', wrt=control_name, rows=arange, cols=arange)
 
         self.declare_partials(of='xdot', wrt='v', rows=arange, cols=arange)
-        self.declare_partials(of='xdot', wrt='theta_rate', rows=arange, cols=arange)
+        self.declare_partials(of='xdot', wrt=control_name, rows=arange, cols=arange)
 
         self.declare_partials(of='ydot', wrt='v', rows=arange, cols=arange)
-        self.declare_partials(of='ydot', wrt='theta_rate', rows=arange, cols=arange)
+        self.declare_partials(of='ydot', wrt=control_name, rows=arange, cols=arange)
 
         self.declare_partials(of='check', wrt='v', rows=arange, cols=arange)
-        self.declare_partials(of='check', wrt='theta_rate', rows=arange, cols=arange)
+        self.declare_partials(of='check', wrt=control_name, rows=arange, cols=arange)
 
         if self.options['static_gravity']:
             c = np.zeros(self.options['num_nodes'])
@@ -81,7 +86,8 @@ class BrachistochroneRateTargetODE(om.ExplicitComponent):
             self.declare_partials(of='vdot', wrt='g', rows=arange, cols=arange)
 
     def compute(self, inputs, outputs):
-        theta = inputs['theta_rate']
+        u_name = self.options['control_name']
+        theta = inputs[u_name]
         cos_theta = np.cos(theta)
         sin_theta = np.sin(theta)
         g = inputs['g']
@@ -93,850 +99,127 @@ class BrachistochroneRateTargetODE(om.ExplicitComponent):
         outputs['check'] = v / sin_theta
 
     def compute_partials(self, inputs, partials):
-        theta = inputs['theta_rate']
+        u_name = self.options['control_name']
+        theta = inputs[u_name]
         cos_theta = np.cos(theta)
         sin_theta = np.sin(theta)
         g = inputs['g']
         v = inputs['v']
 
         partials['vdot', 'g'] = cos_theta
-        partials['vdot', 'theta_rate'] = -g * sin_theta
+        partials['vdot', u_name] = -g * sin_theta
 
         partials['xdot', 'v'] = sin_theta
-        partials['xdot', 'theta_rate'] = v * cos_theta
+        partials['xdot', u_name] = v * cos_theta
 
         partials['ydot', 'v'] = -cos_theta
-        partials['ydot', 'theta_rate'] = v * sin_theta
+        partials['ydot', u_name] = v * sin_theta
 
         partials['check', 'v'] = 1 / sin_theta
-        partials['check', 'theta_rate'] = -v * cos_theta / sin_theta ** 2
+        partials['check', u_name] = -v * cos_theta / sin_theta ** 2
 
 
 @use_tempdirs
 class TestBrachistochroneControlRateTargets(unittest.TestCase):
 
-    @unittest.skipIf(matplotlib is None, "This test requires matplotlib")
-    def test_brachistochrone_control_rate_targets_gauss_lobatto(self):
-        import matplotlib.pyplot as plt
-        plt.switch_backend('Agg')
-        import openmdao.api as om
-        from openmdao.utils.assert_utils import assert_near_equal
-        import dymos as dm
-
-        p = om.Problem(model=om.Group())
-        p.driver = om.ScipyOptimizeDriver()
-        p.driver.declare_coloring()
-
-        phase = dm.Phase(ode_class=BrachistochroneRateTargetODE,
-                         transcription=dm.GaussLobatto(num_segments=10))
-
-        p.model.add_subsystem('phase0', phase)
-
-        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
-
-        phase.add_state('x', rate_source='xdot',
-                        units='m',
-                        fix_initial=True, fix_final=True, solve_segments=False)
-
-        phase.add_state('y', rate_source='ydot',
-                        units='m',
-                        fix_initial=True, fix_final=True, solve_segments=False)
-
-        phase.add_state('v', rate_source='vdot',
-                        units='m/s',
-                        fix_initial=True, fix_final=False, solve_segments=False)
-
-        phase.add_control('theta', lower=0.01, upper=179.9, fix_initial=True)
-
-        phase.add_parameter('g', units='m/s**2', opt=False, val=9.80665)
-
-        # Minimize time at the end of the phase
-        phase.add_objective('time', loc='final', scaler=10)
-
-        p.model.linear_solver = om.DirectSolver()
-
-        p.setup()
-
-        phase.set_time_val(initial=0.0, duration=2.0)
-
-        phase.set_state_val('x', [0, 10])
-        phase.set_state_val('y', [10, 5])
-        phase.set_state_val('v', [0, 9.9])
-        phase.set_control_val('theta', [0, 100])
-
-        # Solve for the optimal trajectory
-        p.run_driver()
-
-        # Test the results
-        assert_near_equal(p.get_val('phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-3)
-
-        # Generate the explicitly simulated trajectory
-        exp_out = phase.simulate()
-
-        fig, ax = plt.subplots()
-        fig.suptitle('Brachistochrone Solution')
-
-        x_imp = p.get_val('phase0.timeseries.x')
-        y_imp = p.get_val('phase0.timeseries.y')
-
-        x_exp = exp_out.get_val('phase0.timeseries.x')
-        y_exp = exp_out.get_val('phase0.timeseries.y')
-
-        ax.plot(x_imp, y_imp, 'ro', label='solution')
-        ax.plot(x_exp, y_exp, 'b-', label='simulated')
-
-        ax.set_xlabel('x (m)')
-        ax.set_ylabel('y (m)')
-        ax.grid(True)
-        ax.legend(loc='upper right')
-
-        fig, ax = plt.subplots()
-
-        t_imp = p.get_val('phase0.timeseries.time')
-        theta_imp = p.get_val('phase0.timeseries.theta_rate')
-        t_exp = exp_out.get_val('phase0.timeseries.time')
-        theta_exp = exp_out.get_val('phase0.timeseries.theta_rate')
-
-        io_meta = p.model.phase0.timeseries.timeseries_comp.get_io_metadata(
-            iotypes=('input', 'output'), get_remote=True)
-        self.assertEqual(io_meta['theta']['units'], 'rad*s')
-
-        ax.plot(t_imp, theta_imp, 'ro', label='solution')
-        ax.plot(t_exp, theta_exp, 'b-', label='simulated')
-
-        ax.set_xlabel('time (s)')
-        ax.set_ylabel(r'$\theta$ (deg)')
-        ax.grid(True)
-        ax.legend(loc='upper right')
-
-        plt.show()
-
-    @unittest.skipIf(matplotlib is None, "This test requires matplotlib")
-    def test_brachistochrone_control_rate_targets_radau(self):
-        import matplotlib.pyplot as plt
-        plt.switch_backend('Agg')
-        import matplotlib.pyplot as plt
-        import openmdao.api as om
-        from openmdao.utils.assert_utils import assert_near_equal
-        import dymos as dm
-
-        p = om.Problem(model=om.Group())
-        p.driver = om.ScipyOptimizeDriver()
-        p.driver.declare_coloring()
-
-        phase = dm.Phase(ode_class=BrachistochroneRateTargetODE,
-                         transcription=dm.Radau(num_segments=10))
-
-        p.model.add_subsystem('phase0', phase)
-
-        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
-
-        phase.add_state('x', rate_source='xdot',
-                        units='m',
-                        fix_initial=True, fix_final=True, solve_segments=False)
-
-        phase.add_state('y', rate_source='ydot',
-                        units='m',
-                        fix_initial=True, fix_final=True, solve_segments=False)
-
-        phase.add_state('v', rate_source='vdot',
-                        units='m/s',
-                        fix_initial=True, fix_final=False, solve_segments=False)
-
-        phase.add_control('theta', units='deg*s', lower=0.01, upper=179.9, fix_initial=True)
-
-        phase.add_parameter('g', units='m/s**2', opt=False, val=9.80665)
-
-        # Minimize time at the end of the phase
-        phase.add_objective('time', loc='final', scaler=10)
-
-        p.model.linear_solver = om.DirectSolver()
-
-        p.setup()
-
-        phase.set_time_val(initial=0.0, duration=2.0)
-
-        phase.set_state_val('x', [0, 10])
-        phase.set_state_val('y', [10, 5])
-        phase.set_state_val('v', [0, 9.9])
-        phase.set_control_val('theta', [0, 100])
-
-        # Solve for the optimal trajectory
-        p.run_driver()
-
-        # Test the results
-        assert_near_equal(p.get_val('phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-3)
-
-        # Generate the explicitly simulated trajectory
-        exp_out = phase.simulate()
-
-        fig, ax = plt.subplots()
-        fig.suptitle('Brachistochrone Solution')
-
-        x_imp = p.get_val('phase0.timeseries.x')
-        y_imp = p.get_val('phase0.timeseries.y')
-
-        x_exp = exp_out.get_val('phase0.timeseries.x')
-        y_exp = exp_out.get_val('phase0.timeseries.y')
-
-        ax.plot(x_imp, y_imp, 'ro', label='solution')
-        ax.plot(x_exp, y_exp, 'b-', label='simulated')
-
-        ax.set_xlabel('x (m)')
-        ax.set_ylabel('y (m)')
-        ax.grid(True)
-        ax.legend(loc='upper right')
-
-        fig, ax = plt.subplots()
-
-        t_imp = p.get_val('phase0.timeseries.time')
-        theta_imp = p.get_val('phase0.timeseries.theta_rate')
-        t_exp = exp_out.get_val('phase0.timeseries.time')
-        theta_exp = exp_out.get_val('phase0.timeseries.theta_rate')
-
-        ax.plot(t_imp, theta_imp, 'ro', label='solution')
-        ax.plot(t_exp, theta_exp, 'b-', label='simulated')
-
-        ax.set_xlabel('time (s)')
-        ax.set_ylabel(r'$\theta$ (deg)')
-        ax.grid(True)
-        ax.legend(loc='upper right')
-
-        plt.show()
-
-
-@use_tempdirs
-class TestBrachistochroneExplicitControlRateTargets(unittest.TestCase):
-
-    @unittest.skipIf(matplotlib is None, "This test requires matplotlib")
-    def test_brachistochrone_control_rate_targets_gauss_lobatto(self):
-        import matplotlib.pyplot as plt
-        plt.switch_backend('Agg')
-        import openmdao.api as om
-        from openmdao.utils.assert_utils import assert_near_equal
-        import dymos as dm
-
-        p = om.Problem(model=om.Group())
-        p.driver = om.ScipyOptimizeDriver()
-        p.driver.declare_coloring()
-
-        phase = dm.Phase(ode_class=BrachistochroneRateTargetODE,
-                         transcription=dm.GaussLobatto(num_segments=10))
-
-        p.model.add_subsystem('phase0', phase)
-
-        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
-
-        phase.add_state('x', rate_source='xdot',
-                        units='m',
-                        fix_initial=True, fix_final=True, solve_segments=False)
-
-        phase.add_state('y', rate_source='ydot',
-                        units='m',
-                        fix_initial=True, fix_final=True, solve_segments=False)
-
-        phase.add_state('v', rate_source='vdot',
-                        units='m/s',
-                        fix_initial=True, fix_final=False, solve_segments=False)
-
-        phase.add_control('theta', units='deg*s', lower=0.01, upper=179.9, fix_initial=True,
-                          rate_targets=['theta_rate'])
-
-        phase.add_parameter('g', units='m/s**2', opt=False, val=9.80665)
-
-        # Minimize time at the end of the phase
-        phase.add_objective('time', loc='final', scaler=10)
-
-        p.model.linear_solver = om.DirectSolver()
-
-        p.setup()
-
-        phase.set_time_val(initial=0.0, duration=2.0)
-
-        phase.set_state_val('x', [0, 10])
-        phase.set_state_val('y', [10, 5])
-        phase.set_state_val('v', [0, 9.9])
-        phase.set_control_val('theta', [0, 100])
-
-        # Solve for the optimal trajectory
-        p.run_driver()
-
-        # Test the results
-        assert_near_equal(p.get_val('phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-3)
-
-        # Generate the explicitly simulated trajectory
-        exp_out = phase.simulate()
-
-        fig, ax = plt.subplots()
-        fig.suptitle('Brachistochrone Solution')
-
-        x_imp = p.get_val('phase0.timeseries.x')
-        y_imp = p.get_val('phase0.timeseries.y')
-
-        x_exp = exp_out.get_val('phase0.timeseries.x')
-        y_exp = exp_out.get_val('phase0.timeseries.y')
-
-        ax.plot(x_imp, y_imp, 'ro', label='solution')
-        ax.plot(x_exp, y_exp, 'b-', label='simulated')
-
-        ax.set_xlabel('x (m)')
-        ax.set_ylabel('y (m)')
-        ax.grid(True)
-        ax.legend(loc='upper right')
-
-        fig, ax = plt.subplots()
-
-        t_imp = p.get_val('phase0.timeseries.time')
-        theta_imp = p.get_val('phase0.timeseries.theta_rate')
-        t_exp = exp_out.get_val('phase0.timeseries.time')
-        theta_exp = exp_out.get_val('phase0.timeseries.theta_rate')
-
-        ax.plot(t_imp, theta_imp, 'ro', label='solution')
-        ax.plot(t_exp, theta_exp, 'b-', label='simulated')
-
-        ax.set_xlabel('time (s)')
-        ax.set_ylabel(r'$\theta$ (deg)')
-        ax.grid(True)
-        ax.legend(loc='upper right')
-
-        plt.show()
-
-    @unittest.skipIf(matplotlib is None, "This test requires matplotlib")
-    def test_brachistochrone_control_rate_targets_radau(self):
-        import matplotlib.pyplot as plt
-        plt.switch_backend('Agg')
-        import matplotlib.pyplot as plt
-        import openmdao.api as om
-        from openmdao.utils.assert_utils import assert_near_equal
-        import dymos as dm
-
-        p = om.Problem(model=om.Group())
-        p.driver = om.ScipyOptimizeDriver()
-        p.driver.declare_coloring()
-
-        phase = dm.Phase(ode_class=BrachistochroneRateTargetODE,
-                         transcription=dm.Radau(num_segments=10))
-
-        p.model.add_subsystem('phase0', phase)
-
-        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
-
-        phase.add_state('x', rate_source='xdot',
-                        units='m',
-                        fix_initial=True, fix_final=True, solve_segments=False)
-
-        phase.add_state('y', rate_source='ydot',
-                        units='m',
-                        fix_initial=True, fix_final=True, solve_segments=False)
-
-        phase.add_state('v', rate_source='vdot',
-                        units='m/s',
-                        fix_initial=True, fix_final=False, solve_segments=False)
-
-        phase.add_control('theta', units='deg*s', lower=0.01, upper=179.9, fix_initial=True,
-                          rate_targets=['theta_rate'])
-
-        phase.add_parameter('g', units='m/s**2', opt=False, val=9.80665)
-
-        # Minimize time at the end of the phase
-        phase.add_objective('time', loc='final', scaler=10)
-
-        p.model.linear_solver = om.DirectSolver()
-
-        p.setup()
-
-        phase.set_time_val(initial=0.0, duration=2.0)
-
-        phase.set_state_val('x', [0, 10])
-        phase.set_state_val('y', [10, 5])
-        phase.set_state_val('v', [0, 9.9])
-        phase.set_control_val('theta', [0, 100])
-
-        # Solve for the optimal trajectory
-        p.run_driver()
-
-        # Test the results
-        assert_near_equal(p.get_val('phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-3)
-
-        # Generate the explicitly simulated trajectory
-        exp_out = phase.simulate()
-
-        fig, ax = plt.subplots()
-        fig.suptitle('Brachistochrone Solution')
-
-        x_imp = p.get_val('phase0.timeseries.x')
-        y_imp = p.get_val('phase0.timeseries.y')
-
-        x_exp = exp_out.get_val('phase0.timeseries.x')
-        y_exp = exp_out.get_val('phase0.timeseries.y')
-
-        ax.plot(x_imp, y_imp, 'ro', label='solution')
-        ax.plot(x_exp, y_exp, 'b-', label='simulated')
-
-        ax.set_xlabel('x (m)')
-        ax.set_ylabel('y (m)')
-        ax.grid(True)
-        ax.legend(loc='upper right')
-
-        fig, ax = plt.subplots()
-
-        t_imp = p.get_val('phase0.timeseries.time')
-        theta_imp = p.get_val('phase0.timeseries.theta_rate')
-        t_exp = exp_out.get_val('phase0.timeseries.time')
-        theta_exp = exp_out.get_val('phase0.timeseries.theta_rate')
+    def test_brachistochrone_control_rate_targets(self):
 
-        ax.plot(t_imp, theta_imp, 'ro', label='solution')
-        ax.plot(t_exp, theta_exp, 'b-', label='simulated')
-
-        ax.set_xlabel('time (s)')
-        ax.set_ylabel(r'$\theta$ (deg)')
-        ax.grid(True)
-        ax.legend(loc='upper right')
-
-        plt.show()
-
-
-@use_tempdirs
-class TestBrachistochronePolynomialControlRateTargets(unittest.TestCase):
-
-    @unittest.skipIf(matplotlib is None, "This test requires matplotlib")
-    def test_brachistochrone_polynomial_control_rate_targets_gauss_lobatto(self):
-        import matplotlib.pyplot as plt
-        plt.switch_backend('Agg')
-        import openmdao.api as om
-        from openmdao.utils.assert_utils import assert_near_equal
-        import dymos as dm
-
-        p = om.Problem(model=om.Group())
-        p.driver = om.ScipyOptimizeDriver()
-        p.driver.declare_coloring()
-
-        phase = dm.Phase(ode_class=BrachistochroneRateTargetODE,
-                         transcription=dm.GaussLobatto(num_segments=10))
-
-        p.model.add_subsystem('phase0', phase)
-
-        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
-
-        phase.add_state('x', rate_source='xdot',
-                        units='m',
-                        fix_initial=True, fix_final=True, solve_segments=False)
-
-        phase.add_state('y', rate_source='ydot',
-                        units='m',
-                        fix_initial=True, fix_final=True, solve_segments=False)
-
-        phase.add_state('v', rate_source='vdot',
-                        units='m/s',
-                        fix_initial=True, fix_final=False, solve_segments=False)
-
-        phase.add_control('theta', order=3, units='deg*s', lower=0.01, upper=179.9,
-                          fix_initial=True, control_type='polynomial')
-
-        phase.add_parameter('g', units='m/s**2', opt=False, val=9.80665)
-
-        # Minimize time at the end of the phase
-        phase.add_objective('time', loc='final', scaler=10)
-
-        p.model.linear_solver = om.DirectSolver()
-
-        p.setup()
-
-        phase.set_time_val(initial=0.0, duration=2.0)
-
-        phase.set_state_val('x', [0, 10])
-        phase.set_state_val('y', [10, 5])
-        phase.set_state_val('v', [0, 9.9])
-        phase.set_control_val('theta', [0, 100])
-
-        # Solve for the optimal trajectory
-        p.run_driver()
-
-        # Test the results
-        assert_near_equal(p.get_val('phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-3)
-
-        # Generate the explicitly simulated trajectory
-        exp_out = phase.simulate()
-
-        fig, ax = plt.subplots()
-        fig.suptitle('Brachistochrone Solution')
-
-        x_imp = p.get_val('phase0.timeseries.x')
-        y_imp = p.get_val('phase0.timeseries.y')
-
-        x_exp = exp_out.get_val('phase0.timeseries.x')
-        y_exp = exp_out.get_val('phase0.timeseries.y')
-
-        ax.plot(x_imp, y_imp, 'ro', label='solution')
-        ax.plot(x_exp, y_exp, 'b-', label='simulated')
-
-        ax.set_xlabel('x (m)')
-        ax.set_ylabel('y (m)')
-        ax.grid(True)
-        ax.legend(loc='upper right')
-
-        fig, ax = plt.subplots()
-
-        t_imp = p.get_val('phase0.timeseries.time')
-        theta_imp = p.get_val('phase0.timeseries.theta_rate')
-        t_exp = exp_out.get_val('phase0.timeseries.time')
-        theta_exp = exp_out.get_val('phase0.timeseries.theta_rate')
-
-        ax.plot(t_imp, theta_imp, 'ro', label='solution')
-        ax.plot(t_exp, theta_exp, 'b-', label='simulated')
-
-        ax.set_xlabel('time (s)')
-        ax.set_ylabel(r'$\theta$ (deg)')
-        ax.grid(True)
-        ax.legend(loc='upper right')
-
-        plt.show()
-
-    @unittest.skipIf(matplotlib is None, "This test requires matplotlib")
-    def test_brachistochrone_polynomial_control_rate_targets_radau(self):
-        import matplotlib.pyplot as plt
-        plt.switch_backend('Agg')
-        import matplotlib.pyplot as plt
-        import openmdao.api as om
-        from openmdao.utils.assert_utils import assert_near_equal
-        import dymos as dm
-
-        p = om.Problem(model=om.Group())
-        p.driver = om.ScipyOptimizeDriver()
-        p.driver.declare_coloring()
-
-        phase = dm.Phase(ode_class=BrachistochroneRateTargetODE,
-                         transcription=dm.Radau(num_segments=10))
-
-        p.model.add_subsystem('phase0', phase)
-
-        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
-
-        phase.add_state('x', rate_source='xdot',
-                        units='m',
-                        fix_initial=True, fix_final=True, solve_segments=False)
-
-        phase.add_state('y', rate_source='ydot',
-                        units='m',
-                        fix_initial=True, fix_final=True, solve_segments=False)
-
-        phase.add_state('v', rate_source='vdot',
-                        units='m/s',
-                        fix_initial=True, fix_final=False, solve_segments=False)
-
-        phase.add_control('theta', order=3, units='deg*s', lower=0.01, upper=179.9,
-                          fix_initial=True, control_type='polynomial')
-
-        phase.add_parameter('g', units='m/s**2', opt=False, val=9.80665)
-
-        # Minimize time at the end of the phase
-        phase.add_objective('time', loc='final', scaler=10)
-
-        p.model.linear_solver = om.DirectSolver()
-
-        p.setup()
-
-        phase.set_time_val(initial=0.0, duration=2.0)
-
-        phase.set_state_val('x', [0, 10])
-        phase.set_state_val('y', [10, 5])
-        phase.set_state_val('v', [0, 9.9])
-        phase.set_control_val('theta', [0, 100])
-
-        # Solve for the optimal trajectory
-        p.run_driver()
-
-        # Test the results
-        assert_near_equal(p.get_val('phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-3)
-
-        # Generate the explicitly simulated trajectory
-        exp_out = phase.simulate()
-
-        fig, ax = plt.subplots()
-        fig.suptitle('Brachistochrone Solution')
-
-        x_imp = p.get_val('phase0.timeseries.x')
-        y_imp = p.get_val('phase0.timeseries.y')
-
-        x_exp = exp_out.get_val('phase0.timeseries.x')
-        y_exp = exp_out.get_val('phase0.timeseries.y')
-
-        ax.plot(x_imp, y_imp, 'ro', label='solution')
-        ax.plot(x_exp, y_exp, 'b-', label='simulated')
-
-        ax.set_xlabel('x (m)')
-        ax.set_ylabel('y (m)')
-        ax.grid(True)
-        ax.legend(loc='upper right')
-
-        fig, ax = plt.subplots()
-
-        t_imp = p.get_val('phase0.timeseries.time')
-        theta_imp = p.get_val('phase0.timeseries.theta_rate')
-        t_exp = exp_out.get_val('phase0.timeseries.time')
-        theta_exp = exp_out.get_val('phase0.timeseries.theta_rate')
-
-        ax.plot(t_imp, theta_imp, 'ro', label='solution')
-        ax.plot(t_exp, theta_exp, 'b-', label='simulated')
-
-        ax.set_xlabel('time (s)')
-        ax.set_ylabel(r'$\theta$ (deg)')
-        ax.grid(True)
-        ax.legend(loc='upper right')
-
-        plt.show()
-
-
-@use_tempdirs
-class TestBrachistochronePolynomialControlExplicitRateTargets(unittest.TestCase):
-
-    @unittest.skipIf(matplotlib is None, "This test requires matplotlib")
-    def test_brachistochrone_polynomial_control_rate_targets_gauss_lobatto(self):
-        import matplotlib.pyplot as plt
-        plt.switch_backend('Agg')
-        import openmdao.api as om
-        from openmdao.utils.assert_utils import assert_near_equal
-        import dymos as dm
-
-        p = om.Problem(model=om.Group())
-        p.driver = om.ScipyOptimizeDriver()
-        p.driver.declare_coloring()
-
-        phase = dm.Phase(ode_class=BrachistochroneRateTargetODE,
-                         transcription=dm.GaussLobatto(num_segments=10))
-
-        p.model.add_subsystem('phase0', phase)
-
-        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
-
-        phase.add_state('x', rate_source='xdot',
-                        units='m',
-                        fix_initial=True, fix_final=True, solve_segments=False)
-
-        phase.add_state('y', rate_source='ydot',
-                        units='m',
-                        fix_initial=True, fix_final=True, solve_segments=False)
-
-        phase.add_state('v', rate_source='vdot',
-                        units='m/s',
-                        fix_initial=True, fix_final=False, solve_segments=False)
-
-        phase.add_control('theta', order=3, units='deg*s', lower=0.01, upper=179.9,
-                          rate_targets=['theta_rate'], fix_initial=True, control_type='polynomial')
-
-        phase.add_parameter('g', units='m/s**2', opt=False, val=9.80665)
-
-        # Minimize time at the end of the phase
-        phase.add_objective('time', loc='final', scaler=10)
-
-        p.model.linear_solver = om.DirectSolver()
-
-        p.setup()
-
-        phase.set_time_val(initial=0.0, duration=2.0)
-
-        phase.set_state_val('x', [0, 10])
-        phase.set_state_val('y', [10, 5])
-        phase.set_state_val('v', [0, 9.9])
-        phase.set_control_val('theta', [0, 100])
-
-        # Solve for the optimal trajectory
-        p.run_driver()
-
-        # Test the results
-        assert_near_equal(p.get_val('phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-3)
-
-        # Generate the explicitly simulated trajectory
-        exp_out = phase.simulate()
-
-        fig, ax = plt.subplots()
-        fig.suptitle('Brachistochrone Solution')
-
-        x_imp = p.get_val('phase0.timeseries.x')
-        y_imp = p.get_val('phase0.timeseries.y')
-
-        x_exp = exp_out.get_val('phase0.timeseries.x')
-        y_exp = exp_out.get_val('phase0.timeseries.y')
-
-        ax.plot(x_imp, y_imp, 'ro', label='solution')
-        ax.plot(x_exp, y_exp, 'b-', label='simulated')
-
-        ax.set_xlabel('x (m)')
-        ax.set_ylabel('y (m)')
-        ax.grid(True)
-        ax.legend(loc='upper right')
-
-        fig, ax = plt.subplots()
-
-        t_imp = p.get_val('phase0.timeseries.time')
-        theta_imp = p.get_val('phase0.timeseries.theta_rate')
-        t_exp = exp_out.get_val('phase0.timeseries.time')
-        theta_exp = exp_out.get_val('phase0.timeseries.theta_rate')
-
-        ax.plot(t_imp, theta_imp, 'ro', label='solution')
-        ax.plot(t_exp, theta_exp, 'b-', label='simulated')
-
-        ax.set_xlabel('time (s)')
-        ax.set_ylabel(r'$\theta$ (deg)')
-        ax.grid(True)
-        ax.legend(loc='upper right')
-
-        plt.show()
-
-    @unittest.skipIf(matplotlib is None, "This test requires matplotlib")
-    def test_brachistochrone_polynomial_control_rate_targets_radau(self):
-        import matplotlib.pyplot as plt
-        plt.switch_backend('Agg')
-        import matplotlib.pyplot as plt
         import openmdao.api as om
         from openmdao.utils.assert_utils import assert_near_equal
         import dymos as dm
 
-        p = om.Problem(model=om.Group())
-        p.driver = om.ScipyOptimizeDriver()
-        p.driver.declare_coloring()
-
-        phase = dm.Phase(ode_class=BrachistochroneRateTargetODE,
-                         transcription=dm.Radau(num_segments=10))
-
-        p.model.add_subsystem('phase0', phase)
-
-        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
-
-        phase.add_state('x', rate_source='xdot',
-                        units='m',
-                        fix_initial=True, fix_final=True, solve_segments=False)
-
-        phase.add_state('y', rate_source='ydot',
-                        units='m',
-                        fix_initial=True, fix_final=True, solve_segments=False)
-
-        phase.add_state('v', rate_source='vdot',
-                        units='m/s',
-                        fix_initial=True, fix_final=False, solve_segments=False)
-
-        phase.add_control('theta', order=3, units='deg*s', lower=0.01, upper=179.9,
-                          rate_targets=['theta_rate'], fix_initial=True, control_type='polynomial')
+        transcriptions = {'gauss-lobatto': dm.GaussLobatto(num_segments=10),
+                          'radau': dm.Radau(num_segments=10),
+                          'birkhoff': dm.Birkhoff(num_nodes=30)}
 
-        phase.add_parameter('g', units='m/s**2', opt=False, val=9.80665)
+        for tx_name, tx in transcriptions.items():
+            for control_target_method in ('implicit', 'explicit'):
+                for control_type in ('full', 'polynomial'):
 
-        # Minimize time at the end of the phase
-        phase.add_objective('time', loc='final', scaler=10)
+                    with self.subTest(f'{tx_name=} {control_target_method=} {control_type=}'):
 
-        p.model.linear_solver = om.DirectSolver()
+                        p = om.Problem(model=om.Group())
+                        p.driver = om.ScipyOptimizeDriver()
+                        p.driver.declare_coloring()
 
-        p.setup()
+                        traj = dm.Trajectory()
 
-        phase.set_time_val(initial=0.0, duration=2.0)
+                        ode_kwargs = {'control_name': 'int_theta_rate' if control_target_method == 'implicit' else 'theta'}
 
-        phase.set_state_val('x', [0, 10])
-        phase.set_state_val('y', [10, 5])
-        phase.set_state_val('v', [0, 9.9])
-        phase.set_control_val('theta', [0, 100])
+                        phase = dm.Phase(ode_class=BrachistochroneRateTargetODE,
+                                         ode_init_kwargs=ode_kwargs,
+                                         transcription=tx)
 
-        # Solve for the optimal trajectory
-        p.run_driver()
+                        traj.add_phase('phase0', phase)
 
-        # Test the results
-        assert_near_equal(p.get_val('phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-3)
+                        p.model.add_subsystem('traj0', traj)
 
-        # Generate the explicitly simulated trajectory
-        exp_out = phase.simulate()
+                        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
 
-        fig, ax = plt.subplots()
-        fig.suptitle('Brachistochrone Solution')
+                        phase.add_state('x', rate_source='xdot',
+                                        units='m',
+                                        fix_initial=True, fix_final=True, solve_segments=False)
 
-        x_imp = p.get_val('phase0.timeseries.x')
-        y_imp = p.get_val('phase0.timeseries.y')
-
-        x_exp = exp_out.get_val('phase0.timeseries.x')
-        y_exp = exp_out.get_val('phase0.timeseries.y')
-
-        ax.plot(x_imp, y_imp, 'ro', label='solution')
-        ax.plot(x_exp, y_exp, 'b-', label='simulated')
-
-        ax.set_xlabel('x (m)')
-        ax.set_ylabel('y (m)')
-        ax.grid(True)
-        ax.legend(loc='upper right')
-
-        fig, ax = plt.subplots()
-
-        t_imp = p.get_val('phase0.timeseries.time')
-        theta_imp = p.get_val('phase0.timeseries.theta_rate')
-        t_exp = exp_out.get_val('phase0.timeseries.time')
-        theta_exp = exp_out.get_val('phase0.timeseries.theta_rate')
-
-        ax.plot(t_imp, theta_imp, 'ro', label='solution')
-        ax.plot(t_exp, theta_exp, 'b-', label='simulated')
-
-        ax.set_xlabel('time (s)')
-        ax.set_ylabel(r'$\theta$ (deg)')
-        ax.grid(True)
-        ax.legend(loc='upper right')
-
-        plt.show()
-
-
-@use_tempdirs
-class TestBrachistochronePolynomialControlExplicitRate2Targets(unittest.TestCase):
-
-    @unittest.skipIf(matplotlib is None, "This test requires matplotlib")
-    def test_brachistochrone_polynomial_control_rate_targets_gauss_lobatto(self):
-        import matplotlib.pyplot as plt
-        plt.switch_backend('Agg')
-        import openmdao.api as om
-        from openmdao.utils.assert_utils import assert_near_equal
-        import dymos as dm
+                        phase.add_state('y', rate_source='ydot',
+                                        units='m',
+                                        fix_initial=True, fix_final=True, solve_segments=False)
 
-        p = om.Problem(model=om.Group())
-        p.driver = om.ScipyOptimizeDriver()
-        p.driver.declare_coloring()
+                        phase.add_state('v', rate_source='vdot',
+                                        units='m/s',
+                                        fix_initial=True, fix_final=False, solve_segments=False)
 
-        phase = dm.Phase(ode_class=BrachistochroneRateTargetODE,
-                         transcription=dm.GaussLobatto(num_segments=10))
+                        phase.add_control('int_theta', lower=0.0, upper=None, fix_initial=True,
+                                          rate_targets=_unspecified if control_target_method == 'implicit' else ['theta'],
+                                          control_type=control_type, order=7)
 
-        p.model.add_subsystem('phase0', phase)
+                        phase.add_parameter('g', units='m/s**2', opt=False, val=9.80665)
 
-        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+                        # Minimize time at the end of the phase
+                        phase.add_objective('time', loc='final', scaler=10)
 
-        phase.add_state('x', rate_source='xdot',
-                        units='m',
-                        fix_initial=True, fix_final=True, solve_segments=False)
+                        p.model.linear_solver = om.DirectSolver()
 
-        phase.add_state('y', rate_source='ydot',
-                        units='m',
-                        fix_initial=True, fix_final=True, solve_segments=False)
+                        phase.timeseries_options['include_control_rates'] = True
 
-        phase.add_state('v', rate_source='vdot',
-                        units='m/s',
-                        fix_initial=True, fix_final=False, solve_segments=False)
+                        p.setup()
 
-        phase.add_control('theta', order=5, units='deg*s**2', lower=0.01, upper=179.9,
-                          rate_targets=None, rate2_targets=['theta_rate'], fix_initial=True, control_type='polynomial')
+                        phase.set_time_val(initial=0.0, duration=2.0)
 
-        phase.add_parameter('g', units='m/s**2', opt=False, val=9.80665)
+                        phase.set_state_val('x', [0, 10])
+                        phase.set_state_val('y', [10, 5])
+                        phase.set_state_val('v', [0, 9.9])
+                        phase.set_control_val('int_theta', [0, 100], units='deg*s')
 
-        # Minimize time at the end of the phase
-        phase.add_objective('time', loc='final', scaler=10)
+                        p.run_model()
 
-        p.model.linear_solver = om.DirectSolver()
+                        # Solve for the optimal trajectory
+                        dm.run_problem(p, simulate=True, make_plots=True)
 
-        p.setup()
+                        sol_db = 'dymos_solution.db'
+                        sim_db = 'dymos_simulation.db'
+                        if om_version()[0] > (3, 34, 2):
+                            sol_db = p.get_outputs_dir() / sol_db
+                            sim_db = traj.sim_prob.get_outputs_dir() / sim_db
 
-        phase.set_time_val(initial=0.0, duration=2.0)
+                        sol_case = om.CaseReader(sol_db).get_case('final')
+                        sim_case = om.CaseReader(sim_db).get_case('final')
 
-        phase.set_state_val('x', [0, 10])
-        phase.set_state_val('y', [10, 5])
-        phase.set_state_val('v', [0, 9.9])
-        phase.set_control_val('theta', [0, 10, 40, 60, 80, 100],
-                              time_vals=[0, 0.4, 0.8, 1.2, 1.6, 2.0])
+                        t_sol = sol_case.get_val('traj0.phase0.timeseries.time')
+                        x_sol = sol_case.get_val('traj0.phase0.timeseries.x')
+                        y_sol = sol_case.get_val('traj0.phase0.timeseries.y')
+                        v_sol = sol_case.get_val('traj0.phase0.timeseries.v')
+                        theta_sol = sol_case.get_val('traj0.phase0.timeseries.int_theta_rate')
 
-        # Solve for the optimal trajectory
-        dm.run_problem(p, simulate=True)
+                        t_sim = sim_case.get_val('traj0.phase0.timeseries.time')
+                        x_sim = sim_case.get_val('traj0.phase0.timeseries.x')
+                        y_sim = sim_case.get_val('traj0.phase0.timeseries.y')
+                        v_sim = sim_case.get_val('traj0.phase0.timeseries.v')
+                        theta_sim = sim_case.get_val('traj0.phase0.timeseries.int_theta_rate')
 
-        # Test the results
-        assert_near_equal(p.get_val('phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-3)
+                        assert_timeseries_near_equal(t_sol, x_sol, t_sim, x_sim, rel_tolerance=5.0E-3, abs_tolerance=0.01)
+                        assert_timeseries_near_equal(t_sol, y_sol, t_sim, y_sim, rel_tolerance=5.0E-3, abs_tolerance=0.01)
+                        assert_timeseries_near_equal(t_sol, v_sol, t_sim, v_sim, rel_tolerance=5.0E-3, abs_tolerance=0.01)
+                        assert_timeseries_near_equal(t_sol, theta_sol, t_sim, theta_sim, rel_tolerance=5.0E-3, abs_tolerance=0.01)
 
 
 if __name__ == '__main__':  # pragma: no cover