Repository: OpenMDAO/dymos (python)

Issue:
Polynomial controls don't work with Birkhoff transcription
### Description

Attempting to use a polynomial control causes an error during configure. This occurs only in the Birkhoff transcription and only with polynomial controls
`AttributeError: 'Phase' object has no attribute 'ode_all'`

### Example

https://github.com/OpenMDAO/dymos/blob/master/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py

### Dymos Version

1.9.2.-dev

### Relevant environment information

_No response_

Pull request description:
Birkhoff Polynomial Controls
### Summary

Fixed the configure function of polynomial controls in the birkhoff transcription to connect the right values to ode_all and boundary_vals

### Related Issues

- Resolves #1019 

### Backwards incompatibilities

None

### New Dependencies

None


Patch:
diff --git a/dymos/transcriptions/pseudospectral/birkhoff.py b/dymos/transcriptions/pseudospectral/birkhoff.py
index 08d3aabb..df914ead 100644
--- a/dymos/transcriptions/pseudospectral/birkhoff.py
+++ b/dymos/transcriptions/pseudospectral/birkhoff.py
@@ -211,34 +211,24 @@ class Birkhoff(TranscriptionBase):
         """
         super(Birkhoff, self).configure_polynomial_controls(phase)
 
-        ode_inputs = get_promoted_vars(self._get_ode(phase), 'input')
-
         for name, options in phase.polynomial_control_options.items():
-            targets = get_targets(ode=ode_inputs, name=name, user_targets=options['targets'])
-            if targets:
-                phase.connect(f'polynomial_control_values:{name}',
-                              [f'ode_all.{t}' for t in targets])
-                phase.connect(f'polynomial_control_values:{name}',
-                              [f'boundary_vals.{t}' for t in targets],
+            if options['targets']:
+                phase.connect(f'polynomial_control_values:{name}', [f'ode_all.{t}' for t in options['targets']])
+                phase.connect(f'polynomial_control_values:{name}', [f'boundary_vals.{t}' for t in options['targets']],
                               src_indices=om.slicer[[0, -1], ...])
 
-            targets = get_targets(ode=phase.ode_all, name=f'{name}_rate',
-                                  user_targets=options['rate_targets'])
-            if targets:
+            if options['rate_targets']:
                 phase.connect(f'polynomial_control_rates:{name}_rate',
-                              [f'ode_all.{t}' for t in targets])
+                              [f'ode_all.{t}' for t in options['rate_targets']])
                 phase.connect(f'polynomial_control_rates:{name}_rate',
-                              [f'boundary_vals.{t}' for t in targets],
+                              [f'boundary_vals.{t}' for t in options['rate_targets']],
                               src_indices=om.slicer[[0, -1], ...])
 
-            targets = get_targets(ode=phase.ode_all, name=f'{name}_rate2',
-                                  user_targets=options['rate2_targets'],
-                                  src_indices=om.slicer[[0, -1], ...])
-            if targets:
+            if options['rate2_targets']:
                 phase.connect(f'polynomial_control_rates:{name}_rate2',
-                              [f'ode_all.{t}' for t in targets])
+                              [f'ode_all.{t}' for t in options['rate2_targets']])
                 phase.connect(f'polynomial_control_rates:{name}_rate2',
-                              [f'boundary_vals.{t}' for t in targets],
+                              [f'boundary_vals.{t}' for t in options['rate2_targets']],
                               src_indices=om.slicer[[0, -1], ...])
 
     def setup_ode(self, phase):


Test patch:
diff --git a/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py b/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py
index d5a87e48..e3047f54 100644
--- a/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py
+++ b/dymos/examples/brachistochrone/test/test_doc_brachistochrone_polynomial_controls.py
@@ -177,6 +177,99 @@ class TestBrachistochronePolynomialControl(unittest.TestCase):
 
         plt.show()
 
+    def test_brachistochrone_polynomial_control_birkhoff(self):
+        import numpy as np
+        import matplotlib
+        matplotlib.use('Agg')
+        import matplotlib.pyplot as plt
+        import openmdao.api as om
+        from openmdao.utils.assert_utils import assert_near_equal
+        import dymos as dm
+        from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE
+
+        p = om.Problem(model=om.Group())
+        p.driver = om.ScipyOptimizeDriver()
+        p.driver.declare_coloring()
+
+        phase = dm.Phase(ode_class=BrachistochroneODE,
+                         transcription=dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=1, nodes_per_seg=15)))
+
+        p.model.add_subsystem('phase0', phase)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+
+        phase.add_state('x', fix_initial=True, fix_final=True)
+
+        phase.add_state('y', fix_initial=True, fix_final=True)
+
+        phase.add_state('v', fix_initial=True, fix_final=False)
+
+        phase.add_polynomial_control('theta', order=1, units='deg', lower=0.01, upper=179.9)
+
+        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()
+
+        p['phase0.t_initial'] = 0.0
+        p['phase0.t_duration'] = 2.0
+
+        p.set_val('phase0.initial_states:x', 0.0)
+        p.set_val('phase0.initial_states:y', 10.0)
+        p.set_val('phase0.initial_states:v', 0.0)
+        p.set_val('phase0.final_states:x', 10.0)
+        p.set_val('phase0.final_states:y', 5.0)
+        p.set_val('phase0.final_states:v', 9.9)
+
+        p.set_val('phase0.states:x', phase.interp('x', [0, 10]))
+        p.set_val('phase0.states:y', phase.interp('y', [10, 5]))
+        p.set_val('phase0.states:v', phase.interp('v', [0, 9.9]))
+        p.set_val('phase0.polynomial_controls:theta', phase.interp('theta', [5, 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')
+
+        ax.plot(t_imp, theta_imp, 'ro', label='solution')
+
+        ax.set_xlabel('time (s)')
+        ax.set_ylabel(r'$\theta$ (deg)')
+        ax.grid(True)
+        ax.legend(loc='upper right')
+
+        plt.show()
+
 
 @use_tempdirs
 class TestBrachistochronePolynomialControlBoundaryConstrained(unittest.TestCase):
@@ -357,6 +450,102 @@ class TestBrachistochronePolynomialControlBoundaryConstrained(unittest.TestCase)
 
         plt.show()
 
+    def test_brachistochrone_polynomial_control_birkhoff(self):
+        import numpy as np
+        import matplotlib
+        matplotlib.use('Agg')
+        import matplotlib.pyplot as plt
+        import openmdao.api as om
+        from openmdao.utils.assert_utils import assert_near_equal
+        import dymos as dm
+        from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE
+
+        p = om.Problem(model=om.Group())
+        p.driver = om.ScipyOptimizeDriver()
+        p.driver.declare_coloring()
+
+        phase = dm.Phase(ode_class=BrachistochroneODE,
+                         transcription=dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=3, nodes_per_seg=5)))
+
+        p.model.add_subsystem('phase0', phase)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+
+        phase.add_state('x', fix_initial=True, fix_final=True)
+
+        phase.add_state('y', fix_initial=True, fix_final=True)
+
+        phase.add_state('v', fix_initial=True, fix_final=False)
+
+        phase.add_polynomial_control('theta', order=1, units='deg', lower=0.01, upper=179.9)
+
+        phase.add_parameter('g', units='m/s**2', opt=False, val=9.80665)
+
+        phase.add_boundary_constraint('theta', loc='initial', lower=0, upper=1.0)
+        phase.add_boundary_constraint('theta', loc='final', upper=105.0)
+
+        # Minimize time at the end of the phase
+        phase.add_objective('time', loc='final', scaler=10)
+
+        p.model.linear_solver = om.DirectSolver()
+
+        p.setup()
+
+        p['phase0.t_initial'] = 0.0
+        p['phase0.t_duration'] = 2.0
+
+        p.set_val('phase0.initial_states:x', 0.0)
+        p.set_val('phase0.initial_states:y', 10.0)
+        p.set_val('phase0.initial_states:v', 0.0)
+        p.set_val('phase0.final_states:x', 10.0)
+        p.set_val('phase0.final_states:y', 5.0)
+        p.set_val('phase0.final_states:v', 9.9)
+
+        p.set_val('phase0.states:x', phase.interp('x', [0, 10]))
+        p.set_val('phase0.states:y', phase.interp('y', [10, 5]))
+        p.set_val('phase0.states:v', phase.interp('v', [0, 9.9]))
+        p.set_val('phase0.polynomial_controls:theta', phase.interp('theta', [1, 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')
+
+        ax.plot(t_imp, theta_imp, 'ro', label='solution')
+
+        ax.set_xlabel('time (s)')
+        ax.set_ylabel(r'$\theta$ (deg)')
+        ax.grid(True)
+        ax.legend(loc='upper right')
+
+        plt.show()
+
 
 @use_tempdirs
 class TestBrachistochronePolynomialControlPathConstrained(unittest.TestCase):
@@ -533,6 +722,101 @@ class TestBrachistochronePolynomialControlPathConstrained(unittest.TestCase):
 
         plt.show()
 
+    def test_brachistochrone_polynomial_control_birkhoff(self):
+        import numpy as np
+        import matplotlib
+        matplotlib.use('Agg')
+        import matplotlib.pyplot as plt
+        import openmdao.api as om
+        from openmdao.utils.assert_utils import assert_near_equal
+        import dymos as dm
+        from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE
+
+        p = om.Problem(model=om.Group())
+        p.driver = om.ScipyOptimizeDriver()
+        p.driver.declare_coloring()
+
+        phase = dm.Phase(ode_class=BrachistochroneODE,
+                         transcription=dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=1, nodes_per_seg=15)))
+
+        p.model.add_subsystem('phase0', phase)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+
+        phase.add_state('x', fix_initial=True, fix_final=True)
+
+        phase.add_state('y', fix_initial=True, fix_final=True)
+
+        phase.add_state('v', fix_initial=True, fix_final=False)
+
+        phase.add_polynomial_control('theta', order=1, units='deg', lower=0.01, upper=179.9)
+
+        phase.add_parameter('g', units='m/s**2', opt=False, val=9.80665)
+
+        phase.add_path_constraint('theta', lower=1, upper=120)
+
+        # Minimize time at the end of the phase
+        phase.add_objective('time', loc='final', scaler=10)
+
+        p.model.linear_solver = om.DirectSolver()
+
+        p.setup()
+
+        p['phase0.t_initial'] = 0.0
+        p['phase0.t_duration'] = 2.0
+
+        p.set_val('phase0.initial_states:x', 0.0)
+        p.set_val('phase0.initial_states:y', 10.0)
+        p.set_val('phase0.initial_states:v', 0.0)
+        p.set_val('phase0.final_states:x', 10.0)
+        p.set_val('phase0.final_states:y', 5.0)
+        p.set_val('phase0.final_states:v', 9.9)
+
+        p.set_val('phase0.states:x', phase.interp('x', [0, 10]))
+        p.set_val('phase0.states:y', phase.interp('y', [10, 5]))
+        p.set_val('phase0.states:v', phase.interp('v', [0, 9.9]))
+        p.set_val('phase0.polynomial_controls:theta', phase.interp('theta', [5, 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')
+
+        ax.plot(t_imp, theta_imp, 'ro', label='solution')
+
+        ax.set_xlabel('time (s)')
+        ax.set_ylabel(r'$\theta$ (deg)')
+        ax.grid(True)
+        ax.legend(loc='upper right')
+
+        plt.show()
+
 
 @use_tempdirs
 class TestBrachistochronePolynomialControlRatePathConstrained(unittest.TestCase):
@@ -713,6 +997,101 @@ class TestBrachistochronePolynomialControlRatePathConstrained(unittest.TestCase)
 
         plt.show()
 
+    def test_brachistochrone_polynomial_control_birkhoff(self):
+        import numpy as np
+        import matplotlib
+        matplotlib.use('Agg')
+        import matplotlib.pyplot as plt
+        import openmdao.api as om
+        from openmdao.utils.assert_utils import assert_near_equal
+        import dymos as dm
+        from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE
+
+        p = om.Problem(model=om.Group())
+        p.driver = om.ScipyOptimizeDriver()
+        p.driver.declare_coloring()
+
+        phase = dm.Phase(ode_class=BrachistochroneODE,
+                         transcription=dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=1, nodes_per_seg=15)))
+
+        p.model.add_subsystem('phase0', phase)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+
+        phase.add_state('x', fix_initial=True, fix_final=True)
+
+        phase.add_state('y', fix_initial=True, fix_final=True)
+
+        phase.add_state('v', fix_initial=True, fix_final=False)
+
+        phase.add_polynomial_control('theta', order=1, units='deg', lower=0.01, upper=179.9)
+
+        phase.add_parameter('g', units='m/s**2', opt=False, val=9.80665)
+
+        phase.add_path_constraint('theta_rate', lower=0, upper=120)
+
+        # Minimize time at the end of the phase
+        phase.add_objective('time', loc='final', scaler=10)
+
+        p.model.linear_solver = om.DirectSolver()
+
+        p.setup()
+
+        p['phase0.t_initial'] = 0.0
+        p['phase0.t_duration'] = 2.0
+
+        p.set_val('phase0.initial_states:x', 0.0)
+        p.set_val('phase0.initial_states:y', 10.0)
+        p.set_val('phase0.initial_states:v', 0.0)
+        p.set_val('phase0.final_states:x', 10.0)
+        p.set_val('phase0.final_states:y', 5.0)
+        p.set_val('phase0.final_states:v', 9.9)
+
+        p.set_val('phase0.states:x', phase.interp('x', [0, 10]))
+        p.set_val('phase0.states:y', phase.interp('y', [10, 5]))
+        p.set_val('phase0.states:v', phase.interp('v', [0, 9.9]))
+        p.set_val('phase0.polynomial_controls:theta', phase.interp('theta', [5, 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')
+
+        ax.plot(t_imp, theta_imp, 'ro', label='solution')
+
+        ax.set_xlabel('time (s)')
+        ax.set_ylabel(r'$\theta$ (deg)')
+        ax.grid(True)
+        ax.legend(loc='upper right')
+
+        plt.show()
+
 
 @use_tempdirs
 class TestBrachistochronePolynomialControlRate2PathConstrained(unittest.TestCase):
@@ -893,6 +1272,101 @@ class TestBrachistochronePolynomialControlRate2PathConstrained(unittest.TestCase
 
         plt.show()
 
+    def test_brachistochrone_polynomial_control_birkhoff(self):
+        import numpy as np
+        import matplotlib
+        matplotlib.use('Agg')
+        import matplotlib.pyplot as plt
+        import openmdao.api as om
+        from openmdao.utils.assert_utils import assert_near_equal
+        import dymos as dm
+        from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE
+
+        p = om.Problem(model=om.Group())
+        p.driver = om.ScipyOptimizeDriver()
+        p.driver.declare_coloring()
+
+        phase = dm.Phase(ode_class=BrachistochroneODE,
+                         transcription=dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=1, nodes_per_seg=15)))
+
+        p.model.add_subsystem('phase0', phase)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+
+        phase.add_state('x', fix_initial=True, fix_final=True)
+
+        phase.add_state('y', fix_initial=True, fix_final=True)
+
+        phase.add_state('v', fix_initial=True, fix_final=False)
+
+        phase.add_polynomial_control('theta', order=1, units='deg', lower=0.01, upper=179.9)
+
+        phase.add_parameter('g', units='m/s**2', opt=False, val=9.80665)
+
+        phase.add_path_constraint('theta_rate2', lower=-0.01, upper=0.01)
+
+        # Minimize time at the end of the phase
+        phase.add_objective('time', loc='final', scaler=10)
+
+        p.model.linear_solver = om.DirectSolver()
+
+        p.setup()
+
+        p['phase0.t_initial'] = 0.0
+        p['phase0.t_duration'] = 2.0
+
+        p.set_val('phase0.initial_states:x', 0.0)
+        p.set_val('phase0.initial_states:y', 10.0)
+        p.set_val('phase0.initial_states:v', 0.0)
+        p.set_val('phase0.final_states:x', 10.0)
+        p.set_val('phase0.final_states:y', 5.0)
+        p.set_val('phase0.final_states:v', 9.9)
+
+        p.set_val('phase0.states:x', phase.interp('x', [0, 10]))
+        p.set_val('phase0.states:y', phase.interp('y', [10, 5]))
+        p.set_val('phase0.states:v', phase.interp('v', [0, 9.9]))
+        p.set_val('phase0.polynomial_controls:theta', phase.interp('theta', [5, 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')
+
+        ax.plot(t_imp, theta_imp, 'ro', label='solution')
+
+        ax.set_xlabel('time (s)')
+        ax.set_ylabel(r'$\theta$ (deg)')
+        ax.grid(True)
+        ax.legend(loc='upper right')
+
+        plt.show()
+
 
 @use_tempdirs
 class TestBrachistochronePolynomialControlSimulation(unittest.TestCase):
@@ -1011,6 +1485,74 @@ class TestBrachistochronePolynomialControlSimulation(unittest.TestCase):
         assert_near_equal(theta_exp[0], theta_imp[0])
         assert_near_equal(theta_exp[-1], theta_imp[-1])
 
+    def test_brachistochrone_polynomial_control_birkhoff(self):
+        import numpy as np
+        import matplotlib
+        matplotlib.use('Agg')
+        import matplotlib.pyplot as plt
+        import openmdao.api as om
+        from openmdao.utils.assert_utils import assert_near_equal
+        import dymos as dm
+        from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE
+
+        p = om.Problem(model=om.Group())
+        p.driver = om.ScipyOptimizeDriver()
+        p.driver.declare_coloring()
+
+        phase = dm.Phase(ode_class=BrachistochroneODE,
+                         transcription=dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=1, nodes_per_seg=15)))
+
+        p.model.add_subsystem('phase0', phase)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+
+        phase.add_state('x', fix_initial=True, fix_final=True)
+
+        phase.add_state('y', fix_initial=True, fix_final=True)
+
+        phase.add_state('v', fix_initial=True, fix_final=False)
+
+        phase.add_polynomial_control('theta', order=1, units='deg', lower=0.01, upper=179.9)
+
+        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()
+
+        p['phase0.t_initial'] = 0.0
+        p['phase0.t_duration'] = 2.0
+
+        p.set_val('phase0.initial_states:x', 0.0)
+        p.set_val('phase0.initial_states:y', 10.0)
+        p.set_val('phase0.initial_states:v', 0.0)
+        p.set_val('phase0.final_states:x', 10.0)
+        p.set_val('phase0.final_states:y', 5.0)
+        p.set_val('phase0.final_states:v', 9.9)
+
+        p.set_val('phase0.states:x', phase.interp('x', [0, 10]))
+        p.set_val('phase0.states:y', phase.interp('y', [10, 5]))
+        p.set_val('phase0.states:v', phase.interp('v', [0, 9.9]))
+        p.set_val('phase0.polynomial_controls:theta', phase.interp('theta', [5, 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()
+
+        theta_imp = p.get_val('phase0.timeseries.theta')
+        theta_exp = exp_out.get_val('phase0.timeseries.theta')
+
+        assert_near_equal(theta_exp[0], theta_imp[0])
+        assert_near_equal(theta_exp[-1], theta_imp[-1])
+
 
 if __name__ == '__main__':  # pragma: no cover
     unittest.main()