Repository: OpenMDAO/dymos
Pull request #1008 (2023-10-17 18:15:33)

Birkhoff Transcription Take 2
### Summary

Added the new birkhoff transcription as laid out by Mike Ross in https://arxiv.org/abs/2308.01400v2
This is currently experimental as we get it cleaned up with a few years worth of new knowledge in how to better implement pseudospectral systems under our belts.

This PR supersedes #1000.

### Related Issues

- Closes #1000

### Backwards incompatibilities

OpenMDAO minimum version is now 3.27.0

### New Dependencies

None


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

coverage: 92.195% (-0.6%) from 92.779% when pulling **a5c9a1603abb5b3e3bcb4ee66c956008d0f220bb on robfalck:kaushik_birkhoff** into **45ab30cac3cf658fe40715d081dcdb586f00da25 on OpenMDAO:master**.


Patch:
diff --git a/dymos/__init__.py b/dymos/__init__.py
index c8c900b0..1f09b912 100644
--- a/dymos/__init__.py
+++ b/dymos/__init__.py
@@ -2,8 +2,8 @@ __version__ = '1.9.2-dev'
 
 
 from .phase import Phase, AnalyticPhase
-from .transcriptions import GaussLobatto, Radau, ExplicitShooting, Analytic
-from .transcriptions.grid_data import GaussLobattoGrid, RadauGrid, UniformGrid
+from .transcriptions import GaussLobatto, Radau, ExplicitShooting, Analytic, Birkhoff
+from .transcriptions.grid_data import GaussLobattoGrid, RadauGrid, UniformGrid, BirkhoffGrid
 from .trajectory.trajectory import Trajectory
 from .run_problem import run_problem
 from .load_case import load_case
diff --git a/dymos/examples/double_integrator/breakwell_ode.py b/dymos/examples/double_integrator/breakwell_ode.py
new file mode 100644
index 00000000..4a42c650
--- /dev/null
+++ b/dymos/examples/double_integrator/breakwell_ode.py
@@ -0,0 +1,25 @@
+import openmdao.api as om
+import numpy as np
+
+
+class BreakwellODE(om.ExplicitComponent):
+
+    def initialize(self):
+        self.options.declare('num_nodes', types=int)
+
+    def setup(self):
+        nn = self.options['num_nodes']
+
+        self.add_input('u', val=np.zeros(nn))
+        self.add_output('J_dot', val=np.zeros(nn))
+
+        ar = np.arange(nn)
+        self.declare_partials(of='J_dot', wrt='u', rows=ar, cols=ar)
+
+    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
+        u = inputs['u']
+        outputs['J_dot'] = 0.5*u**2
+
+    def compute_partials(self, inputs, partials, discrete_inputs=None):
+        u = inputs['u']
+        partials['J_dot', 'u'] = u
diff --git a/dymos/examples/goddard_rocket_problem/__init__.py b/dymos/examples/goddard_rocket_problem/__init__.py
new file mode 100644
index 00000000..bb466fa9
--- /dev/null
+++ b/dymos/examples/goddard_rocket_problem/__init__.py
@@ -0,0 +1,1 @@
+from .rocket_ode import RocketODE
diff --git a/dymos/examples/goddard_rocket_problem/rocket_ode.py b/dymos/examples/goddard_rocket_problem/rocket_ode.py
new file mode 100644
index 00000000..4acb4591
--- /dev/null
+++ b/dymos/examples/goddard_rocket_problem/rocket_ode.py
@@ -0,0 +1,57 @@
+import openmdao.api as om
+import numpy as np
+
+
+class RocketODE(om.ExplicitComponent):
+    def initialize(self):
+        self.options.declare('num_nodes', types=int)
+
+    def setup(self):
+        nn = self.options['num_nodes']
+
+        self.add_input('h', val=np.zeros(nn), units='ft')
+        self.add_input('v', val=np.zeros(nn), units='ft/s')
+        self.add_input('m', val=np.zeros(nn), units='slug')
+        self.add_input('T', val=np.zeros(nn), units='lbf')
+
+        self.add_output('h_dot', val=np.zeros(nn), units='ft/s')
+        self.add_output('v_dot', val=np.zeros(nn), units='ft/s**2')
+        self.add_output('m_dot', val=np.zeros(nn), units='slug/s')
+
+        ar = np.arange(nn)
+        self.declare_partials(of='h_dot', wrt='v', rows=ar, cols=ar, val=1.0)
+        self.declare_partials(of='v_dot', wrt='h', rows=ar, cols=ar)
+        self.declare_partials(of='v_dot', wrt='v', rows=ar, cols=ar)
+        self.declare_partials(of='v_dot', wrt='m', rows=ar, cols=ar)
+        self.declare_partials(of='v_dot', wrt='T', rows=ar, cols=ar)
+
+        self.declare_partials(of='m_dot', wrt='T', rows=ar, cols=ar, val=-1/1580.9425)
+
+    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
+        h = inputs['h']
+        v = inputs['v']
+        m = inputs['m']
+        T = inputs['T']
+
+        sigma = 5.4915e-5
+        g = 32.174
+        c = 1580.9425
+        h0 = 23800
+
+        outputs['h_dot'] = v
+        outputs['v_dot'] = (T - sigma * v**2 * np.exp(-h/h0))/m - g
+        outputs['m_dot'] = -T / c
+
+    def compute_partials(self, inputs, partials, discrete_inputs=None):
+        h = inputs['h']
+        v = inputs['v']
+        m = inputs['m']
+        T = inputs['T']
+
+        sigma = 5.4915e-5
+        h0 = 23800
+
+        partials['v_dot', 'h'] = (sigma * v**2 * np.exp(-h/h0))/(m * h0)
+        partials['v_dot', 'v'] = -(2 * sigma * v * np.exp(-h/h0)) / m
+        partials['v_dot', 'm'] = -(T - sigma * v**2 * np.exp(-h/h0)) / (m**2)
+        partials['v_dot', 'T'] = 1 / m
diff --git a/dymos/examples/low_thrust_spiral/__init__.py b/dymos/examples/low_thrust_spiral/__init__.py
new file mode 100644
index 00000000..654871c8
--- /dev/null
+++ b/dymos/examples/low_thrust_spiral/__init__.py
@@ -0,0 +1,1 @@
+from .low_thrust_ode import LowThrustODE
diff --git a/dymos/examples/low_thrust_spiral/low_thrust_ode.py b/dymos/examples/low_thrust_spiral/low_thrust_ode.py
new file mode 100644
index 00000000..89d0e1cb
--- /dev/null
+++ b/dymos/examples/low_thrust_spiral/low_thrust_ode.py
@@ -0,0 +1,67 @@
+import openmdao.api as om
+import numpy as np
+import openmdao.utils.units as units
+
+
+class LowThrustODE(om.ExplicitComponent):
+    def initialize(self):
+        self.options.declare('num_nodes', types=int)
+
+    def setup(self):
+        nn = self.options['num_nodes']
+
+        self.add_input('r', val=np.zeros(nn), units=None)
+        self.add_input('theta', val=np.zeros(nn), units='rad')
+        self.add_input('vr', val=np.zeros(nn), units='1/s')
+        self.add_input('vt', val=np.zeros(nn), units='1/s')
+        self.add_input('alpha', val=np.zeros(nn), units='rad')
+
+        self.add_output('theta_dot', val=np.zeros(nn), units='rad/s')
+        self.add_output('vr_dot', val=np.zeros(nn), units='1/s**2')
+        self.add_output('vt_dot', val=np.zeros(nn), units='1/s**2')
+
+        ar = np.arange(nn)
+
+        self.declare_partials(of='theta_dot', wrt='r', rows=ar, cols=ar)
+        self.declare_partials(of='theta_dot', wrt='vt', rows=ar, cols=ar)
+
+        self.declare_partials(of='vr_dot', wrt='r', rows=ar, cols=ar)
+        self.declare_partials(of='vr_dot', wrt='vt', rows=ar, cols=ar)
+        self.declare_partials(of='vr_dot', wrt='alpha', rows=ar, cols=ar)
+
+        self.declare_partials(of='vt_dot', wrt='r', rows=ar, cols=ar)
+        self.declare_partials(of='vt_dot', wrt='vr', rows=ar, cols=ar)
+        self.declare_partials(of='vt_dot', wrt='vt', rows=ar, cols=ar)
+        self.declare_partials(of='vt_dot', wrt='alpha', rows=ar, cols=ar)
+
+    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
+        r = inputs['r']
+        vr = inputs['vr']
+        vt = inputs['vt']
+        alpha = inputs['alpha']
+
+        T = 1e-3
+
+        outputs['theta_dot'] = vt/r
+        outputs['vr_dot'] = vt**2/r - 1/r**2 + T*np.sin(alpha)
+        outputs['vt_dot'] = -vr*vt/r + T*np.cos(alpha)
+
+    def compute_partials(self, inputs, partials, discrete_inputs=None):
+        r = inputs['r']
+        vr = inputs['vr']
+        vt = inputs['vt']
+        alpha = inputs['alpha']
+
+        T = 1e-3
+
+        partials['theta_dot', 'r'] = -vt/r**2
+        partials['theta_dot', 'vt'] = 1/r
+
+        partials['vr_dot', 'r'] = -(vt/r)**2 + 2/r**3
+        partials['vr_dot', 'vt'] = 2*vt/r
+        partials['vr_dot', 'alpha'] = T * np.cos(alpha)
+
+        partials['vt_dot', 'r'] = vt*vr/r**2
+        partials['vt_dot', 'vr'] = -vt/r
+        partials['vt_dot', 'vt'] = -vr/r
+        partials['vt_dot', 'alpha'] = -T * np.sin(alpha)
diff --git a/dymos/examples/moon_landing/__init__.py b/dymos/examples/moon_landing/__init__.py
new file mode 100644
index 00000000..2921dda0
--- /dev/null
+++ b/dymos/examples/moon_landing/__init__.py
@@ -0,0 +1,1 @@
+from .moon_landing_ode import MoonLandingProblemODE
diff --git a/dymos/examples/moon_landing/moon_landing_ode.py b/dymos/examples/moon_landing/moon_landing_ode.py
new file mode 100644
index 00000000..731448ba
--- /dev/null
+++ b/dymos/examples/moon_landing/moon_landing_ode.py
@@ -0,0 +1,45 @@
+import numpy as np
+import openmdao.api as om
+
+
+class MoonLandingProblemODE(om.ExplicitComponent):
+    def initialize(self):
+        self.options.declare('num_nodes', types=int)
+
+    def setup(self):
+        nn = self.options['num_nodes']
+
+        # inputs
+        self.add_input('h', val=np.ones(nn), units=None, desc='Altitude')
+        self.add_input('v', val=np.ones(nn), units='1/s', desc='Velocity')
+        self.add_input('m', val=np.ones(nn), units=None, desc='Mass')
+        self.add_input('T', val=np.ones(nn), units=None, desc='Thrust')
+
+        # outputs
+        self.add_output('h_dot', val=np.ones(nn), units='1/s', desc='Rate of change of Altitude')
+        self.add_output('v_dot', val=np.ones(nn), units='1/s**2', desc='Rate of change of Velocity')
+        self.add_output('m_dot', val=np.ones(nn), units='1/s', desc='Rate of change of Mass')
+
+        # partials
+        ar = np.arange(nn)
+        self.declare_partials(of='h_dot', wrt='v', rows=ar, cols=ar, val=1.0)
+        self.declare_partials(of='v_dot', wrt='m', rows=ar, cols=ar)
+        self.declare_partials(of='v_dot', wrt='T', rows=ar, cols=ar)
+        self.declare_partials(of='m_dot', wrt='T', rows=ar, cols=ar, val=-1/2.349)
+        self.declare_partials(of='m_dot', wrt='T', rows=ar, cols=ar, val=-1/2.349)
+
+    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
+        v = inputs['v']
+        m = inputs['m']
+        T = inputs['T']
+
+        outputs['h_dot'] = v
+        outputs['v_dot'] = -1 + T/m
+        outputs['m_dot'] = -T/2.349
+
+    def compute_partials(self, inputs, partials, discrete_inputs=None):
+        m = inputs['m']
+        T = inputs['T']
+
+        partials['v_dot', 'T'] = 1/m
+        partials['v_dot', 'm'] = -T/m**2
diff --git a/dymos/phase/options.py b/dymos/phase/options.py
index ebc32eb7..d6e723ae 100644
--- a/dymos/phase/options.py
+++ b/dymos/phase/options.py
@@ -496,8 +496,10 @@ class StateOptionsDictionary(om.OptionsDictionary):
                      deprecation="State option 'connected_initial' is deprecated. Please use 'input_initial' instead.")
 
         self.declare(name='input_initial', default=False, types=bool,
-                     desc='Whether an input is created to pass in the initial state. This may be '
-                          'set by a trajectory that links phases.')
+                     desc='Whether the initial value of the state is expected to be connected to an exterior value.')
+
+        self.declare(name='input_final', default=False, types=bool,
+                     desc='Whether the final value of the state is expected to be connected to an exterior value.')
 
 
 class TimeOptionsDictionary(om.OptionsDictionary):
diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index 18bac0e6..73a5f8b1 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -1,4 +1,4 @@
-from collections.abc import Iterable, Callable
+from collections.abc import Iterable, Callable, Sequence
 import inspect
 import warnings
 
@@ -1471,7 +1471,11 @@ class Phase(om.Group):
         rate : bool
             If True, add the rate of change of the named variable to the timeseries outputs of the
             phase.  The rate variable will be named f'{name}_rate'.  Defaults to False.
-        expr :
+        expr : bool
+            True if the given name is an expression for an ExecComp.
+        expr_kwargs : dict
+            Keyword arguments for the expression.
+
 
         Returns
         -------
diff --git a/dymos/phase/simulation_phase.py b/dymos/phase/simulation_phase.py
index f200dc69..ee2c4360 100644
--- a/dymos/phase/simulation_phase.py
+++ b/dymos/phase/simulation_phase.py
@@ -2,7 +2,7 @@ from openmdao.utils.mpi import MPI
 
 from .phase import Phase
 from .._options import options as dymos_options
-from ..transcriptions.grid_data import GaussLobattoGrid, RadauGrid, UniformGrid
+from ..transcriptions.grid_data import BirkhoffGrid, GaussLobattoGrid, RadauGrid, UniformGrid
 from ..transcriptions import ExplicitShooting, GaussLobatto, Radau
 
 from ..utils.misc import _unspecified
@@ -63,7 +63,7 @@ class SimulationPhase(Phase):
             grid = RadauGrid(num_segments=num_seg, nodes_per_seg=seg_order + 1, segment_ends=seg_ends,
                              compressed=compressed)
         elif isinstance(phase_tx.grid_data, GaussLobattoGrid) or \
-                isinstance(phase_tx.grid_data, RadauGrid):
+                isinstance(phase_tx.grid_data, RadauGrid) or isinstance(phase_tx.grid_data, BirkhoffGrid):
             grid = phase_tx.grid_data
         else:
             raise RuntimeError(f'Unexpected grid class for {phase_tx.grid_data}. Only phases with GaussLobatto '
diff --git a/dymos/trajectory/trajectory.py b/dymos/trajectory/trajectory.py
index 6ba9570f..a57cf783 100644
--- a/dymos/trajectory/trajectory.py
+++ b/dymos/trajectory/trajectory.py
@@ -760,7 +760,7 @@ class Trajectory(om.Group):
     def _configure_linkages(self):
 
         if self.options['sim_mode']:
-            # If this is a simulation trajectory, theres no need to link the phases.
+            # If this is a simulation trajectory, there's no need to link the phases.
             return
 
         connected_linkage_inputs = []
diff --git a/dymos/transcriptions/__init__.py b/dymos/transcriptions/__init__.py
index ed0600a7..14426a39 100644
--- a/dymos/transcriptions/__init__.py
+++ b/dymos/transcriptions/__init__.py
@@ -2,4 +2,5 @@ from .analytic.analytic import Analytic
 from .explicit_shooting import ExplicitShooting
 from .pseudospectral.gauss_lobatto import GaussLobatto
 from .pseudospectral.radau_pseudospectral import Radau
+from .pseudospectral.birkhoff import Birkhoff
 from .solve_ivp.solve_ivp import SolveIVP
diff --git a/dymos/transcriptions/explicit_shooting/explicit_shooting.py b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
index 59b6dc8d..1c7d2824 100644
--- a/dymos/transcriptions/explicit_shooting/explicit_shooting.py
+++ b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
@@ -8,7 +8,7 @@ from openmdao.utils.om_warnings import warn_deprecation
 from ..common.timeseries_output_comp import TimeseriesOutputComp
 from .explicit_shooting_continuity_comp import ExplicitShootingContinuityComp
 from ..transcription_base import TranscriptionBase
-from ..grid_data import GaussLobattoGrid, RadauGrid, UniformGrid
+from ..grid_data import BirkhoffGrid, GaussLobattoGrid, RadauGrid, UniformGrid
 from .ode_integration_comp import ODEIntegrationComp
 from ..._options import options as dymos_options
 from ...utils.misc import get_rate_units, CoerceDesvar
@@ -62,10 +62,10 @@ class ExplicitShooting(TranscriptionBase):
                                   'setting this option to False should result in faster execution.')
         self.options.declare('subprob_reports', default=False,
                              desc='Controls the reports made when running the subproblems for ExplicitShooting')
-        self.options.declare('grid', types=(GaussLobattoGrid, RadauGrid, str), allow_none=True, default=None,
+        self.options.declare('grid', types=(GaussLobattoGrid, RadauGrid, BirkhoffGrid, str), allow_none=True, default=None,
                              desc='The grid distribution used to layout the control inputs and provide the default '
                                   'output nodes.')
-        self.options.declare('output_grid', types=(GaussLobattoGrid, RadauGrid, UniformGrid), allow_none=True,
+        self.options.declare('output_grid', types=(GaussLobattoGrid, RadauGrid, UniformGrid, BirkhoffGrid), allow_none=True,
                              default=None,
                              desc='The grid distribution determining the location of the output nodes. The default '
                                   'value of None will result in the use of the grid for outputs. This is useful '
diff --git a/dymos/transcriptions/grid_data.py b/dymos/transcriptions/grid_data.py
index 05403ff7..f35ff6e0 100644
--- a/dymos/transcriptions/grid_data.py
+++ b/dymos/transcriptions/grid_data.py
@@ -1,3 +1,4 @@
+import functools
 import numpy as np
 
 from scipy.linalg import block_diag
@@ -5,6 +6,7 @@ import scipy.sparse as sp
 
 from dymos.utils.lgl import lgl
 from dymos.utils.lgr import lgr
+from dymos.utils.cgl import cgl
 from dymos.utils.hermite import hermite_matrices
 from dymos.utils.lagrange import lagrange_matrices
 
@@ -109,6 +111,68 @@ def radau_pseudospectral_subsets_and_nodes(n, seg_idx, compressed=False):
     return subsets, lgr(n, include_endpoint=True)[0]
 
 
+def birkhoff_subsets_and_nodes(n, grid, seg_idx, compressed=False):
+    """
+    Provides node information and the location of the nodes for n Radau nodes on the range [-1, 1].
+
+    Parameters
+    ----------
+    n : int
+        The total number of nodes in the segment.
+    grid : str
+        The type of Gaussian grid used in the transcription.
+    seg_idx : int
+        The index of this segment within its phase.
+    compressed : bool
+        True if the subset requested is for a phase with compressed transcription.
+
+    Returns
+    -------
+    dict
+        A dictionary with the following keys:
+        'state_disc' gives the indices of the state discretization nodes.
+        'state_input' gives the indices of the state input nodes.
+        'control_disc' gives the indices of the control discretization nodes.
+        'control_input' gives the indices of the control input nodes.
+        'segment_ends' gives the indices of the nodes at the start (even) and end (odd) of a segment.
+        'col' gives the indices of the collocation nodes.
+        'all' gives all node indices.
+    np.array
+        The location of all nodes on [-1, 1].
+
+    Notes
+    -----
+    Subset 'state_input' is the same as subset 'state_disc' if `compressed == False` or
+    `first_seg == True`.  For Radau-Pseudospectral transcription, subset 'control_input' is always
+    the same as subset 'control_disc'.
+    """
+    acceptable_grids = {'lgl', 'lgr', 'cgl'}
+
+    subsets = {
+        'state_disc': np.arange(n, dtype=int),
+        'state_input': np.arange(n, dtype=int),
+        'control_disc': np.arange(n, dtype=int),
+        'control_input': np.arange(n, dtype=int),
+        'segment_ends': np.array([0, n], dtype=int),
+        'col': np.arange(n, dtype=int),
+        'solution': np.arange(n, dtype=int),
+    }
+
+    if grid == 'lgl':
+        nodes = lgl(n)[0]
+        subsets['all'] = np.arange(n, dtype=int)
+    elif grid == 'lgr':
+        nodes = lgr(n, include_endpoint=False)[0]
+        subsets['all'] = np.arange(n, dtype=int)
+    elif grid == 'cgl':
+        nodes = cgl(n)[0]
+        subsets['all'] = np.arange(n, dtype=int)
+    else:
+        raise ValueError(f'Unrecognized grid. Acceptable values are one of {acceptable_grids}')
+
+    return subsets, nodes
+
+
 def uniform_subsets_and_nodes(n, *args, **kwargs):
     """
     Provides a dict of node info and locations for a uniformly distributed set of n nodes on the range [-1, 1].
@@ -304,6 +368,8 @@ class GridData(object):
             self.transcription = 'radau-ps'
         elif transcription.lower() in ['gausslobatto', 'gauss-lobatto', 'lgl']:
             self.transcription = 'gauss-lobatto'
+        elif transcription.lower() in ['birkhoff']:
+            self.transcription = 'birkhoff'
         elif transcription.lower() in ['uniform']:
             self.transcription = 'uniform'
         else:
@@ -316,6 +382,8 @@ class GridData(object):
             get_subsets_and_nodes = radau_pseudospectral_subsets_and_nodes
         elif self.transcription == 'uniform':
             get_subsets_and_nodes = uniform_subsets_and_nodes
+        elif self.transcription == 'birkhoff':
+            get_subsets_and_nodes = functools.partial(birkhoff_subsets_and_nodes, grid=self.grid_type)
 
         # Make sure transcription_order is a vector
         if isinstance(transcription_order, str):
@@ -610,6 +678,34 @@ class GaussLobattoGrid(GridData):
                          segment_ends=segment_ends, compressed=compressed)
 
 
+class BirkhoffGrid(GridData):
+    """
+    A GridData object that provides the node information for the Birkhoff transcription.
+
+    Parameters
+    ----------
+    num_segments : int
+        The number of segments in the phase.
+    nodes_per_seg : int or iterable
+        The number of nodes in each segment. As an integer, it applies to each segment. If a sequence, its length
+        must be equal to num_segments.
+    segment_ends : Iterable[num_segments + 1] or None
+        The segments nodes on some arbitrary interval.
+        This will be normalized to the interval [-1, 1].
+    compressed : bool
+        If the transcription is compressed, then states and controls at shared
+        nodes of adjacent segments are only specified once, and then broadcast
+        to the appropriate indices.
+    grid_type : str
+        The type of Gaussian grid used for the transcription. May be 'lgl' or 'cgl'.
+    """
+    def __init__(self, num_segments, nodes_per_seg, segment_ends=None, compressed=False, grid_type='lgl'):
+        self.grid_type = grid_type
+        super().__init__(num_segments=num_segments, transcription='birkhoff',
+                         transcription_order=np.asarray(nodes_per_seg, dtype=int),
+                         segment_ends=segment_ends, compressed=compressed)
+
+
 class RadauGrid(GridData):
     """
     A GridData object that provides the node information for a Radau distribution.
diff --git a/dymos/transcriptions/pseudospectral/birkhoff.py b/dymos/transcriptions/pseudospectral/birkhoff.py
new file mode 100644
index 00000000..86d7f401
--- /dev/null
+++ b/dymos/transcriptions/pseudospectral/birkhoff.py
@@ -0,0 +1,871 @@
+import numpy as np
+
+import openmdao.api as om
+
+from ..transcription_base import TranscriptionBase
+from ..common import TimeComp, TimeseriesOutputGroup, TimeseriesOutputComp
+from .components import BirkhoffIterGroup, BirkhoffBoundaryGroup
+
+from ..grid_data import BirkhoffGrid
+from dymos.utils.misc import get_rate_units, reshape_val
+from dymos.utils.introspection import get_promoted_vars, get_source_metadata, get_targets
+from dymos.utils.indexing import get_constraint_flat_idxs, get_src_indices_by_row
+
+
+class Birkhoff(TranscriptionBase):
+    """
+        Birkhoff Pseudospectral Transcription.
+
+        Parameters
+        ----------
+        **kwargs : dict
+            Dictionary of optional arguments.
+
+        References
+        ----------
+        I. M. Ross, "A Universeal Birkhoff Theory for Fast Trajectory Optimization"
+        https://arxiv.org/abs/2308.01400v2
+        """
+    def __init__(self, **kwargs):
+        super(Birkhoff, self).__init__(**kwargs)
+        self._rhs_source = 'ode_iter_group.ode_all'
+
+    def initialize(self):
+        """
+        Declare transcription options.
+        """
+        self.options.declare('grid', types=(BirkhoffGrid, str),
+                             allow_none=True, default=None,
+                             desc='The grid distribution used to layout the control inputs and provide the default '
+                                  'output nodes.')
+
+    def init_grid(self):
+        """
+        Setup the GridData object for the Transcription.
+        """
+        if self.options['grid'] in ('gauss-lobatto', None):
+            self.grid_data = BirkhoffGrid(num_segments=self.options['num_segments'],
+                                          nodes_per_seg=self.options['order'],
+                                          segment_ends=self.options['segment_ends'],
+                                          compressed=self.options['compressed'])
+        else:
+            self.grid_data = self.options['grid']
+
+    def setup_time(self, phase):
+        """
+        Setup the time component.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        grid_data = self.grid_data
+
+        super(Birkhoff, self).setup_time(phase)
+
+        time_comp = TimeComp(num_nodes=grid_data.num_nodes, node_ptau=grid_data.node_ptau,
+                             node_dptau_dstau=grid_data.node_dptau_dstau,
+                             units=phase.time_options['units'],
+                             initial_val=phase.time_options['initial_val'],
+                             duration_val=phase.time_options['duration_val'])
+
+        phase.add_subsystem('time', time_comp, promotes_inputs=['*'], promotes_outputs=['*'])
+
+    def configure_time(self, phase):
+        """
+        Configure the inputs/outputs on the time component.
+
+        This method assumes that target introspection has already been performed by the phase and thus
+        options['targets'], options['time_phase_targets'], options['t_initial_targets'],
+        and options['t_duration_targets'] are all correctly populated.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        super(Birkhoff, self).configure_time(phase)
+        phase.time.configure_io()
+        options = phase.time_options
+        ode = phase._get_subsystem(self._rhs_source)
+        ode_inputs = get_promoted_vars(ode, iotypes='input')
+
+        # The tuples here are (name, user_specified_targets, dynamic)
+        for name, targets in [('t', options['targets']),
+                              ('t_phase', options['time_phase_targets'])]:
+            if targets:
+                src_idxs = self.grid_data.subset_node_indices['all']
+                phase.connect(name, [f'ode_all.{t}' for t in targets], src_indices=src_idxs,
+                              flat_src_indices=True)
+                src_idxs = om.slicer[[0, -1], ...]
+                phase.connect(name, [f'boundary_vals.{t}' for t in targets], src_indices=src_idxs)
+
+        for name, targets in [('t_initial', options['t_initial_targets']),
+                              ('t_duration', options['t_duration_targets'])]:
+            for t in targets:
+                shape = ode_inputs[t]['shape']
+
+                if shape == (1,):
+                    src_idxs = endpoint_src_idxs = None
+                    flat_src_idxs = None
+                    src_shape = None
+                else:
+                    src_idxs = np.zeros(self.grid_data.subset_num_nodes['all'], dtype=int)
+                    endpoint_src_idxs = np.zeros(2, dtype=int)
+                    flat_src_idxs = True
+                    src_shape = (1,)
+
+                phase.promotes('ode_all', inputs=[(t, name)], src_indices=src_idxs,
+                               flat_src_indices=flat_src_idxs, src_shape=src_shape)
+
+                phase.promotes('boundary_vals', inputs=[(t, name)], src_indices=endpoint_src_idxs,
+                               flat_src_indices=flat_src_idxs, src_shape=src_shape)
+            if targets:
+                phase.set_input_defaults(name=name,
+                                         val=np.ones((1,)),
+                                         units=options['units'])
+
+    def setup_states(self, phase):
+        """
+        Setup the states for this transcription.
+
+        In the Birkhoff transcription, everything typically done in this
+        method is instead done by the BirkhoffIterGroup.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        pass
+
+    def configure_controls(self, phase):
+        """
+        Configure the inputs/outputs for the controls.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        super().configure_controls(phase)
+
+        if phase.control_options:
+            phase.control_group.configure_io()
+            phase.promotes('control_group',
+                           any=['controls:*', 'control_values:*', 'control_rates:*'])
+
+            phase.connect('dt_dstau', 'control_group.dt_dstau')
+
+        for name, options in phase.control_options.items():
+            if options['targets']:
+                phase.connect(f'control_values:{name}', [f'ode_all.{t}' for t in options['targets']])
+                phase.connect(f'control_values:{name}', [f'boundary_vals.{t}' for t in options['targets']],
+                              src_indices=om.slicer[[0, -1], ...])
+
+            if options['rate_targets']:
+                phase.connect(f'control_rates:{name}_rate',
+                              [f'ode_all.{t}' for t in options['rate_targets']])
+                phase.connect(f'control_rates:{name}_rate',
+                              [f'boundary_vals.{t}' for t in options['rate_targets']],
+                              src_indices=om.slicer[[0, -1], ...])
+
+            if options['rate2_targets']:
+                phase.connect(f'control_rates:{name}_rate2',
+                              [f'ode_all.{t}' for t in options['rate2_targets']])
+                phase.connect(f'control_rates:{name}_rate2',
+                              [f'boundary_vals.{t}' for t in options['rate2_targets']],
+                              src_indices=om.slicer[[0, -1], ...])
+
+    def configure_polynomial_controls(self, phase):
+        """
+        Configure the inputs/outputs for the polynomial controls.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        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],
+                              src_indices=om.slicer[[0, -1], ...])
+
+            targets = get_targets(ode=phase.ode_all, name=f'{name}_rate',
+                                  user_targets=options['rate_targets'])
+            if targets:
+                phase.connect(f'polynomial_control_rates:{name}_rate',
+                              [f'ode_all.{t}' for t in targets])
+                phase.connect(f'polynomial_control_rates:{name}_rate',
+                              [f'boundary_vals.{t}' for t in 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:
+                phase.connect(f'polynomial_control_rates:{name}_rate2',
+                              [f'ode_all.{t}' for t in targets])
+                phase.connect(f'polynomial_control_rates:{name}_rate2',
+                              [f'boundary_vals.{t}' for t in targets],
+                              src_indices=om.slicer[[0, -1], ...])
+
+    def setup_ode(self, phase):
+        """
+        Setup the ode for this transcription.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+
+        ODEClass = phase.options['ode_class']
+        grid_data = self.grid_data
+
+        ode_init_kwargs = phase.options['ode_init_kwargs']
+        ibcs = phase._initial_boundary_constraints
+        fbcs = phase._final_boundary_constraints
+
+        phase.add_subsystem('ode_iter_group',
+                            subsys=BirkhoffIterGroup(grid_data=grid_data, state_options=phase.state_options,
+                                                     time_options=phase.time_options,
+                                                     ode_class=ODEClass,
+                                                     ode_init_kwargs=ode_init_kwargs),
+                            promotes=['*'])
+
+        phase.add_subsystem('boundary_vals',
+                            subsys=BirkhoffBoundaryGroup(grid_data=grid_data,
+                                                         ode_class=ODEClass,
+                                                         ode_init_kwargs=ode_init_kwargs,
+                                                         initial_boundary_constraints=ibcs,
+                                                         final_boundary_constraints=fbcs,
+                                                         objectives=phase._objectives))
+
+    def configure_ode(self, phase):
+        """
+        Create connections to the introspected states.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        grid_data = self.grid_data
+        nn = grid_data.subset_num_nodes['all']
+        map_input_indices_to_disc = grid_data.input_maps['state_input_to_disc']
+        ode = phase._get_subsystem(self._rhs_source)
+        ode_inputs = get_promoted_vars(ode, 'input')
+
+        phase._get_subsystem('boundary_vals').configure_io(phase)
+        phase._get_subsystem('ode_iter_group').configure_io(phase)
+
+    def setup_defects(self, phase):
+        """
+        Create the continuity_comp to house the defects.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        pass
+
+    def configure_defects(self, phase):
+        """
+        Configure the continuity_comp and connect the collocation constraints.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        grid_data = self.grid_data
+
+        for name, options in phase.state_options.items():
+            rate_source_type = phase.classify_var(options['rate_source'])
+            rate_src_path = self._get_rate_source_path(name, phase)
+            if rate_source_type not in ('state', 'ode'):
+                phase.connect(rate_src_path, f'f_computed:{name}')
+
+    def setup_duration_balance(self, phase):
+        """
+        Setup the implicit computation of the phase duration.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        pass
+
+    def configure_duration_balance(self, phase):
+        """
+        Configure the implicit computation of the phase duration.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        pass
+
+    def setup_solvers(self, phase):
+        """
+        Setup the solvers.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        pass
+
+    def configure_solvers(self, phase):
+        """
+        Configure the solvers.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        pass
+
+    def configure_timeseries_outputs(self, phase):
+        """
+        Create connections from time series to all post-introspection sources.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        for timeseries_name, timeseries_options in phase._timeseries.items():
+            timeseries_comp = phase._get_subsystem(f'{timeseries_name}.timeseries_comp')
+            ts_inputs_to_promote = []
+            for ts_output_name, ts_output in timeseries_options['outputs'].items():
+                name = ts_output['output_name'] if ts_output['output_name'] is not None else ts_output['name']
+                units = ts_output['units']
+                shape = ts_output['shape']
+                src = ts_output['src']
+                is_rate = ts_output['is_rate']
+
+                added_src = timeseries_comp._add_output_configure(name,
+                                                                  shape=shape,
+                                                                  units=units,
+                                                                  desc='',
+                                                                  src=src,
+                                                                  rate=is_rate)
+
+                if src.startswith('states:'):
+                    var_name = ts_output["output_name"]
+                    state_name = var_name.split(':')[-1]
+                    ts_inputs_to_promote.append((f'input_values:{var_name}',
+                                                 f'states:{state_name}'))
+                elif added_src:
+                    phase.connect(src_name=src,
+                                  tgt_name=f'{timeseries_name}.input_values:{name}',
+                                  src_indices=ts_output['src_idxs'])
+
+            # print(ts_inputs_to_promote)
+            phase.promotes(timeseries_name, inputs=ts_inputs_to_promote)
+
+    def setup_timeseries_outputs(self, phase):
+        """
+        Setup the timeseries for this transcription.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        gd = self.grid_data
+
+        for name, options in phase._timeseries.items():
+            has_expr = False
+            for _, output_options in options['outputs'].items():
+                if output_options['is_expr']:
+                    has_expr = True
+                    break
+
+            if options['transcription'] is None:
+                ogd = None
+            else:
+                ogd = options['transcription'].grid_data
+
+            timeseries_comp = TimeseriesOutputComp(input_grid_data=gd,
+                                                   output_grid_data=ogd,
+                                                   output_subset=options['subset'],
+                                                   time_units=phase.time_options['units'])
+            timeseries_group = TimeseriesOutputGroup(has_expr=has_expr, timeseries_output_comp=timeseries_comp)
+            phase.add_subsystem(name, subsys=timeseries_group)
+
+            phase.connect('dt_dstau', f'{name}.dt_dstau', flat_src_indices=True)
+
+    def _get_constraint_kwargs(self, constraint_type, options, phase):
+        """
+        Given the constraint options provide the keyword arguments for the OpenMDAO add_constraint method.
+
+        Parameters
+        ----------
+        constraint_type : str
+            One of 'initial', 'final', or 'path'.
+        options : dict
+            The constraint options.
+        phase : Phase
+            The dymos phase to which the constraint applies.
+
+        Returns
+        -------
+        con_output : str
+            The phase-relative path being constrained.
+        constraint_kwargs : dict
+            Keyword arguments for the OpenMDAO add_constraint method.
+        """
+        num_nodes = self._get_num_timeseries_nodes()
+
+        constraint_kwargs = {key: options for key, options in options.items()}
+        con_name = constraint_kwargs.pop('constraint_name')
+
+        # Determine the path to the variable which we will be constraining
+        var = con_name if options['is_expr'] else options['name']
+        var_type = phase.classify_var(var)
+
+        # These are the flat indices at a single point in time used
+        # in either initial, final, or path constraints.
+        idxs_in_initial = phase._indices_in_constraints(var, 'initial')
+        idxs_in_final = phase._indices_in_constraints(var, 'final')
+        idxs_in_path = phase._indices_in_constraints(var, 'path')
+
+        size = np.prod(options['shape'], dtype=int)
+
+        flat_idxs = get_constraint_flat_idxs(options)
+
+        # Now we need to convert the indices given by the user at any given point
+        # to flat indices to be given to OpenMDAO as flat indices spanning the phase.
+        if var_type == 'parameter':
+            if any([idxs_in_initial.intersection(idxs_in_final),
+                    idxs_in_initial.intersection(idxs_in_path),
+                    idxs_in_final.intersection(idxs_in_path)]):
+                raise RuntimeError(f'In phase {phase.pathname}, parameter `{var}` is subject to multiple boundary '
+                                   f'or path constraints.\nParameters are single values that do not change in '
+                                   f'time, and may only be used in a single boundary or path constraint.')
+            constraint_kwargs['indices'] = flat_idxs
+        else:
+            if constraint_type == 'initial':
+                constraint_kwargs['indices'] = flat_idxs
+            elif constraint_type == 'final':
+                constraint_kwargs['indices'] = flat_idxs
+            else:
+                # This is a path constraint.
+                # Remove any flat indices involved in an initial constraint from the path constraint
+                flat_idxs_set = set(flat_idxs)
+                idxs_not_in_initial = list(flat_idxs_set - idxs_in_initial)
+
+                # Remove any flat indices involved in the final constraint from the path constraint
+                idxs_not_in_final = list(flat_idxs_set - idxs_in_final)
+                idxs_not_in_final = (size * (num_nodes - 1) + np.asarray(idxs_not_in_final)).tolist()
+
+                intermediate_idxs = []
+                for i in range(1, num_nodes - 1):
+                    intermediate_idxs.extend(size * i + flat_idxs)
+
+                constraint_kwargs['indices'] = idxs_not_in_initial + intermediate_idxs + idxs_not_in_final
+
+        alias_map = {'path': 'path_constraint',
+                     'initial': 'initial_boundary_constraint',
+                     'final': 'final_boundary_constraint'}
+
+        str_idxs = '' if options['indices'] is None else f'{options["indices"]}'
+
+        constraint_kwargs['alias'] = f'{phase.pathname}->{alias_map[constraint_type]}->{con_name}{str_idxs}'
+        constraint_kwargs.pop('name')
+        con_path = constraint_kwargs.pop('constraint_path')
+        constraint_kwargs.pop('shape')
+        constraint_kwargs['flat_indices'] = True
+        constraint_kwargs.pop('is_expr')
+
+        return con_path, constraint_kwargs
+
+    def _get_objective_src(self, var, loc, phase, ode_outputs=None):
+        """
+        Return the path to the variable that will be used as the objective.
+
+        Parameters
+        ----------
+        var : str
+            Name of the variable to be used as the objective.
+        loc : str
+            The location of the objective in the phase ['initial', 'final'].
+        phase : dymos.Phase
+            Phase object containing in which the objective resides.
+        ode_outputs : dict or None
+            A dictionary of ODE outputs as returned by get_promoted_vars.
+
+        Returns
+        -------
+        obj_path : str
+            Path to the source.
+        shape : tuple
+            Source shape.
+        units : str
+            Source units.
+        linear : bool
+            True if the objective quantity1 is linear.
+        """
+        time_units = phase.time_options['units']
+        var_type = phase.classify_var(var)
+
+        if ode_outputs is None:
+            ode_outputs = get_promoted_vars(phase._get_subsystem(self._rhs_source), 'output')
+
+        if var_type == 't':
+            shape = (1,)
+            units = time_units
+            linear = True
+            constraint_path = 't'
+        elif var_type == 't_phase':
+            shape = (1,)
+            units = time_units
+            linear = True
+            constraint_path = 't_phase'
+        elif var_type == 'state':
+            shape = phase.state_options[var]['shape']
+            units = phase.state_options[var]['units']
+            solve_segments = phase.state_options[var]['solve_segments']
+            linear = False
+            constraint_path = f'{loc}_states:{var}'
+        elif var_type == 'indep_control':
+            shape = phase.control_options[var]['shape']
+            units = phase.control_options[var]['units']
+            linear = True
+            constraint_path = f'control_values:{var}'
+        elif var_type == 'input_control':
+            shape = phase.control_options[var]['shape']
+            units = phase.control_options[var]['units']
+            linear = False
+            constraint_path = f'control_values:{var}'
+        elif var_type == 'indep_polynomial_control':
+            shape = phase.polynomial_control_options[var]['shape']
+            units = phase.polynomial_control_options[var]['units']
+            linear = True
+            constraint_path = f'polynomial_control_values:{var}'
+        elif var_type == 'input_polynomial_control':
+            shape = phase.polynomial_control_options[var]['shape']
+            units = phase.polynomial_control_options[var]['units']
+            linear = False
+            constraint_path = f'polynomial_control_values:{var}'
+        elif var_type == 'parameter':
+            shape = phase.parameter_options[var]['shape']
+            units = phase.parameter_options[var]['units']
+            linear = True
+            constraint_path = f'parameter_vals:{var}'
+        elif var_type in ('control_rate', 'control_rate2'):
+            control_var = var[:-5] if var_type == 'control_rate' else var[:-6]
+            shape = phase.control_options[control_var]['shape']
+            control_units = phase.control_options[control_var]['units']
+            d = 2 if var_type == 'control_rate2' else 1
+            control_rate_units = get_rate_units(control_units, time_units, deriv=d)
+            units = control_rate_units
+            linear = False
+            constraint_path = f'control_rates:{var}'
+        elif var_type in ('polynomial_control_rate', 'polynomial_control_rate2'):
+            control_var = var[:-5]
+            shape = phase.polynomial_control_options[control_var]['shape']
+            control_units = phase.polynomial_control_options[control_var]['units']
+            d = 2 if var_type == 'polynomial_control_rate2' else 1
+            control_rate_units = get_rate_units(control_units, time_units, deriv=d)
+            units = control_rate_units
+            linear = False
+            constraint_path = f'polynomial_control_rates:{var}'
+        elif var_type == 'timeseries_exec_comp_output':
+            shape = (1,)
+            units = None
+            constraint_path = f'timeseries.timeseries_exec_comp.{var}'
+            linear = False
+        else:
+            # Failed to find variable, assume it is in the ODE. This requires introspection.
+            constraint_path = f'boundary_val.{var}'
+            meta = get_source_metadata(ode_outputs, var, user_units=None, user_shape=None)
+            shape = meta['shape']
+            units = meta['units']
+            linear = False
+
+        return constraint_path, shape, units, linear
+
+    def _get_num_timeseries_nodes(self):
+        """
+        Returns the number of nodes in the default timeseries for this transcription.
+
+        Returns
+        -------
+        int
+            The number of nodes in the default timeseries for this transcription.
+        """
+        return self.grid_data.num_nodes
+
+    def _get_rate_source_path(self, state_name, phase):
+        """
+        Return the rate source location and indices for a given state name.
+
+        Parameters
+        ----------
+        state_name : str
+            Name of the state.
+        phase : dymos.Phase
+            Phase object containing the rate source.
+
+        Returns
+        -------
+        str
+            Path to the rate source.
+        ndarray
+            Array of source indices.
+        """
+        gd = self.grid_data
+        try:
+            var = phase.state_options[state_name]['rate_source']
+        except RuntimeError:
+            raise ValueError(f"state '{state_name}' in phase '{phase.name}' was not given a "
+                             "rate_source")
+
+        # Note the rate source must be shape-compatible with the state
+        var_type = phase.classify_var(var)
+
+        # Determine the path to the variable
+        if var_type == 't':
+            rate_path = 't'
+        elif var_type == 't_phase':
+            rate_path = 't_phase'
+        elif var_type == 'state':
+            rate_path = f'states:{var}'
+        elif var_type == 'indep_control':
+            rate_path = f'control_values:{var}'
+        elif var_type == 'input_control':
+            rate_path = f'control_values:{var}'
+        elif var_type == 'control_rate':
+            control_name = var[:-5]
+            rate_path = f'control_rates:{control_name}_rate'
+        elif var_type == 'control_rate2':
+            control_name = var[:-6]
+            rate_path = f'control_rates:{control_name}_rate2'
+        elif var_type == 'indep_polynomial_control':
+            rate_path = f'polynomial_control_values:{var}'
+        elif var_type == 'input_polynomial_control':
+            rate_path = f'polynomial_control_values:{var}'
+        elif var_type == 'polynomial_control_rate':
+            control_name = var[:-5]
+            rate_path = f'polynomial_control_rates:{control_name}_rate'
+        elif var_type == 'polynomial_control_rate2':
+            control_name = var[:-6]
+            rate_path = f'polynomial_control_rates:{control_name}_rate2'
+        elif var_type == 'parameter':
+            rate_path = f'parameter_vals:{var}'
+            dynamic = not phase.parameter_options[var]['static_target']
+            if dynamic:
+                node_idxs = np.zeros(gd.subset_num_nodes['col'], dtype=int)
+            else:
+                node_idxs = np.zeros(1, dtype=int)
+        else:
+            # Failed to find variable, assume it is in the ODE
+            rate_path = f'ode_all.{var}'
+
+        return rate_path
+
+    def _get_timeseries_var_source(self, var, output_name, phase):
+        """
+        Return the source path and indices for a given variable to be connected to a timeseries.
+
+        Parameters
+        ----------
+        var : str
+            Name of the timeseries variable whose source is desired.
+        output_name : str
+            Name of the timeseries output whose source is desired.
+        phase : dymos.Phase
+            Phase object containing the variable, either as state, time, control, etc., or as an ODE output.
+
+        Returns
+        -------
+        meta : dict
+            Metadata pertaining to the variable at the given path. This dict contains 'src' (the path to the
+            timeseries source), 'src_idxs' (an array of the
+            source indices), 'units' (the units of the source variable), and 'shape' (the shape of the variable at
+            a given node).
+        """
+        gd = self.grid_data
+        var_type = phase.classify_var(var)
+        time_units = phase.time_options['units']
+
+        transcription = phase.options['transcription']
+        ode = transcription._get_ode(phase)
+        ode_outputs = get_promoted_vars(ode, 'output')
+
+        # The default for node_idxs, applies to everything except states and parameters.
+        node_idxs = gd.subset_node_indices['all']
+
+        meta = {}
+
+        # Determine the path to the variable
+        if var_type == 't':
+            path = 't'
+            src_units = time_units
+            src_shape = (1,)
+        elif var_type == 't_phase':
+            path = 't_phase'
+            src_units = time_units
+            src_shape = (1,)
+        elif var_type == 'state':
+            path = f'states:{var}'
+            src_units = phase.state_options[var]['units']
+            src_shape = phase.state_options[var]['shape']
+
+            if gd.transcription == 'radau-ps':
+                # Find the state_input indices which occur at segment endpoints, and repeat them twice
+                state_input_idxs = gd.subset_node_indices['state_input']
+                repeat_idxs = np.ones_like(state_input_idxs)
+                if self.options['compressed']:
+                    segment_end_idxs = gd.subset_node_indices['segment_ends'][1:-1]
+                    # Repeat nodes that are on segment bounds (but not the first or last nodes in the phase)
+                    nodes_to_repeat = list(set(state_input_idxs).intersection(set(segment_end_idxs)))
+                    # Now find these nodes in the state input indices
+                    idxs_of_ntr_in_state_inputs = np.where(np.in1d(state_input_idxs, nodes_to_repeat))[0]
+                    # All state input nodes are used once, but nodes_to_repeat are used twice
+                    repeat_idxs[idxs_of_ntr_in_state_inputs] = 2
+                # Now we have a way of mapping the state input indices to all nodes
+                map_input_node_idxs_to_all = np.repeat(np.arange(gd.subset_num_nodes['state_input'],
+                                                                 dtype=int), repeats=repeat_idxs)
+                # Now select the subset of nodes we want to use.
+                node_idxs = map_input_node_idxs_to_all[gd.subset_node_indices['all']]
+        elif var_type in ['indep_control', 'input_control']:
+            path = f'control_values:{var}'
+            src_units = phase.control_options[var]['units']
+            src_shape = phase.control_options[var]['shape']
+        elif var_type == 'control_rate':
+            control_name = var[:-5]
+            path = f'control_rates:{control_name}_rate'
+            control_name = var[:-5]
+            src_units = get_rate_units(phase.control_options[control_name]['units'], time_units, deriv=1)
+            src_shape = phase.control_options[control_name]['shape']
+        elif var_type == 'control_rate2':
+            control_name = var[:-6]
+            path = f'control_rates:{control_name}_rate2'
+            src_units = get_rate_units(phase.control_options[control_name]['units'], time_units, deriv=2)
+            src_shape = phase.control_options[control_name]['shape']
+        elif var_type in ['indep_polynomial_control', 'input_polynomial_control']:
+            path = f'polynomial_control_values:{var}'
+            src_units = phase.polynomial_control_options[var]['units']
+            src_shape = phase.polynomial_control_options[var]['shape']
+        elif var_type == 'polynomial_control_rate':
+            control_name = var[:-5]
+            path = f'polynomial_control_rates:{control_name}_rate'
+            control = phase.polynomial_control_options[control_name]
+            src_units = get_rate_units(control['units'], time_units, deriv=1)
+            src_shape = control['shape']
+        elif var_type == 'polynomial_control_rate2':
+            control_name = var[:-6]
+            path = f'polynomial_control_rates:{control_name}_rate2'
+            control = phase.polynomial_control_options[control_name]
+            src_units = get_rate_units(control['units'], time_units, deriv=2)
+            src_shape = control['shape']
+        elif var_type == 'parameter':
+            path = f'parameter_vals:{var}'
+            # Timeseries are never a static_target
+            node_idxs = np.zeros(gd.subset_num_nodes['all'], dtype=int)
+            src_units = phase.parameter_options[var]['units']
+            src_shape = phase.parameter_options[var]['shape']
+        else:
+            # Failed to find variable, assume it is in the ODE
+            path = f'ode_all.{var}'
+            meta = get_source_metadata(ode_outputs, src=var)
+            src_shape = meta['shape']
+            src_units = meta['units']
+            src_tags = meta['tags']
+            if 'dymos.static_output' in src_tags:
+                raise RuntimeError(
+                    f'ODE output {var} is tagged with "dymos.static_output" and cannot be a timeseries output.')
+
+        src_idxs = om.slicer[node_idxs, ...]
+
+        meta['src'] = path
+        meta['src_idxs'] = src_idxs
+        meta['units'] = src_units
+        meta['shape'] = src_shape
+
+        return meta
+
+    def get_parameter_connections(self, name, phase):
+        """
+        Returns info about a parameter's target connections in the phase.
+
+        Parameters
+        ----------
+        name : str
+            Parameter name.
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+
+        Returns
+        -------
+        list of (paths, indices)
+            A list containing a tuple of target paths and corresponding src_indices to which the
+            given design variable is to be connected.
+        """
+        connection_info = []
+        if name in phase.parameter_options:
+            options = phase.parameter_options[name]
+            for tgt in options['targets']:
+                if tgt in options['static_targets']:
+                    src_idxs = np.squeeze(get_src_indices_by_row([0], options['shape']), axis=0)
+                    endpoint_src_idxs = om.slicer[:, ...]
+                else:
+                    src_idxs_raw = np.zeros(self.grid_data.subset_num_nodes['all'], dtype=int)
+                    src_idxs = get_src_indices_by_row(src_idxs_raw, options['shape'])
+                    endpoint_src_idxs = om.slicer[[0, -1], ...]
+                    if options['shape'] == (1,):
+                        src_idxs = src_idxs.ravel()
+
+                connection_info.append((f'ode_all.{tgt}', (src_idxs,)))
+                connection_info.append((f'boundary_vals.{tgt}', (endpoint_src_idxs,)))
+
+        return connection_info
+
+    def _requires_continuity_constraints(self, phase):
+        """
+        Tests whether state and/or control and/or control rate continuity are required.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase to which this transcription applies.
+
+        Returns
+        -------
+        any_state_continuity : bool
+            True if any state continuity is required to be enforced.
+        any_control_continuity : bool
+            True if any control value continuity is required to be enforced.
+        any_control_rate_continuity : bool
+            True if any control rate continuity is required to be enforced.
+        """
+        num_seg = self.grid_data.num_segments
+        compressed = self.grid_data.compressed
+
+        any_state_continuity = num_seg > 1 and not compressed
+        any_control_continuity = any([opts['continuity'] for opts in phase.control_options.values()])
+        any_control_continuity = any_control_continuity and num_seg > 1
+        any_rate_continuity = any([opts['rate_continuity'] or opts['rate2_continuity']
+                                   for opts in phase.control_options.values()])
+        any_rate_continuity = any_rate_continuity and num_seg > 1
+
+        return any_state_continuity, any_control_continuity, any_rate_continuity
diff --git a/dymos/transcriptions/pseudospectral/components/__init__.py b/dymos/transcriptions/pseudospectral/components/__init__.py
index 7714eb11..9b9e01c0 100644
--- a/dymos/transcriptions/pseudospectral/components/__init__.py
+++ b/dymos/transcriptions/pseudospectral/components/__init__.py
@@ -3,4 +3,6 @@ from .state_interp_comp import StateInterpComp
 from .state_independents import StateIndependentsComp
 from .control_endpoint_defect_comp import ControlEndpointDefectComp
 from .gauss_lobatto_interleave_comp import GaussLobattoInterleaveComp
-from ...common.timeseries_output_comp import TimeseriesOutputComp
+from .birkhoff_collocation_comp import BirkhoffCollocationComp
+from .birkhoff_iter_group import BirkhoffIterGroup
+from .birkhoff_boundary_group import BirkhoffBoundaryGroup
diff --git a/dymos/transcriptions/pseudospectral/components/birkhoff_boundary_group.py b/dymos/transcriptions/pseudospectral/components/birkhoff_boundary_group.py
new file mode 100644
index 00000000..44fedc96
--- /dev/null
+++ b/dymos/transcriptions/pseudospectral/components/birkhoff_boundary_group.py
@@ -0,0 +1,165 @@
+import numpy as np
+import openmdao.api as om
+
+from .birkhoff_collocation_comp import BirkhoffCollocationComp
+from .birkhoff_state_resid_comp import BirkhoffStateResidComp
+
+from ...grid_data import GridData
+from ....phase.options import TimeOptionsDictionary
+from dymos._options import options as dymos_options
+
+
+class BirkhoffBoundaryMuxComp(om.ExplicitComponent):
+    """
+    Class definition of the BirtkhoffBoundaryMuxComp.
+
+    This component takes the initial and final values of states
+    and muxes them into a single output.
+
+    For a state of a given `shape` in a phase with `num_seg` segments,
+    the shape of `initial_states:{state_name}` and `final_states:{state_name}`
+    are both `(num_seg,) + shape` and the shape of the resulting
+    `states:{state_name}` is `(2,) + shape`.
+
+    Parameters
+    ----------
+    **kwargs : dict
+        Dictionary of optional phase arguments.
+    """
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        self._io_names = {}
+        self._no_check_partials = not dymos_options['include_check_partials']
+
+    def configure_io(self, state_options):
+        """
+        I/O creation is delayed until configure so that we can determine shape and units for the states.
+
+        Parameters
+        ----------
+        state_options : StateOptionsDictionary
+            The phase object to which this transcription instance applies.
+        """
+        self._io_names = {}
+        for state_name, options in state_options.items():
+            shape = options['shape']
+            size = np.prod(shape, dtype=int)
+            units = options['units']
+            self._io_names[state_name] = {'initial': f'initial_states:{state_name}',
+                                          'final': f'final_states:{state_name}',
+                                          'boundary': state_name}
+
+            iname = self._io_names[state_name]['initial']
+            fname = self._io_names[state_name]['final']
+            bname = self._io_names[state_name]['boundary']
+
+            self.add_input(iname, shape=shape, units=units)
+            self.add_input(fname, shape=shape, units=units)
+            self.add_output(bname, shape=(2,) + shape, units=units)
+
+            ar = np.arange(size, dtype=int)
+
+            self.declare_partials(of=bname, wrt=iname, rows=ar, cols=ar, val=1.0)
+            self.declare_partials(of=bname, wrt=fname, rows=ar + size, cols=ar, val=1.0)
+
+    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
+        """
+        Compute component outputs.
+
+        Parameters
+        ----------
+        inputs : Vector
+            Unscaled, dimensional input variables read via inputs[key].
+        outputs : Vector
+            Unscaled, dimensional output variables read via outputs[key].
+        discrete_inputs : dict or None
+            If not None, dict containing discrete input values.
+        discrete_outputs : dict or None
+            If not None, dict containing discrete output values.
+        """
+        for state_name, io_names in self._io_names.items():
+            outputs[io_names['boundary']][0] = inputs[io_names['initial']]
+            outputs[io_names['boundary']][1] = inputs[io_names['final']]
+
+
+class BirkhoffBoundaryGroup(om.Group):
+    """
+    Class definition for the BirkhoffBoundaryEvalGroup.
+
+    This group accepts values for initial and final times, states, controls, and parameters
+    and evaluates the ODE with those in order to compute the boundary values and
+    objectives.
+
+    Note that in the Birkhoff transcription, the initial and final state values are
+    decoupled from the initial and final states in the interpolating polynomial.
+
+    Dymos uses the Birkhoff LGL or CGL approaches so that the control values are provided
+    at the endpoints of the phase without the need for extrapolation (unlike the classical
+    Radau approach in Dymos)
+
+    Parameters
+    ----------
+    **kwargs : dict
+        Dictionary of optional arguments.
+    """
+
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        self._no_check_partials = not dymos_options['include_check_partials']
+
+    def initialize(self):
+        """
+        Declare group options.
+        """
+        self.options.declare('grid_data', types=GridData, desc='Container object for grid info.')
+        self.options.declare('ode_class', default=None,
+                             desc='Callable that instantiates the ODE system.',
+                             recordable=False)
+        self.options.declare('ode_init_kwargs', types=dict, default={}, recordable=False,
+                             desc='Keyword arguments provided when initializing the ODE System')
+        self.options.declare('initial_boundary_constraints', types=list, recordable=False,
+                             desc='Initial boundary constraints from the containing phase.')
+        self.options.declare('final_boundary_constraints', types=list, recordable=False,
+                             desc='Final boundary constraints from the containing phase.')
+        self.options.declare('objectives', types=dict, recordable=False,
+                             desc='Objectives from the containing phase.')
+
+    def setup(self):
+        """
+        Define the structure of the control group.
+        """
+        ode_class = self.options['ode_class']
+        ode_init_kwargs = self.options['ode_init_kwargs']
+        ibcs = self.options['initial_boundary_constraints']
+        fbcs = self.options['final_boundary_constraints']
+        objs = [meta for meta in self.options['objectives'].values()]
+
+        self.add_subsystem('boundary_mux', subsys=BirkhoffBoundaryMuxComp(),
+                           promotes_inputs=['*'], promotes_outputs=['*'])
+
+        self.add_subsystem('boundary_ode', subsys=ode_class(num_nodes=2, **ode_init_kwargs),
+                           promotes_inputs=['*'], promotes_outputs=['*'])
+
+        if any([response['is_expr'] for response in ibcs + fbcs + objs]):
+            print([response['is_expr'] for response in ibcs + fbcs + objs])
+            self.add_subsystem('boundary_constraint_exec_comp', subsys=om.ExecComp(),
+                               promotes_inputs=['*'], promotes_outputs=['*'])
+
+    def configure_io(self, phase):
+        """
+        I/O creation is delayed until configure so that we can determine shape and units for the states.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        self._get_subsystem('boundary_mux').configure_io(state_options=phase.state_options)
+
+        for state_name, options in phase.state_options.items():
+            for tgt in options['targets']:
+                self.promotes('boundary_ode', inputs=[(tgt, state_name)])
+            if options['targets']:
+                self.set_input_defaults(name=state_name,
+                                        val=options['val'],
+                                        units=options['units'])
diff --git a/dymos/transcriptions/pseudospectral/components/birkhoff_collocation_comp.py b/dymos/transcriptions/pseudospectral/components/birkhoff_collocation_comp.py
new file mode 100644
index 00000000..39780899
--- /dev/null
+++ b/dymos/transcriptions/pseudospectral/components/birkhoff_collocation_comp.py
@@ -0,0 +1,415 @@
+import numpy as np
+import scipy
+import scipy.sparse as sp
+
+import openmdao.api as om
+
+from dymos.transcriptions.grid_data import GridData
+from dymos.utils.misc import get_rate_units
+from dymos._options import options as dymos_options
+from dymos.utils.lgl import lgl
+from dymos.utils.cgl import cgl
+from dymos.utils.birkhoff import birkhoff_matrix
+
+
+class BirkhoffCollocationComp(om.ExplicitComponent):
+    """
+    Class definition for the BirkhoffCollocationComp.
+
+    BirkhoffCollocationComp computes the generalized defects of a segment for implicit collocation.
+    There are four defects to be evaluated; state, state rate, initial state, and final state.
+    The state defect is the difference between the state value design variables and the state rate
+    design variables.
+    The state rate defect is the difference between the state rate design variables and the computed
+    state derivative at the collocation nodes.
+    The initial and final state defects are the differences between the initial and final state
+    design variables and the first and last index of the state design variable array.
+
+    Parameters
+    ----------
+    **kwargs : dict
+        Dictionary of optional arguments.
+    """
+
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        self._no_check_partials = not dymos_options['include_check_partials']
+
+    def initialize(self):
+        """
+        Declare component options.
+        """
+        self.options.declare(
+            'grid_data', types=GridData,
+            desc='Container object for grid info')
+
+        self.options.declare(
+            'state_options', types=dict,
+            desc='Dictionary of state names/options for the phase')
+
+        self.options.declare(
+            'time_units', default=None, allow_none=True, types=str,
+            desc='Units of time')
+
+    def configure_io(self, phase):
+        """
+        I/O creation is delayed until configure so we can determine shape and units.
+
+        Parameters
+        ----------
+        phase : Phase
+            The phase object that contains this collocation comp.
+        """
+        gd = self.options['grid_data']
+        num_nodes = gd.subset_num_nodes['col']
+        num_segs = gd.num_segments
+        time_units = self.options['time_units']
+        state_options = self.options['state_options']
+
+        self.add_input('dt_dstau', units=time_units, shape=(gd.subset_num_nodes['col'],))
+        self.var_names = var_names = {}
+        for state_name, options in state_options.items():
+            var_names[state_name] = {
+                'f_value': f'state_rates:{state_name}',
+                'f_computed': f'f_computed:{state_name}',
+                'state_value': f'states:{state_name}',
+                'state_segment_ends': f'state_segment_ends:{state_name}',
+                'state_initial_value': f'initial_states:{state_name}',
+                'state_final_value': f'final_states:{state_name}',
+                'state_defect': f'state_defects:{state_name}',
+                'state_rate_defect': f'state_rate_defects:{state_name}',
+                'initial_state_defect': f'initial_state_defects:{state_name}',
+                'final_state_defect': f'final_state_defects:{state_name}',
+                'state_continuity_defect': f'state_cnty_defects:{state_name}'
+            }
+
+        for state_name, options in state_options.items():
+            shape = options['shape']
+            units = options['units']
+
+            rate_units = get_rate_units(units, time_units)
+            var_names = self.var_names[state_name]
+
+            rate_source = options['rate_source']
+            rate_source_type = phase.classify_var(rate_source)
+
+            self.add_input(
+                name=var_names['f_value'],
+                shape=(num_nodes,) + shape,
+                desc=f'Estimated derivative of state {state_name} at the polynomial nodes',
+                units=units)
+
+            if rate_source_type == 'state':
+                var_names['f_computed'] = f'states:{rate_source}'
+            else:
+                self.add_input(
+                    name=var_names['f_computed'],
+                    shape=(num_nodes,) + shape,
+                    desc=f'Computed derivative of state {state_name} at the polynomial nodes',
+                    units=rate_units)
+
+            self.add_input(
+                name=var_names['state_value'],
+                shape=(num_nodes,) + shape,
+                desc=f'Value of the state {state_name} at the polynomial nodes',
+                units=units
+            )
+
+            self.add_input(
+                name=var_names['state_segment_ends'],
+                shape=(num_segs, 2) + shape,
+                desc=f'Initial and final value of the state {state_name} at each segment end.',
+                units=units
+            )
+
+            self.add_input(
+                name=var_names['state_initial_value'],
+                shape=(1,) + shape,
+                desc=f'Desired initial value of state {state_name}',
+                units=units
+            )
+
+            self.add_input(
+                name=var_names['state_final_value'],
+                shape=(1,) + shape,
+                desc=f'Desired final value of state {state_name}',
+                units=units
+            )
+
+            self.add_output(
+                name=var_names['state_defect'],
+                shape=(num_nodes + num_segs,) + shape,
+                units=units
+            )
+
+            self.add_output(
+                name=var_names['state_rate_defect'],
+                shape=(num_nodes,) + shape,
+                units=units
+            )
+
+            self.add_output(
+                name=var_names['initial_state_defect'],
+                shape=shape,
+                units=units
+            )
+
+            self.add_output(
+                name=var_names['final_state_defect'],
+                shape=shape,
+                units=units
+            )
+
+            if 'defect_ref' in options and options['defect_ref'] is not None:
+                defect_ref = options['defect_ref']
+            elif 'defect_scaler' in options and options['defect_scaler'] is not None:
+                defect_ref = np.divide(1.0, options['defect_scaler'])
+            else:
+                if 'ref' in options and options['ref'] is not None:
+                    defect_ref = options['ref']
+                elif 'scaler' in options and options['scaler'] is not None:
+                    defect_ref = np.divide(1.0, options['scaler'])
+                else:
+                    defect_ref = 1.0
+
+            if not np.isscalar(defect_ref):
+                defect_ref = np.asarray(defect_ref)
+                if defect_ref.shape == shape:
+                    defect_ref_state = np.tile(defect_ref.flatten(), num_nodes+num_segs)
+                    defect_ref_v = np.tile(defect_ref.flatten(), num_nodes)
+                else:
+                    raise ValueError('array-valued scaler/ref must length equal to state-size')
+            else:
+                defect_ref_state = defect_ref
+                defect_ref_v = defect_ref
+
+            if not options['solve_segments']:
+                self.add_constraint(name=var_names['state_defect'],
+                                    equals=0.0,
+                                    ref=defect_ref_state)
+
+                self.add_constraint(name=var_names['state_rate_defect'],
+                                    equals=0.0,
+                                    ref=defect_ref_v)
+
+                self.add_constraint(name=var_names['initial_state_defect'],
+                                    equals=0.0,
+                                    ref=defect_ref)
+
+                self.add_constraint(name=var_names['final_state_defect'],
+                                    equals=0.0,
+                                    ref=defect_ref)
+
+                if gd.num_segments > 1:
+                    self.add_constraint(name=var_names['state_continuity_defect'],
+                                        equals=0.0,
+                                        ref=defect_ref)
+
+        A_blocks = []
+        B_blocks = []
+        C_blocks = []
+
+        # _xv_idxs is a set of indices that arranges the stacked
+        # [[X^T],[V^T]] arrays into a segment-by-segment ordering
+        # instead of all of the state values followed by all of the state rates.
+        self._xv_idxs = []
+
+        idx0 = 0
+
+        for i in range(num_segs):
+            nnps_i = gd.subset_num_nodes_per_segment['all'][0]
+            if gd.grid_type == 'lgl':
+                tau_i, w_i = lgl(nnps_i)
+            elif gd.grid_type == 'cgl':
+                tau_i, w_i = cgl(nnps_i)
+            else:
+                raise ValueError('invalid grid type')
+
+            B_i = birkhoff_matrix(tau_i, w_i, grid_type=gd.grid_type)
+
+            A_i = np.zeros((nnps_i + 1, 2 * nnps_i))
+            A_i[:nnps_i, :nnps_i] = np.eye(nnps_i)
+            A_i[:nnps_i, nnps_i:] = -B_i
+            A_i[-1, nnps_i:] = B_i[-1, :]
+
+            C_i = np.zeros((nnps_i + 1, 2))
+            C_i[:-1, 0] = 1.
+            C_i[-1, :] = [-1, 1]
+
+            A_blocks.append(A_i)
+            B_blocks.append(B_i)
+            C_blocks.append(C_i)
+
+            ar_nnps = np.arange(nnps_i, dtype=int)
+            self._xv_idxs.extend(idx0 + ar_nnps)
+            self._xv_idxs.extend(idx0 + ar_nnps + num_nodes)
+            idx0 += nnps_i
+
+        self._A = scipy.linalg.block_diag(*A_blocks)
+        self._B = scipy.linalg.block_diag(*B_blocks)
+        self._C = scipy.linalg.block_diag(*C_blocks)
+
+        # Setup partials
+        for state_name, options in state_options.items():
+            shape = options['shape']
+            size = np.prod(shape)
+
+            ar1 = np.arange(num_nodes * size)
+            c1 = np.repeat(np.arange(num_nodes), size)
+
+            var_names = self.var_names[state_name]
+
+            # The derivative of the state defect wrt x_ab is [-C].
+            # Take the Kronecker product of that and an identity matrix
+            # of the state's size to get the pattern that includes each
+            # individual element in the state.
+            c_sparse = sp.kron(self._C, sp.eye(size), format='csr')
+            c_rows, c_cols = c_sparse.nonzero()
+            c_data = c_sparse.data.ravel()
+
+            self.declare_partials(of=var_names['state_defect'],
+                                  wrt=var_names['state_segment_ends'],
+                                  rows=c_rows, cols=c_cols, val=-c_data)
+
+            d_state_defect_dXV = self._A
+            dXV_dX = np.vstack((np.eye(num_nodes), np.zeros((num_nodes, num_nodes))))[self._xv_idxs]
+            dXV_dV = np.vstack((np.zeros((num_nodes, num_nodes)), np.eye(num_nodes)))[self._xv_idxs]
+            d_state_defect_dX = sp.csr_matrix(np.kron(np.dot(d_state_defect_dXV, dXV_dX), np.eye(size)))
+            d_state_defect_dV = sp.csr_matrix(np.kron(np.dot(d_state_defect_dXV, dXV_dV), np.eye(size)))
+
+            rs_dX, cs_dX = d_state_defect_dX.nonzero()
+            rs_dV, cs_dV = d_state_defect_dV.nonzero()
+
+            self.declare_partials(of=var_names['state_defect'],
+                                  wrt=var_names['state_value'],
+                                  rows=rs_dX, cols=cs_dX, val=d_state_defect_dX.data.ravel())
+
+            self.declare_partials(of=var_names['state_defect'],
+                                  wrt=var_names['f_value'],
+                                  rows=rs_dV, cols=cs_dV, val=d_state_defect_dV.data.ravel())
+
+            self.declare_partials(of=var_names['state_rate_defect'],
+                                  wrt=var_names['f_value'],
+                                  rows=ar1, cols=ar1, val=1.0)
+
+            self.declare_partials(of=var_names['state_rate_defect'],
+                                  wrt=var_names['f_computed'],
+                                  rows=ar1, cols=ar1)
+
+            self.declare_partials(of=var_names['state_rate_defect'],
+                                  wrt='dt_dstau',
+                                  rows=ar1, cols=c1)
+
+            self.declare_partials(of=var_names['initial_state_defect'],
+                                  wrt=var_names['state_initial_value'],
+                                  rows=np.arange(size, dtype=int),
+                                  cols=np.arange(size, dtype=int),
+                                  val=-1.0)
+
+            self.declare_partials(of=var_names['initial_state_defect'],
+                                  wrt=var_names['state_segment_ends'],
+                                  rows=np.arange(size, dtype=int),
+                                  cols=np.arange(size, dtype=int),
+                                  val=1.0)
+
+            self.declare_partials(of=var_names['final_state_defect'],
+                                  wrt=var_names['state_final_value'],
+                                  rows=np.arange(size, dtype=int),
+                                  cols=np.arange(size, dtype=int),
+                                  val=-1.0)
+
+            self.declare_partials(of=var_names['final_state_defect'],
+                                  wrt=var_names['state_segment_ends'],
+                                  rows=np.arange(size, dtype=int),
+                                  cols=size * (2 * num_segs - 1) + np.arange(size, dtype=int),
+                                  val=1.0)
+
+            if num_segs > 1:
+                self.add_output(
+                    name=var_names['state_continuity_defect'],
+                    shape=(num_segs - 1,) + shape,
+                    units=units
+                )
+
+                # Build the derivative matrix for the derivative of state continuity wrt the segment end vals.
+                cnty_deriv_pattern = np.zeros((num_segs - 1,) + (size,) + (num_segs, 2) + (size,))
+                for of_segbound in range(num_segs - 1):
+                    next_seg_idx = of_segbound + 1
+                    prev_seg_idx = of_segbound
+                    for state_idx in range(size):
+                        cnty_deriv_pattern[of_segbound, state_idx, prev_seg_idx, 1, state_idx] = 1.0
+                        cnty_deriv_pattern[of_segbound, state_idx, next_seg_idx, 0, state_idx] = -1.0
+                cnty_deriv = sp.csr_matrix(cnty_deriv_pattern.reshape((num_segs - 1) * size, num_segs * 2 * size))
+                rs, cs = cnty_deriv.nonzero()
+                val = cnty_deriv.data.ravel()
+
+                self.declare_partials(of=var_names['state_continuity_defect'],
+                                      wrt=var_names['state_segment_ends'],
+                                      rows=rs, cols=cs, val=val)
+
+    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
+        """
+        Compute component outputs.
+
+        Parameters
+        ----------
+        inputs : Vector
+            Unscaled, dimensional input variables read via inputs[key].
+        outputs : Vector
+            Unscaled, dimensional output variables read via outputs[key].
+        discrete_inputs : dict or None
+            If not None, dict containing discrete input values.
+        discrete_outputs : dict or None
+            If not None, dict containing discrete output values.
+        """
+        dt_dstau = np.atleast_2d(inputs['dt_dstau']).T
+        num_segs = self.options['grid_data'].num_segments
+
+        for state_name, options in self.options['state_options'].items():
+            var_names = self.var_names[state_name]
+
+            x_a = inputs[var_names['state_initial_value']]
+            x_b = inputs[var_names['state_final_value']]
+            x_ab = inputs[var_names['state_segment_ends']]
+            X = inputs[var_names['state_value']]
+            V = inputs[var_names['f_value']]
+            f = inputs[var_names['f_computed']]
+            shape = x_a.shape[1:]
+            size = x_a.size
+
+            # Get stacked XV, but sorted in segment-by-segment order
+            XV = np.vstack((X, V))[self._xv_idxs]
+
+            state_defect = np.einsum('ij,jk...->ik...', self._A, XV) - \
+                np.einsum('ij, j...->i...', self._C, x_ab.reshape((2 * num_segs,) + shape))
+
+            outputs[var_names['state_defect']] = state_defect
+            outputs[var_names['state_rate_defect']] = (V - np.einsum('i...,i...->i...', f, dt_dstau))
+            outputs[var_names['initial_state_defect']] = x_ab[0, 0, ...] - x_a
+            outputs[var_names['final_state_defect']] = x_ab[-1, 1, ...] - x_b
+
+            if num_segs > 1:
+                (num_segs - 1) * size
+                outputs[var_names['state_continuity_defect']] = x_ab[:-1, 1, ...] - x_ab[1:, 0, ...]
+
+    def compute_partials(self, inputs, partials):
+        """
+        Compute sub-jacobian parts. The model is assumed to be in an unscaled state.
+
+        Parameters
+        ----------
+        inputs : Vector
+            Unscaled, dimensional input variables read via inputs[key].
+        partials : Jacobian
+            Subjac components written to partials[output_name, input_name].
+        """
+        dt_dstau = inputs['dt_dstau']
+
+        for state_name, options in self.options['state_options'].items():
+            var_names = self.var_names[state_name]
+            shape = options['shape']
+            size = np.prod(shape)
+            f = inputs[var_names['f_computed']]
+
+            partials[var_names['state_rate_defect'], var_names['f_computed']] = np.repeat(-dt_dstau, size)
+            partials[var_names['state_rate_defect'], 'dt_dstau'] = -f.ravel()
diff --git a/dymos/transcriptions/pseudospectral/components/birkhoff_iter_group.py b/dymos/transcriptions/pseudospectral/components/birkhoff_iter_group.py
new file mode 100644
index 00000000..4afa3fb8
--- /dev/null
+++ b/dymos/transcriptions/pseudospectral/components/birkhoff_iter_group.py
@@ -0,0 +1,357 @@
+import numpy as np
+import openmdao.api as om
+
+from .birkhoff_collocation_comp import BirkhoffCollocationComp
+from .birkhoff_state_resid_comp import BirkhoffStateResidComp
+
+from ...grid_data import GridData
+from ....phase.options import TimeOptionsDictionary
+
+
+class BirkhoffIterGroup(om.Group):
+    """
+    Class definition for the BirkhoffIterGroup.
+
+    This group allows for iteration of the state variables and initial _or_ final value of the state
+    depending on the direction of the solve.
+
+    Parameters
+    ----------
+    **kwargs : dict
+        Dictionary of optional arguments.
+    """
+    def initialize(self):
+        """
+        Declare group options.
+        """
+        self.options.declare('state_options', types=dict,
+                             desc='Dictionary of options for the states.')
+        self.options.declare('time_options', types=TimeOptionsDictionary,
+                             desc='Options for time in the phase.')
+        self.options.declare('grid_data', types=GridData, desc='Container object for grid info.')
+        self.options.declare('ode_class', default=None,
+                             desc='Callable that instantiates the ODE system.',
+                             recordable=False)
+        self.options.declare('ode_init_kwargs', types=dict, default={},
+                             desc='Keyword arguments provided when initializing the ODE System')
+
+    def setup(self):
+        """
+        Define the structure of the control group.
+        """
+        gd = self.options['grid_data']
+        nn = gd.subset_num_nodes['all']
+        state_options = self.options['state_options']
+        time_options = self.options['time_options']
+        ode_class = self.options['ode_class']
+        ode_init_kwargs = self.options['ode_init_kwargs']
+
+        self.add_subsystem('ode_all', subsys=ode_class(num_nodes=nn, **ode_init_kwargs))
+
+        self.add_subsystem('collocation_comp',
+                           subsys=BirkhoffCollocationComp(grid_data=gd,
+                                                          state_options=state_options,
+                                                          time_units=time_options['units']),
+                           promotes_inputs=['*'], promotes_outputs=['*'])
+
+        if any([opts['solve_segments'] in ('forward', 'backward') for opts in state_options.values()]):
+            self.add_subsystem('states_balance_comp', subsys=BirkhoffStateResidComp(),
+                               promotes_inputs=['*'], promotes_outputs=['*'])
+
+    def _configure_desvars(self, name, options):
+        state_name = f'states:{name}'
+        initial_state_name = f'initial_states:{name}'
+        final_state_name = f'final_states:{name}'
+        state_rate_name = f'state_rates:{name}'
+        state_segment_ends_name = f'state_segment_ends:{name}'
+
+        solve_segs = options['solve_segments']
+        opt = options['opt']
+
+        num_nodes = self.options['grid_data'].subset_num_nodes['col']
+
+        ib = (None, None) if options['initial_bounds'] is None else options['initial_bounds']
+        fb = (None, None) if options['final_bounds'] is None else options['final_bounds']
+        lower = options['lower']
+        upper = options['upper']
+        scaler = options['scaler']
+        adder = options['adder']
+        ref0 = options['ref0']
+        ref = options['ref']
+        fix_initial = options['fix_initial']
+        fix_final = options['fix_final']
+        input_initial = options['input_initial']
+        input_final = options['input_final']
+        shape = options['shape']
+
+        if not np.isscalar(ref0) and ref0 is not None:
+            ref0 = np.asarray(ref0)
+            if ref0.shape == shape:
+                ref0_state = np.tile(ref0.flatten(), num_nodes)
+                ref0_seg_ends = np.tile(ref0.flatten(), 2)
+            else:
+                raise ValueError('array-valued scaler/ref must length equal to state-size')
+        else:
+            ref0_state = ref0
+            ref0_seg_ends = ref0
+        if not np.isscalar(ref) and ref is not None:
+            ref = np.asarray(ref)
+            if ref.shape == shape:
+                ref_state = np.tile(ref.flatten(), num_nodes)
+                ref_seg_ends = np.tile(ref.flatten(), 2)
+            else:
+                raise ValueError('array-valued scaler/ref must length equal to state-size')
+        else:
+            ref_state = ref
+            ref_seg_ends = ref
+
+        free_vars = {state_name, state_rate_name, state_segment_ends_name, initial_state_name, final_state_name}
+
+        if solve_segs == 'forward':
+            implicit_outputs = {state_name, state_rate_name, final_state_name, state_segment_ends_name}
+        elif solve_segs == 'backward':
+            implicit_outputs = {state_name, state_rate_name, initial_state_name, state_segment_ends_name}
+        else:
+            implicit_outputs = set()
+
+        free_vars = free_vars - implicit_outputs
+
+        if fix_initial or input_initial:
+            free_vars = free_vars - {initial_state_name}
+        if fix_final or input_final:
+            free_vars = free_vars - {final_state_name}
+
+        if opt:
+            # Add design variables for the remaining free variables
+            if state_name in free_vars:
+                self.add_design_var(name=state_name,
+                                    lower=lower,
+                                    upper=upper,
+                                    scaler=scaler,
+                                    adder=adder,
+                                    ref0=ref0_state,
+                                    ref=ref_state)
+
+            if state_segment_ends_name in free_vars:
+                self.add_design_var(name=state_segment_ends_name,
+                                    lower=lower,
+                                    upper=upper,
+                                    scaler=scaler,
+                                    adder=adder,
+                                    ref0=ref0_seg_ends,
+                                    ref=ref_seg_ends)
+
+            if state_rate_name in free_vars:
+                self.add_design_var(name=state_rate_name,
+                                    scaler=scaler,
+                                    adder=adder,
+                                    ref0=ref0_state,
+                                    ref=ref_state)
+
+            if initial_state_name in free_vars:
+                self.add_design_var(name=initial_state_name,
+                                    lower=ib[0],
+                                    upper=ib[1],
+                                    scaler=scaler,
+                                    adder=adder,
+                                    ref0=ref0,
+                                    ref=ref)
+
+            if final_state_name in free_vars:
+                self.add_design_var(name=final_state_name,
+                                    lower=fb[0],
+                                    upper=fb[1],
+                                    scaler=scaler,
+                                    adder=adder,
+                                    ref0=ref0,
+                                    ref=ref)
+
+        return implicit_outputs
+
+    def configure_io(self, phase):
+        """
+        I/O creation is delayed until configure so that we can determine shape and units for the states.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        collocation_comp = self._get_subsystem('collocation_comp')
+        collocation_comp.configure_io(phase)
+
+        gd = self.options['grid_data']
+        nn = gd.subset_num_nodes['all']
+        ns = gd.num_segments
+
+        state_options = self.options['state_options']
+        time_options = self.options['time_options']
+        time_units = time_options['units']
+        states_balance_comp = self._get_subsystem('states_balance_comp')
+
+        for name, options in state_options.items():
+            units = options['units']
+            rate_source = options['rate_source']
+            shape = options['shape']
+
+            for tgt in options['targets']:
+                self.promotes('ode_all', [(tgt, f'states:{name}')])
+                self.set_input_defaults(f'states:{name}', val=1.0, units=units, src_shape=(nn,) + shape)
+
+            implicit_outputs = self._configure_desvars(name, options)
+
+            if f'states:{name}' in implicit_outputs:
+                states_balance_comp.add_output(f'states:{name}',
+                                               shape=(nn,) + shape,
+                                               units=units)
+
+                states_balance_comp.add_output(f'state_segment_ends:{name}',
+                                               shape=(ns, 2) + shape,
+                                               units=units)
+
+                states_balance_comp.add_residual_from_input(f'state_defects:{name}',
+                                                            shape=(nn+ns,) + shape,
+                                                            units=units)
+
+                states_balance_comp.add_residual_from_input(f'initial_state_defects:{name}',
+                                                            shape=(1,) + shape,
+                                                            units=units)
+
+                states_balance_comp.add_residual_from_input(f'final_state_defects:{name}',
+                                                            shape=(1,) + shape,
+                                                            units=units)
+
+            if f'state_rates:{name}' in implicit_outputs:
+                states_balance_comp.add_output(f'state_rates:{name}', shape=(nn,) + shape, units=units)
+                states_balance_comp.add_residual_from_input(f'state_rate_defects:{name}',
+                                                            shape=(nn,) + shape,
+                                                            units=units)
+
+            if f'initial_states:{name}' in implicit_outputs:
+                states_balance_comp.add_output(f'initial_states:{name}', shape=(1,) + shape, units=units)
+
+            if f'final_states:{name}' in implicit_outputs:
+                states_balance_comp.add_output(f'final_states:{name}', shape=(1,) + shape, units=units)
+
+            try:
+                rate_source_var = options['rate_source']
+            except RuntimeError:
+                raise ValueError(f"state '{name}' in phase '{phase.name}' was not given a "
+                                 "rate_source")
+
+            # Note the rate source must be shape-compatible with the state
+            var_type = phase.classify_var(rate_source_var)
+
+            if var_type == 'ode':
+                self.connect(f'ode_all.{rate_source}', f'f_computed:{name}')
+
+            if implicit_outputs:
+                if isinstance(phase.nonlinear_solver, om.NonlinearRunOnce):
+                    newton = phase.nonlinear_solver = om.NewtonSolver()
+                    newton.options['solve_subsystems'] = True
+                    newton.options['maxiter'] = 100
+                    newton.options['iprint'] = 2
+                    newton.options['stall_limit'] = 3
+                    newton.linesearch = om.BoundsEnforceLS()
+
+                if isinstance(phase.linear_solver, om.LinearRunOnce):
+                    phase.linear_solver = om.DirectSolver()
+
+    def _get_rate_source_path(self, state_name, nodes, phase):
+        """
+        Return the rate source location and indices for a given state name.
+
+        Parameters
+        ----------
+        state_name : str
+            Name of the state.
+        nodes : str
+            One of ['col', 'all'].
+        phase : dymos.Phase
+            Phase object containing the rate source.
+
+        Returns
+        -------
+        str
+            Path to the rate source.
+        ndarray
+            Array of source indices.
+        """
+        gd = self.grid_data
+        try:
+            var = phase.state_options[state_name]['rate_source']
+        except RuntimeError:
+            raise ValueError(f"state '{state_name}' in phase '{phase.name}' was not given a "
+                             "rate_source")
+
+        # Note the rate source must be shape-compatible with the state
+        var_type = phase.classify_var(var)
+
+        # Determine the path to the variable
+        if var_type == 't':
+            rate_path = 't'
+            node_idxs = gd.subset_node_indices[nodes]
+        elif var_type == 't_phase':
+            rate_path = 't_phase'
+            node_idxs = gd.subset_node_indices[nodes]
+        elif var_type == 'state':
+            rate_path = f'states:{var}'
+            # Find the state_input indices which occur at segment endpoints, and repeat them twice
+            state_input_idxs = gd.subset_node_indices['state_input']
+            repeat_idxs = np.ones_like(state_input_idxs)
+            if self.options['compressed']:
+                segment_end_idxs = gd.subset_node_indices['segment_ends'][1:-1]
+                # Repeat nodes that are on segment bounds (but not the first or last nodes in the phase)
+                nodes_to_repeat = list(set(state_input_idxs).intersection(segment_end_idxs))
+                # Now find these nodes in the state input indices
+                idxs_of_ntr_in_state_inputs = np.where(np.in1d(state_input_idxs, nodes_to_repeat))[0]
+                # All state input nodes are used once, but nodes_to_repeat are used twice
+                repeat_idxs[idxs_of_ntr_in_state_inputs] = 2
+            # Now we have a way of mapping the state input indices to all nodes
+            map_input_node_idxs_to_all = np.repeat(np.arange(gd.subset_num_nodes['state_input'],
+                                                             dtype=int), repeats=repeat_idxs)
+            # Now select the subset of nodes we want to use.
+            node_idxs = map_input_node_idxs_to_all[gd.subset_node_indices[nodes]]
+        elif var_type == 'indep_control':
+            rate_path = f'control_values:{var}'
+            node_idxs = gd.subset_node_indices[nodes]
+        elif var_type == 'input_control':
+            rate_path = f'control_values:{var}'
+            node_idxs = gd.subset_node_indices[nodes]
+        elif var_type == 'control_rate':
+            control_name = var[:-5]
+            rate_path = f'control_rates:{control_name}_rate'
+            node_idxs = gd.subset_node_indices[nodes]
+        elif var_type == 'control_rate2':
+            control_name = var[:-6]
+            rate_path = f'control_rates:{control_name}_rate2'
+            node_idxs = gd.subset_node_indices[nodes]
+        elif var_type == 'indep_polynomial_control':
+            rate_path = f'polynomial_control_values:{var}'
+            node_idxs = gd.subset_node_indices[nodes]
+        elif var_type == 'input_polynomial_control':
+            rate_path = f'polynomial_control_values:{var}'
+            node_idxs = gd.subset_node_indices[nodes]
+        elif var_type == 'polynomial_control_rate':
+            control_name = var[:-5]
+            rate_path = f'polynomial_control_rates:{control_name}_rate'
+            node_idxs = gd.subset_node_indices[nodes]
+        elif var_type == 'polynomial_control_rate2':
+            control_name = var[:-6]
+            rate_path = f'polynomial_control_rates:{control_name}_rate2'
+            node_idxs = gd.subset_node_indices[nodes]
+        elif var_type == 'parameter':
+            rate_path = f'parameter_vals:{var}'
+            dynamic = not phase.parameter_options[var]['static_target']
+            if dynamic:
+                node_idxs = np.zeros(gd.subset_num_nodes[nodes], dtype=int)
+            else:
+                node_idxs = np.zeros(1, dtype=int)
+        else:
+            # Failed to find variable, assume it is in the ODE
+            rate_path = f'ode_all.{var}'
+            node_idxs = gd.subset_node_indices[nodes]
+
+        src_idxs = om.slicer[node_idxs, ...]
+
+        return rate_path, src_idxs
diff --git a/dymos/transcriptions/pseudospectral/components/birkhoff_state_resid_comp.py b/dymos/transcriptions/pseudospectral/components/birkhoff_state_resid_comp.py
new file mode 100644
index 00000000..fdb3286c
--- /dev/null
+++ b/dymos/transcriptions/pseudospectral/components/birkhoff_state_resid_comp.py
@@ -0,0 +1,67 @@
+import numpy as np
+import openmdao.api as om
+
+from openmdao.utils.general_utils import ensure_compatible
+
+from ...._options import options as dymos_options
+
+
+class BirkhoffStateResidComp(om.ImplicitComponent):
+    """
+    Class definition for the BirkhoffStateResidComp.
+
+    Generates the residuals for any states that are solved for implicitly.
+
+    Parameters
+    ----------
+    **kwargs : dict
+        Dictionary of optional arguments.
+    """
+
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        self._no_check_partials = not dymos_options['include_check_partials']
+        self._io_pairs = []
+
+    def add_residual_from_input(self, name, **kwargs):
+        """
+        Adds a residual whose value is given by resid_input.
+
+        Parameters
+        ----------
+        name : str
+            The name of the input providing the residuals for the given output.
+        **kwargs : dict
+            Additional keyword arguments for add_input and add_residual.
+        """
+        val = kwargs['val'] if 'val' in kwargs else 1.0
+        shape = kwargs['shape'] if 'shape' in kwargs else None
+        val, shape = ensure_compatible(name, value=val, shape=shape)
+        resid_name = 'resid_' + name
+
+        self._io_pairs.append((resid_name, name))
+
+        size = np.prod(shape, dtype=int)
+        ar = np.arange(size, dtype=int)
+
+        self.add_input(name, **kwargs)
+        self.add_residual(resid_name, **kwargs)
+        # self.declare_partials(of=resid_name, wrt=name, rows=ar, cols=ar, val=1.0)
+        self.declare_partials(of='*', wrt='*', method='fd')
+
+    def apply_nonlinear(self, inputs, outputs, residuals):
+        """
+        Compute residuals given inputs and outputs.
+
+        The model is assumed to be in an unscaled state.
+
+        Parameters
+        ----------
+        inputs : Vector
+            Unscaled, dimensional input variables read via inputs[key].
+        outputs : Vector
+            Unscaled, dimensional output variables read via outputs[key].
+        residuals : Vector
+            Unscaled, dimensional residuals written to via residuals[key].
+        """
+        residuals.set_val(inputs.asarray())
diff --git a/dymos/transcriptions/transcription_base.py b/dymos/transcriptions/transcription_base.py
index 3e182d58..db64a0a0 100644
--- a/dymos/transcriptions/transcription_base.py
+++ b/dymos/transcriptions/transcription_base.py
@@ -344,13 +344,13 @@ class TranscriptionBase(object):
             for ts_name, ts_options in phase._timeseries.items():
                 if f'{state_prefix}{name}' not in ts_options['outputs']:
                     phase.add_timeseries_output(name, output_name=f'{state_prefix}{name}',
-                                                timeseries=ts_name)
+                                                timeseries=ts_name, tags='state')
                 if options['rate_source'] and phase.timeseries_options['include_state_rates']:
                     output_name = f'{state_rate_prefix}{name}' if state_rate_prefix else options['rate_source']
                     if output_name not in ts_options['outputs']:
                         phase.add_timeseries_output(name=options['rate_source'],
                                                     output_name=output_name,
-                                                    timeseries=ts_name)
+                                                    timeseries=ts_name, tags='state_rate')
 
     def setup_ode(self, phase):
         """
diff --git a/dymos/utils/birkhoff.py b/dymos/utils/birkhoff.py
new file mode 100644
index 00000000..dbec5bf1
--- /dev/null
+++ b/dymos/utils/birkhoff.py
@@ -0,0 +1,69 @@
+import numpy as np
+import scipy.special as sp
+
+
+def birkhoff_matrix(tau, w, grid_type):
+    """
+    Returns the pseudospectral integration matrix for a Birkhoff polynomial at the given nodes.
+
+    Parameters
+    ----------
+    tau : ndarray[:]
+        Vector of given nodes in the polynomial.
+    w : ndarray[:]
+        Vector of nodes at which the polynomial is evaluated.
+    grid_type : str
+        The type of Gaussian grid used in the transcription.
+
+    Returns
+    -------
+    np.array
+        An num_disc_nodes-1 x num_disc_nodes matrix used for the interpolation of state values
+        at the interior LGL nodes.
+    """
+    N = tau.size - 1
+    end_node = N if tau[-1] == 1 else N + 1
+
+    alpha = np.zeros((N + 1, N + 1))
+    S = np.zeros((N + 1, N + 1))
+
+    if grid_type[0] == 'l':
+        pol = sp.eval_legendre
+    elif grid_type[0] == 'c':
+        pol = sp.eval_chebyt
+    else:
+        raise ValueError('invalid grid type')
+
+    for j in range(0, N + 1):
+        for n in range(0, end_node):
+            alpha[n, j] = w[j] * pol(n, tau[j])
+        if grid_type == 'lgl':
+            alpha[N, j] = N * pol(N, tau[j]) * w[j] / (2 * N + 1)
+        elif grid_type == 'cgl':
+            alpha[N, j] = pol(N, tau[j]) * w[j] / 2
+
+    if grid_type[0] == 'l':
+        for i in range(0, N + 1):
+            S[i, 0] = (tau[i] - tau[0]) / 2
+
+            for n in range(1, N+1):
+                gamma = 2 / (2 * n + 1)
+                int_p = (pol(n+1, tau[i]) - pol(n-1, tau[i])) / (2*n+1)
+                S[i, n] = int_p / gamma
+
+    elif grid_type[0] == 'c':
+        gamma = np.pi / 2
+        for i in range(0, N+1):
+            # chebyshev polynomial of order 0: 1
+            S[i, 0] = (tau[i] - tau[0]) / np.pi
+
+            # chebyshev polynomial of order 1: x
+            S[i, 1] = (tau[i]**2 - tau[0]**2) / np.pi
+
+            for n in range(2, N+1):
+                int_p = pol(n+1, tau[i]) / (2*n+2) - pol(n-1, tau[i]) / (2*n-2) - (-1)**n / (n**2 - 1)
+                S[i, n] = int_p / gamma
+
+    B = S @ alpha
+
+    return B
diff --git a/dymos/utils/cgl.py b/dymos/utils/cgl.py
new file mode 100644
index 00000000..db992241
--- /dev/null
+++ b/dymos/utils/cgl.py
@@ -0,0 +1,95 @@
+import numpy as np
+import scipy.special as sp
+
+_cgl_cache = {}
+""" Cache for the LGL nodes and weights, keyed by n. """
+
+
+def clenshaw_curtis(n):
+    """
+    Returns the Chebyshev-Gauss-Lobatto nodes and weights for a Jacobi Polynomial with n abscissae.
+
+    The nodes are on the range [-1, 1].
+
+    Parameters
+    ----------
+    n : int
+        The number of LGL nodes requested.  The order of the polynomial is n-1.
+
+    Returns
+    -------
+    x : numpy.array
+        An array of the LGL nodes for a polynomial of the given order.
+
+    w : numpy.array
+        An array of the corresponding CGL weights at the nodes in x.
+    """
+    x = -np.cos(np.pi * np.arange(n) / (n-1))
+    w = np.zeros(n)
+    even = False if n % 2 else True
+
+    N_2 = int(n/2) if even else int((n+1)/2)
+    j_lim = N_2 if even else int((n-1)/2)
+
+    for k in range(n):
+        c = 1 if k == 0 else 2
+        s = 0
+        for j in range(1, N_2):
+            b = 1 if j == j_lim else 2
+            s += b * np.cos(2*j*k*np.pi/(n-1)) / (4*j**2 - 1)
+        w[k] = c * (1 - s) / (n-1)
+
+    if even:
+        w[N_2:] = np.flip(w[:N_2])
+    else:
+        w[N_2:] = np.flip(w[:(N_2-1)])
+
+    return x, w
+
+
+def _cgl(n):
+    """
+    Returns the Chebyshev-Gauss-Lobatto nodes and weights for a Jacobi Polynomial with n abscissae.
+
+    The nodes are on the range [-1, 1].
+
+    Parameters
+    ----------
+    n : int
+        The number of CGL nodes requested.  The order of the polynomial is n-1.
+
+    Returns
+    -------
+    x : numpy.array
+        An array of the CGL nodes for a polynomial of the given order.
+
+    w : numpy.array
+        An array of the corresponding CGL weights at the nodes in x.
+    """
+    x = -np.cos(np.pi * np.arange(n) / (n-1))
+    w = np.pi / (n-1) * np.ones(n)
+    w[0] *= 0.5
+    w[-1] *= 0.5
+
+    return x, w
+
+
+def cgl(n):
+    """
+    Retrieve the cgl nodes and weights for n nodes.
+
+    Results are cached to avoid repeated calculation of nodes and weights for a given n.
+
+    Parameters
+    ----------
+    n : int
+        Node number.
+
+    Returns
+    -------
+    float
+        Tuple with cgl nodes and weights.
+    """
+    if n not in _cgl_cache:
+        _cgl_cache[n] = _cgl(n)
+    return _cgl_cache[n]
diff --git a/dymos/utils/introspection.py b/dymos/utils/introspection.py
index e76087be..b24cc1d5 100644
--- a/dymos/utils/introspection.py
+++ b/dymos/utils/introspection.py
@@ -173,6 +173,9 @@ def _configure_constraint_introspection(phase):
     phase : Phase
         The phase object whose boundary and path constraints are to be introspected.
     """
+    from ..transcriptions import Birkhoff
+    birkhoff = isinstance(phase.options['transcription'], Birkhoff)
+
     for constraint_type, constraints in [('initial', phase._initial_boundary_constraints),
                                          ('final', phase._final_boundary_constraints),
                                          ('path', phase._path_constraints)]:
@@ -205,7 +208,10 @@ def _configure_constraint_introspection(phase):
                 state_units = phase.state_options[var]['units']
                 con['shape'] = state_shape
                 con['units'] = state_units if con['units'] is None else con['units']
-                con['constraint_path'] = f'timeseries.{prefix}{var}'
+                if birkhoff and constraint_type in ('initial', 'final'):
+                    con['constraint_path'] = f'{constraint_type}_states:{var}'
+                else:
+                    con['constraint_path'] = f'timeseries.{prefix}{var}'
 
             elif var_type == 'parameter':
                 param_shape = phase.parameter_options[var]['shape']
diff --git a/dymos/utils/misc.py b/dymos/utils/misc.py
index 1576205c..fd55f736 100644
--- a/dymos/utils/misc.py
+++ b/dymos/utils/misc.py
@@ -144,7 +144,7 @@ class CoerceDesvar(object):
                              'as states ({1})'.format(option, self.options['shape']))
 
 
-def CompWrapperConfig(comp_class):
+def CompWrapperConfig(comp_class, config_io_args=None):
     """
     Returns a wrapped comp_class that calls its configure_io method at the end of setup.
 
@@ -154,7 +154,9 @@ def CompWrapperConfig(comp_class):
     Parameters
     ----------
     comp_class : Component class
-       Class that we would like to wrap.
+        Class that we would like to wrap.
+    config_io_args : list
+        Arguments to be passed to config_io.
 
     Returns
     -------
@@ -168,13 +170,14 @@ def CompWrapperConfig(comp_class):
             Appends a call to configure_io after setup.
             """
             super(WrappedClass, self).setup()
-            self.configure_io()
+            args = [] if config_io_args is None else config_io_args
+            self.configure_io(*args)
 
     return WrappedClass
 
 
 # Modify class so we can run it standalone.
-def GroupWrapperConfig(comp_class):
+def GroupWrapperConfig(comp_class, config_io_args=None):
     """
     Returns a wrapped group_class that calls its configure_io method at the end of setup.
 
@@ -185,6 +188,8 @@ def GroupWrapperConfig(comp_class):
     ----------
     comp_class : Group class
        Class that we would like to wrap.
+    config_io_args : list
+        Arguments to be passed to config_io.
 
     Returns
     -------
@@ -203,6 +208,7 @@ def GroupWrapperConfig(comp_class):
             """
             Call configure_io during configure.
             """
-            self.configure_io()
+            args = [] if config_io_args is None else config_io_args
+            self.configure_io(*args)
 
     return WrappedClass
diff --git a/setup.py b/setup.py
index 47f6f6f5..07cfbc0e 100644
--- a/setup.py
+++ b/setup.py
@@ -73,7 +73,7 @@ The software has two primary objectives:
     packages=find_packages(),
     python_requires=">=3.8",
     install_requires=[
-        'openmdao>=3.26.0',
+        'openmdao>=3.27.0',
         'numpy',
         'scipy'
     ],


Test patch:
diff --git a/.github/workflows/dymos_tests_workflow.yml b/.github/workflows/dymos_tests_workflow.yml
index a75df368..3dd1c75a 100644
--- a/.github/workflows/dymos_tests_workflow.yml
+++ b/.github/workflows/dymos_tests_workflow.yml
@@ -88,7 +88,7 @@ jobs:
             PETSc: 3.13
             PYOPTSPARSE: 'v2.6.1'
             SNOPT: 7.2
-            OPENMDAO: 3.26.0
+            OPENMDAO: 3.27.0
             OPTIONAL: '[test]'
 
     steps:
@@ -312,10 +312,14 @@ jobs:
           testflo -n 1 joss/test --pre_announce
           testflo -b benchmark --pre_announce
           cd $HOME
-          testflo dymos -n 2 --show_skipped --coverage --durations 20 --coverpkg dymos
+          if [[ "${{ matrix.NAME }}" != "latest" ]]; then
+            testflo dymos -n 2 --show_skipped --durations 20 --coverage --coverpkg dymos
+          else
+            testflo dymos -n 2 --show_skipped --durations 20
+          fi
 
       - name: Submit coverage
-        if: (github.event_name != 'workflow_dispatch')
+        if: ((github.event_name != 'workflow_dispatch') && (matrix.NAME != 'latest'))
         shell: bash -l {0}
         env:
           COVERALLS_REPO_TOKEN: ${{ secrets.GITHUB_TOKEN }}
diff --git a/dymos/examples/aircraft_steady_flight/test/test_aircraft_cruise.py b/dymos/examples/aircraft_steady_flight/test/test_aircraft_cruise.py
index ca08722c..fc1e9ca6 100644
--- a/dymos/examples/aircraft_steady_flight/test/test_aircraft_cruise.py
+++ b/dymos/examples/aircraft_steady_flight/test/test_aircraft_cruise.py
@@ -208,6 +208,106 @@ class TestAircraftCruise(unittest.TestCase):
 
         assert_near_equal(range, tas*time, tolerance=1.0E-4)
 
+    def test_cruise_results_birkhoff(self):
+        p = om.Problem(model=om.Group())
+        if optimizer == 'SNOPT':
+            p.driver = om.pyOptSparseDriver()
+            p.driver.options['optimizer'] = optimizer
+            p.driver.declare_coloring()
+            p.driver.opt_settings['Major iterations limit'] = 100
+            p.driver.opt_settings['Major step limit'] = 0.05
+            p.driver.opt_settings['Major feasibility tolerance'] = 1.0E-6
+            p.driver.opt_settings['Major optimality tolerance'] = 1.0E-6
+            p.driver.opt_settings["Linesearch tolerance"] = 0.10
+            p.driver.opt_settings['iSumm'] = 6
+            p.driver.opt_settings['Verify level'] = 3
+        else:
+            p.driver = om.ScipyOptimizeDriver()
+            p.driver.declare_coloring()
+
+        # Pass Reference Area from an external source
+        assumptions = p.model.add_subsystem('assumptions', om.IndepVarComp())
+        assumptions.add_output('S', val=427.8, units='m**2')
+        assumptions.add_output('mass_empty', val=1.0, units='kg')
+        assumptions.add_output('mass_payload', val=1.0, units='kg')
+
+        grid = dm.BirkhoffGrid(num_segments=1, nodes_per_seg=10, grid_type='lgl')
+
+        transcription = dm.Birkhoff(grid=grid)
+        phase = dm.Phase(ode_class=AircraftODE, transcription=transcription)
+        p.model.add_subsystem('phase0', phase)
+
+        phase.set_time_options(fix_initial=True,
+                               duration_bounds=(3600, 3600),
+                               duration_ref=3600)
+
+        phase.set_time_options(fix_initial=True,
+                               duration_bounds=(3600, 3600),
+                               duration_ref=3600)
+
+        phase.add_state('range', units='km', fix_initial=True, fix_final=False, scaler=0.01,
+                        rate_source='range_rate_comp.dXdt:range',
+                        defect_scaler=0.01)
+        phase.add_state('mass_fuel', units='kg', fix_final=True, upper=20000.0, lower=0.0,
+                        rate_source='propulsion.dXdt:mass_fuel',
+                        scaler=1.0E-4, defect_scaler=1.0E-2)
+        phase.add_state('alt',
+                        rate_source='climb_rate',
+                        units='km', fix_initial=True)
+
+        phase.add_control('mach', targets=['tas_comp.mach', 'aero.mach'], opt=False)
+
+        phase.add_control('climb_rate', targets=['gam_comp.climb_rate'], units='m/s', opt=False)
+
+        phase.add_parameter('S',
+                            targets=['aero.S', 'flight_equilibrium.S', 'propulsion.S'],
+                            units='m**2')
+
+        phase.add_parameter('mass_empty', targets=['mass_comp.mass_empty'], units='kg')
+        phase.add_parameter('mass_payload', targets=['mass_comp.mass_payload'], units='kg')
+
+        phase.add_path_constraint('propulsion.tau', lower=0.01, upper=1.0)
+
+        phase.add_timeseries_output('tas_comp.TAS')
+
+        p.model.connect('assumptions.S', 'phase0.parameters:S')
+        p.model.connect('assumptions.mass_empty', 'phase0.parameters:mass_empty')
+        p.model.connect('assumptions.mass_payload', 'phase0.parameters:mass_payload')
+
+        phase.add_objective('time', loc='final', ref=3600)
+
+        p.model.linear_solver = om.DirectSolver()
+
+        p.setup()
+
+        p['phase0.t_initial'] = 0.0
+        p['phase0.t_duration'] = 1.0 * 3600.0
+        p['phase0.initial_states:range'] = 0.0
+        p['phase0.final_states:mass_fuel'] = 0.0
+        p['phase0.initial_states:alt'] = 5.0
+        p['phase0.controls:mach'] = 0.8
+        p['phase0.controls:climb_rate'] = 0.0
+
+        p['assumptions.S'] = 427.8
+        p['assumptions.mass_empty'] = 0.15E6
+        p['assumptions.mass_payload'] = 84.02869 * 400
+
+        dm.run_problem(p)
+
+        time = p.get_val('phase0.timeseries.time')
+        tas = p.get_val('phase0.timeseries.TAS', units='km/s')
+        range = p.get_val('phase0.timeseries.range')
+
+        assert_near_equal(range, tas*time, tolerance=1.0E-4)
+
+        exp_out = phase.simulate()
+
+        time = exp_out.get_val('phase0.timeseries.time')
+        tas = exp_out.get_val('phase0.timeseries.TAS', units='km/s')
+        range = exp_out.get_val('phase0.timeseries.range')
+
+        assert_near_equal(range, tas*time, tolerance=1.0E-4)
+
 
 if __name__ == '__main__':  # pragma: no cover
     unittest.main()
diff --git a/dymos/examples/brachistochrone/test/ex_brachistochrone.py b/dymos/examples/brachistochrone/test/ex_brachistochrone.py
index 3fcc5d4f..6b461377 100644
--- a/dymos/examples/brachistochrone/test/ex_brachistochrone.py
+++ b/dymos/examples/brachistochrone/test/ex_brachistochrone.py
@@ -1,4 +1,5 @@
 import matplotlib
+import numpy as np
 
 import openmdao.api as om
 from openmdao.utils.testing_utils import require_pyoptsparse
@@ -39,13 +40,19 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran
                             nodes_per_seg=transcription_order + 1,
                             compressed=compressed)
         t = dm.ExplicitShooting(grid=grid)
+    elif transcription == 'birkhoff':
+        grid = dm.BirkhoffGrid(num_segments=num_segments,
+                               nodes_per_seg=transcription_order + 1,
+                               compressed=compressed, grid_type='cgl')
+        t = dm.Birkhoff(grid=grid)
+        # phase = dm.ImplicitPhase(ode_class=BrachistochroneODE, num_nodes=11)
 
     traj = dm.Trajectory()
     phase = dm.Phase(ode_class=BrachistochroneODE, transcription=t)
     p.model.add_subsystem('traj0', traj)
     traj.add_phase('phase0', phase)
 
-    phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+    phase.set_time_options(fix_initial=True, duration_bounds=(0.5, 10))
 
     phase.add_state('x', fix_initial=True, fix_final=False, solve_segments=solve_segments)
     phase.add_state('y', fix_initial=True, fix_final=False, solve_segments=solve_segments)
@@ -68,9 +75,12 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran
 
     phase.add_boundary_constraint('x', loc='final', equals=10)
     phase.add_boundary_constraint('y', loc='final', equals=5)
+    phase.add_path_constraint('y', lower=0, upper=20)
     # Minimize time at the end of the phase
     phase.add_objective('time_phase', loc='final', scaler=10)
 
+    phase.add_timeseries_output('check')
+
     p.setup(check=['unconnected_inputs'], force_alloc_complex=force_alloc_complex)
 
     phase.set_simulate_options(method='RK23')
@@ -78,7 +88,7 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran
     p['traj0.phase0.t_initial'] = 0.0
     p['traj0.phase0.t_duration'] = 2.0
 
-    if transcription.startswith('shooting'):
+    if transcription.startswith('shooting') or transcription == 'birkhoff':
         p['traj0.phase0.initial_states:x'] = 0
         p['traj0.phase0.initial_states:y'] = 10
         p['traj0.phase0.initial_states:v'] = 0
@@ -90,12 +100,15 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran
     p['traj0.phase0.controls:theta'] = phase.interp('theta', [5, 100])
     p['traj0.phase0.parameters:g'] = 9.80665
 
-    dm.run_problem(p, run_driver=run_driver, simulate=True, make_plots=True)
+    dm.run_problem(p, run_driver=run_driver, simulate=True, make_plots=True,
+                   simulate_kwargs={'times_per_seg': 100})
 
     return p
 
 
 if __name__ == '__main__':
-    p = brachistochrone_min_time(transcription='gauss-lobatto', num_segments=5, run_driver=True,
-                                 transcription_order=5, compressed=False, optimizer='IPOPT',
-                                 solve_segments=False, force_alloc_complex=True)
+
+    with dm.options.temporary(include_check_partials=True):
+        p = brachistochrone_min_time(transcription='birkhoff', num_segments=1, run_driver=True,
+                                     transcription_order=19, compressed=False, optimizer='SLSQP',
+                                     solve_segments=False, force_alloc_complex=True)
diff --git a/dymos/examples/brachistochrone/test/ex_brachistochrone_vector_states.py b/dymos/examples/brachistochrone/test/ex_brachistochrone_vector_states.py
index 714fcf3a..d3fa8cff 100644
--- a/dymos/examples/brachistochrone/test/ex_brachistochrone_vector_states.py
+++ b/dymos/examples/brachistochrone/test/ex_brachistochrone_vector_states.py
@@ -5,7 +5,7 @@ from dymos.examples.brachistochrone.brachistochrone_vector_states_ode \
 
 
 def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, transcription_order=3,
-                             compressed=True, optimizer='SLSQP',
+                             grid_type='lgl', compressed=True, optimizer='SLSQP',
                              dynamic_simul_derivs=True, force_alloc_complex=False,
                              solve_segments=False, run_driver=True):
     p = om.Problem(model=om.Group())
@@ -33,6 +33,11 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran
                                  order=transcription_order,
                                  compressed=compressed)
         fix_final = not solve_segments
+    elif transcription == 'birkhoff':
+        gd = dm.BirkhoffGrid(num_segments=num_segments, nodes_per_seg=transcription_order,
+                             grid_type=grid_type)
+        transcription = dm.Birkhoff(grid=gd)
+        fix_final = not solve_segments
 
     traj = dm.Trajectory()
     phase = dm.Phase(ode_class=BrachistochroneVectorStatesODE,
@@ -45,7 +50,8 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran
 
     # can't fix final position if you're solving the segments
 
-    phase.add_state('pos', fix_initial=True, fix_final=fix_final, solve_segments=solve_segments, ref=[1, 1])
+    phase.add_state('pos', fix_initial=True, fix_final=fix_final,
+                    solve_segments=solve_segments, ref=[1, 1])
     #
     phase.add_state('v', fix_initial=True, fix_final=False, solve_segments=solve_segments)
     #
@@ -74,6 +80,11 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran
     p['traj0.phase0.states:v'] = phase.interp('v', [0, 9.9])
     p['traj0.phase0.controls:theta'] = phase.interp('theta', [5, 100])
     p['traj0.phase0.parameters:g'] = 9.80665
+    if isinstance(transcription, dm.Birkhoff):
+        p['traj0.phase0.initial_states:pos'] = pos0
+        p['traj0.phase0.initial_states:v'] = 0.0
+        p['traj0.phase0.final_states:pos'] = posf
+        p['traj0.phase0.final_states:v'] = 9.9
 
     p.run_model()
     if run_driver:
@@ -85,4 +96,5 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran
 if __name__ == '__main__':
     p = brachistochrone_min_time(transcription='radau-ps', num_segments=5, run_driver=True,
                                  transcription_order=5, compressed=False, optimizer='SLSQP',
-                                 solve_segments=False, force_alloc_complex=True, dynamic_simul_derivs=True)
+                                 solve_segments=False, force_alloc_complex=True,
+                                 dynamic_simul_derivs=True)
diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_solve_segments.py b/dymos/examples/brachistochrone/test/test_brachistochrone_solve_segments.py
index 065d16ea..f0909489 100644
--- a/dymos/examples/brachistochrone/test/test_brachistochrone_solve_segments.py
+++ b/dymos/examples/brachistochrone/test/test_brachistochrone_solve_segments.py
@@ -204,6 +204,33 @@ class TestBrachistochroneVectorStatesExampleSolveSegments(unittest.TestCase):
                                                            transcription_order=11)
         self.assert_results(p)
 
+    def test_ex_brachistochrone_vs_birkhoff_lgl(self):
+        ex_brachistochrone_vs.SHOW_PLOTS = False
+        p = ex_brachistochrone_vs.brachistochrone_min_time(transcription='birkhoff',
+                                                           force_alloc_complex=True,
+                                                           solve_segments=False,
+                                                           grid_type='lgl',
+                                                           num_segments=1,
+                                                           transcription_order=15,
+                                                           optimizer='SLSQP')
+
+        cpd = p.check_partials(compact_print=True, method='cs')
+        from dymos.utils.testing_utils import assert_check_partials
+        assert_check_partials(cpd)
+
+        self.assert_results(p)
+
+    def test_ex_brachistochrone_vs_birkhoff_cgl(self):
+        ex_brachistochrone_vs.SHOW_PLOTS = False
+        p = ex_brachistochrone_vs.brachistochrone_min_time(transcription='birkhoff',
+                                                           compressed=False,
+                                                           force_alloc_complex=True,
+                                                           solve_segments='forward',
+                                                           grid_type='cgl',
+                                                           num_segments=1,
+                                                           transcription_order=11)
+        self.assert_results(p)
+
 
 @require_pyoptsparse(optimizer='IPOPT')
 @use_tempdirs
diff --git a/dymos/examples/cannonball/test/test_cannonball_matrix_state.py b/dymos/examples/cannonball/test/test_cannonball_matrix_state.py
index cf1cfd6b..b7aeec00 100644
--- a/dymos/examples/cannonball/test/test_cannonball_matrix_state.py
+++ b/dymos/examples/cannonball/test/test_cannonball_matrix_state.py
@@ -1,5 +1,6 @@
 import openmdao.api as om
 import dymos as dm
+import numpy as np
 
 import unittest
 from openmdao.utils.assert_utils import assert_near_equal
@@ -56,6 +57,9 @@ class TestCannonballMatrixState(unittest.TestCase):
 
         p.set_val('traj.phase.t_initial', 0)
         p.set_val('traj.phase.t_duration', 5)
+        if isinstance(tx, dm.Birkhoff):
+            p.set_val('traj.phase.initial_states:z', np.array([[0, 0], [10, 10]]))
+            p.set_val('traj.phase.final_states:z', np.array([[10, 0], [10, -10]]))
         p.set_val('traj.phase.states:z', phase.interp('z', [[[0, 0], [10, 10]], [[10, 0], [10, -10]]]))
 
         return p
@@ -99,6 +103,42 @@ class TestCannonballMatrixState(unittest.TestCase):
         assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
         assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
 
+    @require_pyoptsparse(optimizer='IPOPT')
+    def test_cannonball_matrix_state_birkhoff_lgl(self):
+        tx = dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=1, nodes_per_seg=20, grid_type='lgl'))
+
+        p = self._make_problem(tx)
+
+        dm.run_problem(p, simulate=True)
+
+        assert_near_equal(p.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
+        assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
+        assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
+
+        c = om.CaseReader('dymos_simulation.db').get_case('final')
+
+        assert_near_equal(c.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
+        assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
+        assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
+
+    @require_pyoptsparse(optimizer='IPOPT')
+    def test_cannonball_matrix_state_birkhoff_cgl(self):
+        tx = dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=1, nodes_per_seg=20, grid_type='cgl'))
+
+        p = self._make_problem(tx)
+
+        dm.run_problem(p, simulate=True)
+
+        assert_near_equal(p.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
+        assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
+        assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
+
+        c = om.CaseReader('dymos_simulation.db').get_case('final')
+
+        assert_near_equal(c.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
+        assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
+        assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
+
     @require_pyoptsparse(optimizer='IPOPT')
     def test_cannonball_matrix_state_radau_solve_segments(self):
 
diff --git a/dymos/examples/cannonball/test/test_two_phase_cannonball_birkhoff.py b/dymos/examples/cannonball/test/test_two_phase_cannonball_birkhoff.py
new file mode 100644
index 00000000..17ecf915
--- /dev/null
+++ b/dymos/examples/cannonball/test/test_two_phase_cannonball_birkhoff.py
@@ -0,0 +1,270 @@
+import unittest
+
+import matplotlib.pyplot as plt
+
+from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
+
+#
+# plt.switch_backend('Agg')
+import numpy as np
+from scipy.interpolate import interp1d
+
+import openmdao.api as om
+from openmdao.utils.assert_utils import assert_near_equal
+
+import dymos as dm
+from dymos.models.atmosphere.atmos_1976 import USatm1976Data
+
+
+class CannonballSizeComp(om.ExplicitComponent):
+    """
+    Compute the area and mass of a cannonball with a given radius and density.
+
+    Notes
+    -----
+    This component is not vectorized with 'num_nodes' as is the usual way
+    with Dymos, but is instead intended to compute a scalar mass and reference
+    area from scalar radius and density inputs. This component does not reside
+    in the ODE but instead its outputs are connected to the trajectory via
+    input design parameters.
+    """
+
+    def setup(self):
+        self.add_input(name='radius', val=1.0, desc='cannonball radius', units='m')
+        self.add_input(name='dens', val=7870., desc='cannonball density', units='kg/m**3')
+
+        self.add_output(name='mass', shape=(1,), desc='cannonball mass', units='kg')
+        self.add_output(name='S', shape=(1,), desc='aerodynamic reference area', units='m**2')
+
+        self.declare_partials(of='mass', wrt='dens')
+        self.declare_partials(of='mass', wrt='radius')
+
+        self.declare_partials(of='S', wrt='radius')
+
+    def compute(self, inputs, outputs):
+        radius = inputs['radius']
+        dens = inputs['dens']
+
+        outputs['mass'] = (4 / 3.) * dens * np.pi * radius ** 3
+        outputs['S'] = np.pi * radius ** 2
+
+    def compute_partials(self, inputs, partials):
+        radius = inputs['radius']
+        dens = inputs['dens']
+
+        partials['mass', 'dens'] = (4 / 3.) * np.pi * radius ** 3
+        partials['mass', 'radius'] = 4. * dens * np.pi * radius ** 2
+
+        partials['S', 'radius'] = 2 * np.pi * radius
+
+
+class CannonballODE(om.ExplicitComponent):
+    """
+    Cannonball ODE assuming flat earth and accounting for air resistance
+    """
+
+    def initialize(self):
+        self.options.declare('num_nodes', types=int)
+
+    def setup(self):
+        nn = self.options['num_nodes']
+
+        # static parameters
+        self.add_input('m', units='kg')
+        self.add_input('S', units='m**2')
+        # 0.5 good assumption for a sphere
+        self.add_input('CD', 0.5)
+
+        # 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', tags=['dymos.state_rate_source:v'])
+        self.add_output('gam_dot', shape=nn, units='rad/s', tags=['dymos.state_rate_source:gam'])
+        self.add_output('h_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:h'])
+        self.add_output('r_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:r'])
+        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, and use
+        # a graph coloring algorithm to automatically detect the sparsity pattern.
+        self.declare_coloring(wrt='*', method='cs')
+
+        alt_data = USatm1976Data.alt * om.unit_conversion('ft', 'm')[0]
+        rho_data = USatm1976Data.rho * om.unit_conversion('slug/ft**3', 'kg/m**3')[0]
+        self.rho_interp = interp1d(np.array(alt_data, dtype=complex),
+                                   np.array(rho_data, dtype=complex),
+                                   kind='linear')
+
+    def compute(self, inputs, outputs):
+
+        gam = inputs['gam']
+        v = inputs['v']
+        h = inputs['h']
+        m = inputs['m']
+        S = inputs['S']
+        CD = inputs['CD']
+
+        GRAVITY = 9.80665  # m/s**2
+
+        # handle complex-step gracefully from the interpolant
+        if np.iscomplexobj(h):
+            rho = self.rho_interp(inputs['h'])
+        else:
+            rho = self.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 TestTwoPhaseCannonballForDocs(unittest.TestCase):
+
+    @require_pyoptsparse(optimizer='SLSQP')
+    def make_problem(self, connected=False):
+
+        #############################################
+        # Setup the Dymos problem
+        #############################################
+
+        p = om.Problem(model=om.Group())
+
+        p.driver = om.pyOptSparseDriver()
+        p.driver.options['optimizer'] = 'SLSQP'
+        p.driver.declare_coloring()
+
+        p.model.add_subsystem('size_comp', CannonballSizeComp(),
+                              promotes_inputs=['radius', 'dens'])
+        p.model.set_input_defaults('dens', val=7.87, units='g/cm**3')
+        p.model.add_design_var('radius', lower=0.01, upper=0.10,
+                               ref0=0.01, ref=0.10, units='m')
+
+        traj = p.model.add_subsystem('traj', dm.Trajectory())
+
+        transcription = dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=1, nodes_per_seg=10))
+        ascent = dm.Phase(ode_class=CannonballODE, transcription=transcription)
+
+        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).
+        # The output of the ODE which provides the rate source for each state
+        # is obtained from the tags used on those outputs in the ODE.
+        # The units of the states are automatically inferred by multiplying the units
+        # of those rates by the time units.
+        ascent.set_time_options(fix_initial=True, duration_bounds=(1, 100),
+                                duration_ref=100, units='s')
+        ascent.set_state_options('r', fix_initial=True, fix_final=False)
+        ascent.set_state_options('h', fix_initial=True, fix_final=False)
+        ascent.set_state_options('gam', fix_initial=False, fix_final=True)
+        ascent.set_state_options('v', fix_initial=False, fix_final=False)
+
+        ascent.add_parameter('S', units='m**2', static_target=True)
+        ascent.add_parameter('m', units='kg', static_target=True)
+
+        # Limit the muzzle energy
+        ascent.add_boundary_constraint('ke', loc='initial',
+                                       upper=400000, lower=0, ref=100000)
+
+        # Second Phase (descent)
+        transcription = dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=1, nodes_per_seg=10))
+        descent = dm.Phase(ode_class=CannonballODE, transcription=transcription)
+
+        traj.add_phase('descent', descent)
+
+        # All initial states and time are free, since
+        #    they will be linked to the final states of ascent.
+        # Final altitude is fixed, because we will set
+        #    it to zero so that the phase ends at ground impact)
+        descent.set_time_options(duration_bounds=(.5, 100),
+                                 duration_ref=100, units='s')
+        descent.add_state('r')
+        descent.add_state('h', fix_initial=False, fix_final=True)
+        descent.add_state('gam', fix_initial=False, fix_final=False)
+        descent.add_state('v', fix_initial=False, fix_final=False)
+
+        descent.add_parameter('S', units='m**2', static_target=True)
+        descent.add_parameter('m', units='kg', static_target=True)
+
+        descent.add_objective('r', loc='final', scaler=-1.0)
+
+        # Add internally-managed design parameters to the trajectory.
+        traj.add_parameter('CD',
+                           targets={'ascent': ['CD'], 'descent': ['CD']},
+                           val=0.5, units=None, opt=False, static_target=True)
+
+        # Add externally-provided design parameters to the trajectory.
+        # In this case, we connect 'm' to pre-existing input parameters
+        # named 'm' in each phase.
+        traj.add_parameter('m', units='kg', val=1.0,
+                           targets={'ascent': 'm', 'descent': 'm'}, 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=['*'], connected=connected)
+
+        # Issue Connections
+        p.model.connect('size_comp.mass', 'traj.parameters:m')
+        p.model.connect('size_comp.S', 'traj.parameters:S')
+
+        # A linear solver at the top level can improve performance.
+        p.model.linear_solver = om.DirectSolver()
+
+        # Finish Problem Setup
+        p.setup()
+
+        #############################################
+        # Set constants and initial guesses
+        #############################################
+        p.set_val('radius', 0.05, units='m')
+        p.set_val('dens', 7.87, units='g/cm**3')
+
+        p.set_val('traj.parameters: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.initial_states:r', 0.0)
+        p.set_val('traj.ascent.initial_states:h', 0.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.ascent.final_states:gam', 0.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')
+        p.set_val('traj.descent.final_states:h', 0.0)
+
+        dm.run_problem(p)
+        assert_near_equal(p.get_val('traj.descent.states:r')[-1],
+                          3183.25, tolerance=1.0E-4)
+
+    def test_two_phase_cannonball_birkhoff_unconnected(self):
+        self.make_problem(connected=False)
+
+    def test_two_phase_cannonball_birkhoff_connected(self):
+        self.make_problem(connected=True)
+
+
+if __name__ == '__main__':  # pragma: no cover
+    unittest.main()
diff --git a/dymos/examples/double_integrator/test/test_double_integrator.py b/dymos/examples/double_integrator/test/test_double_integrator.py
index 1d8792cf..07a7f78c 100644
--- a/dymos/examples/double_integrator/test/test_double_integrator.py
+++ b/dymos/examples/double_integrator/test/test_double_integrator.py
@@ -10,17 +10,26 @@ import dymos as dm
 from dymos.examples.double_integrator.double_integrator_ode import DoubleIntegratorODE
 
 
-@require_pyoptsparse(optimizer='SLSQP')
-def double_integrator_direct_collocation(transcription='gauss-lobatto', compressed=True):
+@require_pyoptsparse(optimizer='IPOPT')
+def double_integrator(transcription='gauss-lobatto', compressed=True, grid_type='lgl', optimizer='IPOPT'):
 
     p = om.Problem(model=om.Group())
-    p.driver = om.pyOptSparseDriver()
+    p.driver = om.pyOptSparseDriver(optimizer=optimizer)
     p.driver.declare_coloring()
 
+    if optimizer == 'IPOPT':
+        p.driver.opt_settings['max_iter'] = 500
+        p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'
+        p.driver.opt_settings['print_level'] = 0
+        p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'
+        p.driver.opt_settings['tol'] = 1.0E-7
+
     if transcription == 'gauss-lobatto':
         t = dm.GaussLobatto(num_segments=30, order=3, compressed=compressed)
     elif transcription == "radau-ps":
         t = dm.Radau(num_segments=30, order=3, compressed=compressed)
+    elif transcription == 'birkhoff':
+        t = dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=1, nodes_per_seg=100, grid_type=grid_type))
     else:
         raise ValueError('invalid transcription')
 
@@ -33,7 +42,7 @@ def double_integrator_direct_collocation(transcription='gauss-lobatto', compress
     phase.add_state('v', fix_initial=True, fix_final=True, rate_source='u', units='m/s')
     phase.add_state('x', fix_initial=True, rate_source='v', units='m', shape=(1, ))
 
-    phase.add_control('u', units='m/s**2', scaler=0.01, continuity=False, rate_continuity=False,
+    phase.add_control('u', units='m/s**2', scaler=1.0, continuity=False, rate_continuity=False,
                       rate2_continuity=False, shape=(1, ), lower=-1.0, upper=1.0)
 
     # Maximize distance travelled in one second.
@@ -50,7 +59,13 @@ def double_integrator_direct_collocation(transcription='gauss-lobatto', compress
     p['traj.phase0.states:v'] = phase.interp('v', [0, 0])
     p['traj.phase0.controls:u'] = phase.interp('u', [1, -1])
 
-    p.run_driver()
+    if transcription == 'birkhoff':
+        p['traj.phase0.initial_states:x'] = 0.0
+        p['traj.phase0.initial_states:v'] = 0.0
+        p['traj.phase0.final_states:v'] = 0.0
+
+    simulate_kwargs = {'times_per_seg': 100 if transcription == 'birkhoff' else 10}
+    dm.run_problem(p, simulate=False, make_plots=True, simulate_kwargs=simulate_kwargs)
 
     return p
 
@@ -79,56 +94,40 @@ class TestDoubleIntegratorExample(unittest.TestCase):
         assert_near_equal(v[-1], 0.0, tolerance=tol)
 
     def test_ex_double_integrator_gl_compressed(self):
-        p = double_integrator_direct_collocation('gauss-lobatto',
-                                                 compressed=True)
-
-        x = p.get_val('traj.phase0.timeseries.x')
-        v = p.get_val('traj.phase0.timeseries.v')
-
-        assert_near_equal(x[0], 0.0, tolerance=1.0E-4)
-        assert_near_equal(x[-1], 0.25, tolerance=1.0E-4)
+        p = double_integrator('gauss-lobatto',
+                              compressed=True)
 
-        assert_near_equal(v[0], 0.0, tolerance=1.0E-4)
-        assert_near_equal(v[-1], 0.0, tolerance=1.0E-4)
+        self._assert_results(p)
 
     def test_ex_double_integrator_gl_uncompressed(self):
-        p = double_integrator_direct_collocation('gauss-lobatto',
-                                                 compressed=False)
+        p = double_integrator('gauss-lobatto',
+                              compressed=False)
 
-        x = p.get_val('traj.phase0.timeseries.x')
-        v = p.get_val('traj.phase0.timeseries.v')
-
-        assert_near_equal(x[0], 0.0, tolerance=1.0E-4)
-        assert_near_equal(x[-1], 0.25, tolerance=1.0E-4)
-
-        assert_near_equal(v[0], 0.0, tolerance=1.0E-4)
-        assert_near_equal(v[-1], 0.0, tolerance=1.0E-4)
+        self._assert_results(p)
 
     def test_ex_double_integrator_radau_compressed(self):
-        p = double_integrator_direct_collocation('radau-ps',
-                                                 compressed=True)
+        p = double_integrator('radau-ps',
+                              compressed=True)
 
-        x = p.get_val('traj.phase0.timeseries.x')
-        v = p.get_val('traj.phase0.timeseries.v')
+        self._assert_results(p)
 
-        assert_near_equal(x[0], 0.0, tolerance=1.0E-4)
-        assert_near_equal(x[-1], 0.25, tolerance=1.0E-4)
+    def test_ex_double_integrator_radau_uncompressed(self):
+        p = double_integrator('radau-ps',
+                              compressed=False)
 
-        assert_near_equal(v[0], 0.0, tolerance=1.0E-4)
-        assert_near_equal(v[-1], 0.0, tolerance=1.0E-4)
+        self._assert_results(p)
 
-    def test_ex_double_integrator_radau_uncompressed(self):
-        p = double_integrator_direct_collocation('radau-ps',
-                                                 compressed=False)
+    def test_ex_double_integrator_birkhoff_lgl(self):
+        p = double_integrator('birkhoff',
+                              grid_type='lgl')
 
-        x = p.get_val('traj.phase0.timeseries.x')
-        v = p.get_val('traj.phase0.timeseries.v')
+        self._assert_results(p)
 
-        assert_near_equal(x[0], 0.0, tolerance=1.0E-4)
-        assert_near_equal(x[-1], 0.25, tolerance=1.0E-4)
+    def test_ex_double_integrator_birkhoff_cgl(self):
+        p = double_integrator('birkhoff',
+                              grid_type='cgl')
 
-        assert_near_equal(v[0], 0.0, tolerance=1.0E-4)
-        assert_near_equal(v[-1], 0.0, tolerance=1.0E-4)
+        self._assert_results(p)
 
     @require_pyoptsparse(optimizer='SLSQP')
     def test_ex_double_integrator_input_times_uncompressed(self):
diff --git a/dymos/examples/double_integrator/test/test_double_integrator_breakwell.py b/dymos/examples/double_integrator/test/test_double_integrator_breakwell.py
new file mode 100644
index 00000000..c46a1b58
--- /dev/null
+++ b/dymos/examples/double_integrator/test/test_double_integrator_breakwell.py
@@ -0,0 +1,140 @@
+import itertools
+import os
+import unittest
+
+import numpy as np
+import openmdao.api as om
+from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials
+from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
+
+import dymos as dm
+from dymos.examples.double_integrator.breakwell_ode import BreakwellODE
+
+
+@require_pyoptsparse(optimizer='IPOPT')
+def double_integrator_direct_collocation(transcription='gauss-lobatto', compressed=True, optimizer='IPOPT'):
+
+    p = om.Problem(model=om.Group())
+    p.driver = om.pyOptSparseDriver()
+    p.driver.declare_coloring()
+    p.driver.options['optimizer'] = optimizer
+
+    if optimizer == 'IPOPT':
+        p.driver.opt_settings['max_iter'] = 5000
+        p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'
+        p.driver.opt_settings['print_level'] = 5
+        p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'
+        p.driver.opt_settings['tol'] = 1.0E-3
+        p.driver.opt_settings['constr_viol_tol'] = 1.0E-6
+    elif optimizer == 'SNOPT':
+        p.driver.opt_settings['Major iterations limit'] = 1000
+        p.driver.opt_settings['Major feasibility tolerance'] = 1.0E-6
+        p.driver.opt_settings['Major optimality tolerance'] = 1.0E-3
+        p.driver.opt_settings['iSumm'] = 6
+
+    if transcription == 'gauss-lobatto':
+        t = dm.GaussLobatto(num_segments=30, order=3, compressed=compressed)
+    elif transcription == "radau-ps":
+        t = dm.Radau(num_segments=30, order=3, compressed=compressed)
+    elif transcription == 'birkhoff':
+        grid = dm.BirkhoffGrid(num_segments=1, nodes_per_seg=51)
+        t = dm.Birkhoff(grid=grid)
+    else:
+        raise ValueError('invalid transcription')
+
+    traj = p.model.add_subsystem('traj', dm.Trajectory())
+
+    phase = traj.add_phase('phase0', dm.Phase(ode_class=BreakwellODE, transcription=t))
+
+    phase.set_time_options(fix_initial=True, fix_duration=True, units='s')
+
+    phase.add_state('v', fix_initial=True, rate_source='u', units='m/s')
+    phase.add_state('x', fix_initial=True, rate_source='v', units='m', shape=(1, ))
+    phase.add_state('J', fix_initial=True, rate_source='J_dot')
+
+    phase.add_control('u', units='m/s**2', scaler=0.1, continuity=False, rate_continuity=False,
+                      rate2_continuity=False, shape=(1, ))
+
+    # Maximize distance travelled in one second.
+    phase.add_objective('J', loc='final')
+
+    phase.add_path_constraint('x', upper=0.1)
+    phase.add_boundary_constraint('x', loc='final', equals=0.0)
+    phase.add_boundary_constraint('v', loc='final', equals=-1.0)
+
+    p.model.linear_solver = om.DirectSolver()
+
+    p.setup(check=True, force_alloc_complex=True)
+
+    p['traj.phase0.t_initial'] = 0.0
+    p['traj.phase0.t_duration'] = 1.0
+    p['traj.phase0.states:x'] = phase.interp('x', [0, 0])
+    p['traj.phase0.states:v'] = phase.interp('v', [1, 0])
+    p['traj.phase0.states:J'] = phase.interp('J', [0, 4])
+    p['traj.phase0.controls:u'] = phase.interp('u', [-5, -5])
+    if transcription == 'birkhoff':
+        p['traj.phase0.initial_states:x'] = 0.0
+        p['traj.phase0.initial_states:v'] = 1.0
+        p['traj.phase0.initial_states:J'] = 0.0
+
+    return p
+
+
+@use_tempdirs
+class TestBreakwellExample(unittest.TestCase):
+
+    @classmethod
+    def tearDownClass(cls):
+        for filename in ['total_coloring.pkl', 'SLSQP.out', 'SNOPT_print.out']:
+            if os.path.exists(filename):
+                os.remove(filename)
+
+    @staticmethod
+    def _assert_results(p, tol=1.0E-4):
+        x = p.get_val('traj.phase0.timeseries.x')
+        v = p.get_val('traj.phase0.timeseries.v')
+        J = p.get_val('traj.phase0.timeseries.J')
+        J_expected = 2*v[0]**3/(9*(0.1 - x[0])) - 2*v[-1]**3/(9*(0.1-x[-1]))
+
+        assert_near_equal(x[-1], 0.0, tolerance=tol)
+        assert_near_equal(v[-1], -1.0, tolerance=tol)
+        assert_near_equal(J[-1], J_expected, tolerance=tol)
+
+    def test_ex_double_integrator_gl_compressed(self):
+        p = double_integrator_direct_collocation('gauss-lobatto',
+                                                 compressed=True)
+        dm.run_problem(p)
+        self._assert_results(p)
+
+    def test_ex_double_integrator_gl_uncompressed(self):
+        p = double_integrator_direct_collocation('gauss-lobatto',
+                                                 compressed=False)
+        dm.run_problem(p)
+        self._assert_results(p)
+
+    def test_ex_double_integrator_radau_compressed(self):
+        p = double_integrator_direct_collocation('radau-ps',
+                                                 compressed=True, optimizer='IPOPT')
+        dm.run_problem(p)
+        self._assert_results(p)
+
+    def test_ex_double_integrator_radau_uncompressed(self):
+        p = double_integrator_direct_collocation('radau-ps',
+                                                 compressed=False)
+        dm.run_problem(p)
+        self._assert_results(p)
+
+    def test_ex_double_integrator_birkhoff(self):
+        p = double_integrator_direct_collocation('birkhoff', optimizer='IPOPT')
+        dm.run_problem(p)
+        self._assert_results(p)
+
+    def test_check_partials(self):
+        p = double_integrator_direct_collocation('birkhoff', optimizer='SLSQP')
+        p.run_model()
+        cpd = p.check_partials(method='cs', out_stream=None)
+        assert_check_partials(cpd)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/dymos/examples/double_integrator/test/test_double_integrator_timeseries_units.py b/dymos/examples/double_integrator/test/test_double_integrator_timeseries_units.py
index 488a898f..37699105 100644
--- a/dymos/examples/double_integrator/test/test_double_integrator_timeseries_units.py
+++ b/dymos/examples/double_integrator/test/test_double_integrator_timeseries_units.py
@@ -55,7 +55,7 @@ def double_integrator_direct_collocation(transcription=dm.GaussLobatto, compress
     return p
 
 
-# @use_tempdirs
+@use_tempdirs
 class TestDoubleIntegratorExample(unittest.TestCase):
 
     @classmethod
diff --git a/dymos/examples/goddard_rocket_problem/test/__init__.py b/dymos/examples/goddard_rocket_problem/test/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/dymos/examples/goddard_rocket_problem/test/test_goddard_rocket_problem.py b/dymos/examples/goddard_rocket_problem/test/test_goddard_rocket_problem.py
new file mode 100644
index 00000000..ef627e0d
--- /dev/null
+++ b/dymos/examples/goddard_rocket_problem/test/test_goddard_rocket_problem.py
@@ -0,0 +1,138 @@
+import openmdao.api as om
+import dymos as dm
+import unittest
+import os
+
+from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials
+from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
+
+from dymos.examples.goddard_rocket_problem import RocketODE
+
+SHOW_PLOTS = True
+
+
+@require_pyoptsparse(optimizer='IPOPT')
+def goddard_rocket_direct_collocation(grid_type='lgl'):
+
+    optimizer = 'IPOPT'
+    p = om.Problem(model=om.Group())
+    p.driver = om.pyOptSparseDriver()
+    p.driver.options['optimizer'] = optimizer
+    p.driver.declare_coloring()
+    if optimizer == 'SNOPT':
+        p.driver.opt_settings['Major iterations limit'] = 1000
+        p.driver.opt_settings['Major feasibility tolerance'] = 1.0E-6
+        p.driver.opt_settings['Major optimality tolerance'] = 1.0E-3
+        p.driver.opt_settings['iSumm'] = 6
+
+    t = dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=1, nodes_per_seg=150, grid_type=grid_type))
+
+    traj = p.model.add_subsystem('traj', dm.Trajectory())
+
+    phase = traj.add_phase('phase0', dm.Phase(ode_class=RocketODE, transcription=t))
+
+    phase.set_time_options(fix_initial=True, fix_duration=False, duration_bounds=(20, 80), units='s',
+                           duration_ref=100)
+
+    phase.add_state('h', fix_initial=True, fix_final=False, rate_source='h_dot', units='ft', lower=0, ref=20000)
+    phase.add_state('v', fix_initial=True, fix_final=False, rate_source='v_dot', units='ft/s', ref=1000, lower=0)
+    phase.add_state('m', fix_initial=True, fix_final=True, rate_source='m_dot', units='slug',
+                    lower=1e-3, ref=3, upper=3.0)
+
+    phase.add_control('T', units='lbf', lower=0, upper=200, ref=200)
+
+    # Minimize the control effort
+    phase.add_objective('h', loc='final', scaler=-1)
+
+    p.model.linear_solver = om.DirectSolver()
+
+    p.setup(check=True, force_alloc_complex=True)
+
+    p['traj.phase0.t_initial'] = 0.0
+    p['traj.phase0.t_duration'] = 40.0
+
+    p['traj.phase0.initial_states:h'] = 0.0
+    p['traj.phase0.initial_states:v'] = 0.0
+    p['traj.phase0.initial_states:m'] = 3.0
+    p['traj.phase0.final_states:h'] = 10000.0
+    p['traj.phase0.final_states:v'] = 0.0
+    p['traj.phase0.final_states:m'] = 1.0
+
+    p['traj.phase0.states:h'] = phase.interp('h', [0.0, 10000.0])
+    p['traj.phase0.states:v'] = phase.interp('v', [0.0, 0.0])
+    p['traj.phase0.states:m'] = phase.interp('m', [3.0, 3.0])
+    p['traj.phase0.controls:T'] = phase.interp('T', [0.0, 0.0])
+
+    return p
+
+
+@use_tempdirs
+class TestGoddardRocketProblem(unittest.TestCase):
+
+    @classmethod
+    def tearDownClass(cls):
+        for filename in ['total_coloring.pkl', 'SLSQP.out', 'SNOPT_print.out']:
+            if os.path.exists(filename):
+                os.remove(filename)
+
+    @staticmethod
+    def _assert_results(p, tol=1.0E-4):
+        h = p.get_val('traj.phase0.timeseries.h')
+        assert_near_equal(h[-1], 18745.0, tolerance=tol)
+
+    def test_goddard_rocket_problem_lgl(self):
+        p = goddard_rocket_direct_collocation(grid_type='lgl')
+        dm.run_problem(p)
+        if SHOW_PLOTS:
+            t = p.get_val('traj.phase0.timeseries.time')
+            h = p.get_val('traj.phase0.timeseries.h')
+            v = p.get_val('traj.phase0.timeseries.v')
+            m = p.get_val('traj.phase0.timeseries.m')
+            T = p.get_val('traj.phase0.timeseries.T')
+
+            import matplotlib.pyplot as plt
+            fig, axs = plt.subplots(2, 2)
+            axs[0, 0].plot(t, h)
+            axs[0, 0].set_xlabel('time (s)')
+            axs[0, 0].set_ylabel('altitude (ft)')
+            axs[0, 1].plot(t, v)
+            axs[0, 1].set_xlabel('time (s)')
+            axs[0, 1].set_ylabel('velocity (ft/s)')
+            axs[1, 0].plot(t, m)
+            axs[1, 0].set_xlabel('time (s)')
+            axs[1, 0].set_ylabel('mass (slug)')
+            axs[1, 1].plot(t, T)
+            axs[1, 1].set_xlabel('time (s)')
+            axs[1, 1].set_ylabel('thrust (lb)')
+            plt.show()
+
+        self._assert_results(p)
+
+    def test_goddard_rocket_problem_cgl(self):
+        p = goddard_rocket_direct_collocation(grid_type='cgl')
+        dm.run_problem(p)
+        if SHOW_PLOTS:
+            t = p.get_val('traj.phase0.timeseries.time')
+            h = p.get_val('traj.phase0.timeseries.h')
+            v = p.get_val('traj.phase0.timeseries.v')
+            m = p.get_val('traj.phase0.timeseries.m')
+            T = p.get_val('traj.phase0.timeseries.T')
+
+            import matplotlib.pyplot as plt
+            fig, axs = plt.subplots(2, 2)
+            axs[0, 0].plot(t, h)
+            axs[0, 1].plot(t, v)
+            axs[1, 0].plot(t, m)
+            axs[1, 1].plot(t, T)
+            plt.show()
+        self._assert_results(p)
+
+    def test_check_partials(self):
+        p = goddard_rocket_direct_collocation(grid_type='cgl')
+        p.run_model()
+        cpd = p.check_partials(method='cs', out_stream=None)
+        assert_check_partials(cpd)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/dymos/examples/hyper_sensitive/test/test_hyper_sensitive.py b/dymos/examples/hyper_sensitive/test/test_hyper_sensitive.py
index 579bde36..5ed7dbe8 100644
--- a/dymos/examples/hyper_sensitive/test/test_hyper_sensitive.py
+++ b/dymos/examples/hyper_sensitive/test/test_hyper_sensitive.py
@@ -5,14 +5,13 @@ import unittest
 import warnings
 
 import numpy as np
+import openmdao.api as om
 
-from openmdao.api import Problem, Group, pyOptSparseDriver
 from openmdao.utils.assert_utils import assert_near_equal
 from openmdao.utils.general_utils import printoptions
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 
 import dymos as dm
-from dymos import Trajectory, GaussLobatto, Phase, Radau
 from dymos.examples.hyper_sensitive.hyper_sensitive_ode import HyperSensitiveODE
 
 
@@ -28,29 +27,38 @@ class TestHyperSensitive(unittest.TestCase):
                 os.remove(filename)
 
     @require_pyoptsparse(optimizer='SLSQP')
-    def make_problem(self, transcription=GaussLobatto, optimizer='SLSQP', numseg=30):
-        p = Problem(model=Group())
-        p.driver = pyOptSparseDriver()
+    def make_problem(self, transcription=dm.GaussLobatto, optimizer='SLSQP', numseg=30, order=3):
+        p = om.Problem(model=om.Group())
+        p.driver = om.pyOptSparseDriver()
         p.driver.declare_coloring()
         p.driver.options['optimizer'] = optimizer
 
         if optimizer == 'SNOPT':
             p.driver.declare_coloring()
-            p.driver.opt_settings['Major iterations limit'] = 100
+            p.driver.opt_settings['iSumm'] = 6
+            p.driver.opt_settings['Major iterations limit'] = 500
             p.driver.opt_settings['Major feasibility tolerance'] = 1.0E-6
             p.driver.opt_settings['Major optimality tolerance'] = 1.0E-6
         elif optimizer == 'IPOPT':
-            p.driver.opt_settings['hessian_approximation'] = 'limited-memory'
-            # p.driver.opt_settings['nlp_scaling_method'] = 'user-scaling'
-            p.driver.opt_settings['print_level'] = 0
-            p.driver.opt_settings['linear_solver'] = 'mumps'
+            p.driver.opt_settings['print_level'] = 5
+            p.driver.opt_settings['mu_strategy'] = 'adaptive'
+            p.driver.opt_settings['bound_mult_init_method'] = 'mu-based'
+            p.driver.opt_settings['mu_init'] = 0.01
+            p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'
+            p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'
             p.driver.declare_coloring()
         else:
             p.driver.declare_coloring()
 
-        traj = p.model.add_subsystem('traj', Trajectory())
-        phase0 = traj.add_phase('phase0', Phase(ode_class=HyperSensitiveODE,
-                                                transcription=transcription(num_segments=numseg, order=3)))
+        if transcription == 'gauss-lobatto':
+            t = dm.GaussLobatto(num_segments=numseg, order=3)
+        elif transcription == 'radau-ps':
+            t = dm.Radau(num_segments=numseg, order=3)
+        elif transcription == 'birkhoff':
+            t = dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=5, nodes_per_seg=51, grid_type='lgl'))
+
+        traj = p.model.add_subsystem('traj', dm.Trajectory())
+        phase0 = traj.add_phase('phase0', dm.Phase(ode_class=HyperSensitiveODE, transcription=t))
         phase0.set_time_options(fix_initial=True, fix_duration=True)
         phase0.add_state('x', fix_initial=True, fix_final=False, rate_source='x_dot')
         phase0.add_state('xL', fix_initial=True, fix_final=False, rate_source='L')
@@ -64,15 +72,26 @@ class TestHyperSensitive(unittest.TestCase):
 
         p.setup(check=True)
 
-        p.set_val('traj.phase0.states:x', phase0.interp('x', [1.5, 1]))
-        p.set_val('traj.phase0.states:xL', phase0.interp('xL', [0, 1]))
+        if transcription == 'birkhoff':
+            p.set_val('traj.phase0.initial_states:x', 1.5)
+            p.set_val('traj.phase0.initial_states:xL', 0.0)
+        else:
+            p.set_val('traj.phase0.states:x', phase0.interp('x', [1.5, 1]))
+            p.set_val('traj.phase0.states:xL', phase0.interp('xL', [0, 1]))
         p.set_val('traj.phase0.t_initial', 0)
         p.set_val('traj.phase0.t_duration', tf)
         p.set_val('traj.phase0.controls:u', phase0.interp('u', [-0.6, 2.4]))
 
+        if transcription == dm.Birkhoff:
+            p.set_val('traj.phase0.initial_states:x', 1.5)
+            p.set_val('traj.phase0.final_states:x', 1.0)
+            p.set_val('traj.phase0.initial_states:xL', 0.0)
+            p.set_val('traj.phase0.final_states:xL', 1.0)
+
         return p
 
-    def solution(self):
+    @staticmethod
+    def solution():
         sqrt_two = np.sqrt(2)
         val = sqrt_two * tf
         c1 = (1.5 * np.exp(-val) - 1) / (np.exp(-val) - np.exp(val))
@@ -85,14 +104,14 @@ class TestHyperSensitive(unittest.TestCase):
         return ui, uf, J
 
     def test_partials(self):
-        p = self.make_problem(transcription=Radau, optimizer='SLSQP')
+        p = self.make_problem(transcription='radau-ps', optimizer='SLSQP')
         p.run_model()
         with printoptions(linewidth=1024, edgeitems=100):
             cpd = p.check_partials(method='fd', compact_print=True, out_stream=None)
 
     @require_pyoptsparse(optimizer='IPOPT')
     def test_hyper_sensitive_radau(self):
-        p = self.make_problem(transcription=Radau, optimizer='IPOPT')
+        p = self.make_problem(transcription='radau-ps', optimizer='IPOPT')
         dm.run_problem(p, refine_iteration_limit=5)
         ui, uf, J = self.solution()
 
@@ -108,9 +127,27 @@ class TestHyperSensitive(unittest.TestCase):
                           J,
                           tolerance=1e-6)
 
+    @require_pyoptsparse(optimizer='IPOPT')
+    def test_hyper_sensitive_birkhoff(self):
+        p = self.make_problem(transcription='birkhoff', optimizer='IPOPT')
+        dm.run_problem(p, make_plots=True)
+        ui, uf, J = self.solution()
+
+        assert_near_equal(p.get_val('traj.phase0.timeseries.u')[0],
+                          ui,
+                          tolerance=1e-3)
+
+        assert_near_equal(p.get_val('traj.phase0.timeseries.u')[-1],
+                          uf,
+                          tolerance=1e-3)
+
+        assert_near_equal(p.get_val('traj.phase0.timeseries.xL')[-1],
+                          J,
+                          tolerance=1e-6)
+
     @require_pyoptsparse(optimizer='IPOPT')
     def test_hyper_sensitive_gauss_lobatto(self):
-        p = self.make_problem(transcription=GaussLobatto, optimizer='IPOPT')
+        p = self.make_problem(transcription='gauss-lobatto', optimizer='IPOPT')
         dm.run_problem(p, refine_iteration_limit=5)
 
         ui, uf, J = self.solution()
@@ -129,7 +166,7 @@ class TestHyperSensitive(unittest.TestCase):
 
     @require_pyoptsparse(optimizer='IPOPT')
     def test_refinement_warning(self):
-        p = self.make_problem(transcription=Radau, optimizer='IPOPT')
+        p = self.make_problem(transcription='radau-ps', optimizer='IPOPT')
 
         msg = "Refinement not performed. Set run_driver to True to perform refinement."
 
diff --git a/dymos/examples/low_thrust_spiral/test/__init__.py b/dymos/examples/low_thrust_spiral/test/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/dymos/examples/low_thrust_spiral/test/test_low_thrust_spiral.py b/dymos/examples/low_thrust_spiral/test/test_low_thrust_spiral.py
new file mode 100644
index 00000000..3b8419dd
--- /dev/null
+++ b/dymos/examples/low_thrust_spiral/test/test_low_thrust_spiral.py
@@ -0,0 +1,123 @@
+import openmdao.api as om
+import dymos as dm
+import unittest
+import numpy as np
+
+from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials
+from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
+from dymos.examples.low_thrust_spiral import LowThrustODE
+
+show_plots = False
+
+
+@require_pyoptsparse(optimizer='SNOPT')
+def low_thrust_spiral_direct_collocation(grid_type='lgl'):
+
+    optimizer = 'SNOPT'
+    p = om.Problem(model=om.Group())
+    p.driver = om.pyOptSparseDriver()
+    p.driver.options['optimizer'] = optimizer
+    p.driver.declare_coloring()
+    if optimizer == 'SNOPT':
+        p.driver.opt_settings['Major iterations limit'] = 1000
+        p.driver.opt_settings['Major feasibility tolerance'] = 1.0E-6
+        p.driver.opt_settings['Major optimality tolerance'] = 5.0E-4
+        p.driver.opt_settings['iSumm'] = 6
+
+    t = dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=12, nodes_per_seg=50, grid_type=grid_type))
+
+    traj = p.model.add_subsystem('traj', dm.Trajectory())
+
+    phase = traj.add_phase('phase0', dm.Phase(ode_class=LowThrustODE, transcription=t))
+
+    phase.set_time_options(fix_initial=True, fix_duration=False, units='s', duration_bounds=(100, 500))
+
+    phase.add_state('r', fix_initial=True, fix_final=True, rate_source='vr', lower=0.5, upper=6.5,
+                    defect_ref=1e2)
+    phase.add_state('theta', fix_initial=True, fix_final=False, rate_source='theta_dot', units='rad', defect_ref=1e2)
+    phase.add_state('vr', fix_initial=True, fix_final=True, rate_source='vr_dot', units='1/s', defect_ref=1e2)
+    phase.add_state('vt', fix_initial=True, fix_final=True, rate_source='vt_dot', units='1/s', defect_ref=1e2)
+
+    phase.add_control('alpha', units='rad')
+
+    # Minimize the control effort
+    phase.add_objective('time', loc='final')
+
+    p.model.linear_solver = om.DirectSolver()
+
+    p.setup(check=True, force_alloc_complex=True)
+
+    p['traj.phase0.t_initial'] = 0.0
+    p['traj.phase0.t_duration'] = 300.0
+
+    p['traj.phase0.initial_states:r'] = 1.0
+    p['traj.phase0.initial_states:theta'] = 0.0
+    p['traj.phase0.initial_states:vr'] = 0.0
+    p['traj.phase0.initial_states:vt'] = 1.0
+    p['traj.phase0.final_states:r'] = 6.0
+    p['traj.phase0.final_states:vr'] = 0.0
+    p['traj.phase0.final_states:vt'] = 1/np.sqrt(6)
+
+    p['traj.phase0.states:r'] = phase.interp('r', [1.0, 3.0])
+    p['traj.phase0.states:vr'] = phase.interp('vr', [0.0, 0.0])
+    p['traj.phase0.states:vt'] = phase.interp('vt', [1.0, 1.0])
+    p['traj.phase0.controls:alpha'] = phase.interp('alpha', [0.0, 0.0])
+
+    return p
+
+
+@use_tempdirs
+class TestLowThrustSpiral(unittest.TestCase):
+
+    # @classmethod
+    # def tearDownClass(cls):
+    #     for filename in ['total_coloring.pkl', 'SLSQP.out', 'SNOPT_print.out']:
+    #         if os.path.exists(filename):
+    #             os.remove(filename)
+
+    @staticmethod
+    def _assert_results(p, tol=1.0E-4):
+        # t = p.get_val('traj.phase0.timeseries.time')
+        #
+        # assert_near_equal(t[-1], 228, tolerance=tol)
+        return
+
+    def test_low_thrust_spiral_lgl(self):
+        p = low_thrust_spiral_direct_collocation(grid_type='lgl')
+        dm.run_problem(p)
+        theta = p.get_val('traj.phase0.timeseries.theta')
+        r = p.get_val('traj.phase0.timeseries.r')
+
+        if show_plots:
+            import matplotlib.pyplot as plt
+            fig, ax = plt.subplots(subplot_kw={'projection': 'polar'})
+            ax.plot(theta, r)
+            ax.grid(True)
+            plt.show()
+
+        self._assert_results(p)
+
+    def test_low_thrust_spiral_cgl(self):
+        p = low_thrust_spiral_direct_collocation(grid_type='cgl')
+        dm.run_problem(p)
+        theta = p.get_val('traj.phase0.timeseries.theta')
+        r = p.get_val('traj.phase0.timeseries.r')
+
+        if show_plots:
+            import matplotlib.pyplot as plt
+            fig, ax = plt.subplots(subplot_kw={'projection': 'polar'})
+            ax.plot(theta, r)
+            ax.grid(True)
+            plt.show()
+
+        self._assert_results(p)
+
+    def test_check_partials(self):
+        p = low_thrust_spiral_direct_collocation(grid_type='cgl')
+        p.run_model()
+        cpd = p.check_partials(method='cs', out_stream=None)
+        assert_check_partials(cpd)
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/dymos/examples/moon_landing/test/__init__.py b/dymos/examples/moon_landing/test/__init__.py
new file mode 100644
index 00000000..e69de29b
diff --git a/dymos/examples/moon_landing/test/test_moon_landing_problem.py b/dymos/examples/moon_landing/test/test_moon_landing_problem.py
new file mode 100644
index 00000000..8930de0b
--- /dev/null
+++ b/dymos/examples/moon_landing/test/test_moon_landing_problem.py
@@ -0,0 +1,87 @@
+import numpy as np
+import openmdao.api as om
+import dymos as dm
+import unittest
+
+from openmdao.utils.assert_utils import assert_near_equal
+from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
+from dymos.examples.moon_landing import MoonLandingProblemODE
+
+
+@use_tempdirs
+class TestMoonLandingProblem(unittest.TestCase):
+
+    @require_pyoptsparse(optimizer='IPOPT')
+    def make_problem(self, grid_type):
+        self.p = om.Problem(model=om.Group())
+        self.p.driver = om.pyOptSparseDriver()
+        self.p.driver.declare_coloring()
+        self.p.driver.options['optimizer'] = 'IPOPT'
+        self.p.driver.opt_settings['hessian_approximation'] = 'limited-memory'
+        self.p.driver.opt_settings['print_level'] = 0
+        self.p.driver.opt_settings['linear_solver'] = 'mumps'
+        self.p.driver.declare_coloring()
+
+        t = dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=1, nodes_per_seg=20, grid_type=grid_type))
+
+        traj = self.p.model.add_subsystem('traj', dm.Trajectory())
+        phase = dm.Phase(ode_class=MoonLandingProblemODE,
+                         transcription=t)
+
+        phase.set_time_options(fix_initial=True, fix_duration=False)
+        phase.add_state('h', fix_initial=True, rate_source='h_dot')
+        phase.add_state('v', fix_initial=True, rate_source='v_dot')
+        phase.add_state('m', fix_initial=True, lower=1e-3, rate_source='m_dot')
+        phase.add_control('T', lower=0.0, upper=1.227)
+
+        phase.add_boundary_constraint('h', loc='final', equals=0.0)
+        phase.add_boundary_constraint('v', loc='final', equals=0.0)
+
+        phase.add_objective('m', scaler=-1)
+        phase.set_simulate_options(atol=1.0E-1, rtol=1.0E-2)
+
+        traj.add_phase('phase', phase)
+
+        self.p.setup(check=True, force_alloc_complex=True)
+
+        self.p.set_val('traj.phase.initial_states:h', 1.0)
+        self.p.set_val('traj.phase.initial_states:v', -0.783)
+        self.p.set_val('traj.phase.initial_states:m', 1.0)
+        self.p.set_val('traj.phase.t_initial', 0)
+        self.p.set_val('traj.phase.t_duration', 1.0)
+        self.p.set_val('traj.phase.controls:T', phase.interp('T', [0.0, 1.227]))
+
+    def test_problem_lgl(self):
+
+        with dm.options.temporary(include_check_partials=True):
+            self.make_problem(grid_type='lgl')
+            dm.run_problem(self.p, simulate=True, simulate_kwargs={'times_per_seg': 100},
+                           make_plots=True)
+
+            with np.printoptions(linewidth=1024):
+                self.p.check_partials(compact_print=True, method='cs')
+
+        h = self.p.get_val('traj.phase.timeseries.h')
+        v = self.p.get_val('traj.phase.timeseries.v')
+        m = self.p.get_val('traj.phase.timeseries.m')
+        T = self.p.get_val('traj.phase.timeseries.T')
+
+        assert_near_equal(T[0], 0.0, tolerance=1e-5)
+        assert_near_equal(T[-1], 1.227, tolerance=1e-5)
+        assert_near_equal(h[-1], 0.0, tolerance=1e-5)
+        assert_near_equal(v[-1], 0.0, tolerance=1e-5)
+        assert_near_equal(m[-1], 0.3953, tolerance=1e-3)
+
+    def test_problem_cgl(self):
+        self.make_problem(grid_type='cgl')
+        dm.run_problem(self.p, make_plots=True)
+        h = self.p.get_val('traj.phase.timeseries.h')
+        v = self.p.get_val('traj.phase.timeseries.v')
+        m = self.p.get_val('traj.phase.timeseries.m')
+        T = self.p.get_val('traj.phase.timeseries.T')
+
+        assert_near_equal(T[0], 0.0, tolerance=1e-5)
+        assert_near_equal(T[-1], 1.227, tolerance=1e-5)
+        assert_near_equal(h[-1], 0.0, tolerance=1e-5)
+        assert_near_equal(v[-1], 0.0, tolerance=1e-5)
+        assert_near_equal(m[-1], 0.3953, tolerance=1e-4)
diff --git a/dymos/examples/robot_arm/doc/test_doc_robot_arm.py b/dymos/examples/robot_arm/doc/test_doc_robot_arm.py
index 992f3474..5cca7b62 100644
--- a/dymos/examples/robot_arm/doc/test_doc_robot_arm.py
+++ b/dymos/examples/robot_arm/doc/test_doc_robot_arm.py
@@ -1,17 +1,14 @@
-from __future__ import print_function, division, absolute_import
-
 import os
 import unittest
 
 import numpy as np
+import openmdao.api as om
+import dymos as dm
 
-from openmdao.api import Problem, Group, pyOptSparseDriver
 from openmdao.utils.assert_utils import assert_near_equal
 from openmdao.utils.general_utils import printoptions
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 
-import dymos as dm
-from dymos import Trajectory, GaussLobatto, Phase, Radau
 from dymos.examples.robot_arm.robot_arm_ode import RobotArmODE
 
 
@@ -24,23 +21,21 @@ class TestRobotArm(unittest.TestCase):
                 os.remove(filename)
 
     @require_pyoptsparse(optimizer='SLSQP')
-    def make_problem(self, transcription=Radau, optimizer='SLSQP', numseg=30):
-        p = Problem(model=Group())
-        p.driver = pyOptSparseDriver()
-        p.driver.declare_coloring(tol=1.0E-12)
+    def make_problem(self, tx, optimizer='SLSQP', numseg=30):
+        p = om.Problem(model=om.Group())
+        p.driver = om.pyOptSparseDriver()
+        p.driver.declare_coloring()
         p.driver.options['optimizer'] = optimizer
         if optimizer == 'SNOPT':
             p.driver.opt_settings['iSumm'] = 6
             p.driver.opt_settings['Verify level'] = 3
         elif optimizer == 'IPOPT':
+            p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'
             p.driver.opt_settings['max_iter'] = 500
             p.driver.opt_settings['print_level'] = 0
-            p.driver.opt_settings['tol'] = 1.0E-6
-            p.driver.opt_settings['acceptable_tol'] = 1.0E-5
-        traj = p.model.add_subsystem('traj', Trajectory())
+        traj = p.model.add_subsystem('traj', dm.Trajectory())
 
-        phase = traj.add_phase('phase', Phase(ode_class=RobotArmODE,
-                                              transcription=transcription(num_segments=numseg, order=3)))
+        phase = traj.add_phase('phase', dm.Phase(ode_class=RobotArmODE, transcription=tx))
         phase.set_time_options(fix_initial=True, fix_duration=False)
 
         phase.add_state('x0', fix_initial=True, fix_final=True, rate_source='x0_dot', units='m')
@@ -69,24 +64,39 @@ class TestRobotArm(unittest.TestCase):
 
         p.set_val('traj.phase.t_initial', 0)
         p.set_val('traj.phase.t_duration', 10)
-        p.set_val('traj.phase.states:x0', phase.interp('x0', [4.5, 4.5]))
-        p.set_val('traj.phase.states:x1', phase.interp('x1', [0.0, 2 * np.pi / 3]))
-        p.set_val('traj.phase.states:x2', phase.interp('x2', [np.pi / 4, np.pi / 4]))
-        p.set_val('traj.phase.states:x3', phase.interp('x3', [0.0, 0.0]))
-        p.set_val('traj.phase.states:x4', phase.interp('x4', [0.0, 0.0]))
-        p.set_val('traj.phase.states:x5', phase.interp('x5', [0.0, 0.0]))
+
+        if isinstance(tx, dm.Birkhoff):
+            p.set_val('traj.phase.initial_states:x0', 4.5)
+            p.set_val('traj.phase.initial_states:x1', 0.0)
+            p.set_val('traj.phase.initial_states:x2', np.pi / 4)
+            p.set_val('traj.phase.initial_states:x3', 0.0)
+            p.set_val('traj.phase.initial_states:x4', 0.0)
+            p.set_val('traj.phase.initial_states:x5', 0.0)
+            p.set_val('traj.phase.final_states:x0', 4.5)
+            p.set_val('traj.phase.final_states:x1', 2 * np.pi / 3)
+            p.set_val('traj.phase.final_states:x2', np.pi / 4)
+            p.set_val('traj.phase.final_states:x3', 0.0)
+            p.set_val('traj.phase.final_states:x4', 0.0)
+            p.set_val('traj.phase.final_states:x5', 0.0)
+        else:
+            p.set_val('traj.phase.states:x0', phase.interp('x0', [4.5, 4.5]))
+            p.set_val('traj.phase.states:x1', phase.interp('x1', [0.0, 2 * np.pi / 3]))
+            p.set_val('traj.phase.states:x2', phase.interp('x2', [np.pi / 4, np.pi / 4]))
+            p.set_val('traj.phase.states:x3', phase.interp('x3', [0.0, 0.0]))
+            p.set_val('traj.phase.states:x4', phase.interp('x4', [0.0, 0.0]))
+            p.set_val('traj.phase.states:x5', phase.interp('x5', [0.0, 0.0]))
 
         return p
 
     def test_partials(self):
-        p = self.make_problem(transcription=Radau, optimizer='SLSQP')
+        p = self.make_problem(tx=dm.Radau(num_segments=30, order=3), optimizer='SLSQP')
         p.run_model()
         with printoptions(linewidth=1024, edgeitems=100):
             cpd = p.check_partials(method='fd', compact_print=True, out_stream=None)
 
     @require_pyoptsparse(optimizer='IPOPT')
     def test_robot_arm_radau(self):
-        p = self.make_problem(transcription=Radau, optimizer='IPOPT', numseg=12)
+        p = self.make_problem(tx=dm.Radau(num_segments=12, order=3), optimizer='IPOPT', numseg=12)
         dm.run_problem(p, refine_iteration_limit=5)
 
         t = p.get_val('traj.phase.timeseries.time')
@@ -94,9 +104,27 @@ class TestRobotArm(unittest.TestCase):
 
     @require_pyoptsparse(optimizer='IPOPT')
     def test_robot_arm_gl(self):
-        p = self.make_problem(transcription=GaussLobatto, optimizer='IPOPT', numseg=12)
+        p = self.make_problem(tx=dm.GaussLobatto(num_segments=20, order=3), optimizer='IPOPT')
         dm.run_problem(p, refine_iteration_limit=3)
 
+        t = p.get_val('traj.phase.timeseries.time')
+        assert_near_equal(t[-1], 9.14138, tolerance=1e-3)
+
+    @require_pyoptsparse(optimizer='IPOPT')
+    def test_robot_arm_birkhoff_lgl(self):
+        tx = dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=1, nodes_per_seg=30, grid_type='lgl'))
+        p = self.make_problem(tx=tx, optimizer='IPOPT')
+        dm.run_problem(p)
+
+        t = p.get_val('traj.phase.timeseries.time')
+        assert_near_equal(t[-1], 9.14138, tolerance=1e-3)
+
+    @require_pyoptsparse(optimizer='IPOPT')
+    def test_robot_arm_birkhoff_cgl(self):
+        tx = dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=1, nodes_per_seg=30, grid_type='cgl'))
+        p = self.make_problem(tx=tx, optimizer='IPOPT')
+        dm.run_problem(p)
+
         t = p.get_val('traj.phase.timeseries.time')
 
         assert_near_equal(t[-1], 9.14138, tolerance=1e-3)
diff --git a/dymos/examples/robot_arm/test/test_robot_arm.py b/dymos/examples/robot_arm/test/test_robot_arm.py
index 59a87621..5c595e2b 100644
--- a/dymos/examples/robot_arm/test/test_robot_arm.py
+++ b/dymos/examples/robot_arm/test/test_robot_arm.py
@@ -3,14 +3,13 @@ import unittest
 
 import numpy as np
 import matplotlib.pyplot as plt
+import openmdao.api as om
+import dymos as dm
 
-from openmdao.api import Problem, Group, pyOptSparseDriver
 from openmdao.utils.assert_utils import assert_near_equal
 from openmdao.utils.general_utils import printoptions
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 
-import dymos as dm
-from dymos import Trajectory, GaussLobatto, Phase, Radau
 from dymos.examples.robot_arm.robot_arm_ode import RobotArmODE
 
 plt.switch_backend('Agg')
@@ -27,9 +26,9 @@ class TestRobotArm(unittest.TestCase):
                 os.remove(filename)
 
     @require_pyoptsparse(optimizer='SLSQP')
-    def make_problem(self, transcription=Radau, optimizer='SLSQP', numseg=30):
-        p = Problem(model=Group())
-        p.driver = pyOptSparseDriver()
+    def make_problem(self, tx, optimizer='SLSQP', numseg=30):
+        p = om.Problem(model=om.Group())
+        p.driver = om.pyOptSparseDriver()
         p.driver.declare_coloring()
         p.driver.options['optimizer'] = optimizer
         if optimizer == 'SNOPT':
@@ -39,10 +38,9 @@ class TestRobotArm(unittest.TestCase):
             p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'
             p.driver.opt_settings['max_iter'] = 500
             p.driver.opt_settings['print_level'] = 0
-        traj = p.model.add_subsystem('traj', Trajectory())
+        traj = p.model.add_subsystem('traj', dm.Trajectory())
 
-        phase = traj.add_phase('phase', Phase(ode_class=RobotArmODE,
-                                              transcription=transcription(num_segments=numseg, order=3)))
+        phase = traj.add_phase('phase', dm.Phase(ode_class=RobotArmODE, transcription=tx))
         phase.set_time_options(fix_initial=True, fix_duration=False)
 
         phase.add_state('x0', fix_initial=True, fix_final=True, rate_source='x0_dot', units='m')
@@ -71,24 +69,39 @@ class TestRobotArm(unittest.TestCase):
 
         p.set_val('traj.phase.t_initial', 0)
         p.set_val('traj.phase.t_duration', 10)
-        p.set_val('traj.phase.states:x0', phase.interp('x0', [4.5, 4.5]))
-        p.set_val('traj.phase.states:x1', phase.interp('x1', [0.0, 2 * np.pi / 3]))
-        p.set_val('traj.phase.states:x2', phase.interp('x2', [np.pi / 4, np.pi / 4]))
-        p.set_val('traj.phase.states:x3', phase.interp('x3', [0.0, 0.0]))
-        p.set_val('traj.phase.states:x4', phase.interp('x4', [0.0, 0.0]))
-        p.set_val('traj.phase.states:x5', phase.interp('x5', [0.0, 0.0]))
+
+        if isinstance(tx, dm.Birkhoff):
+            p.set_val('traj.phase.initial_states:x0', 4.5)
+            p.set_val('traj.phase.initial_states:x1', 0.0)
+            p.set_val('traj.phase.initial_states:x2', np.pi / 4)
+            p.set_val('traj.phase.initial_states:x3', 0.0)
+            p.set_val('traj.phase.initial_states:x4', 0.0)
+            p.set_val('traj.phase.initial_states:x5', 0.0)
+            p.set_val('traj.phase.final_states:x0', 4.5)
+            p.set_val('traj.phase.final_states:x1', 2 * np.pi / 3)
+            p.set_val('traj.phase.final_states:x2', np.pi / 4)
+            p.set_val('traj.phase.final_states:x3', 0.0)
+            p.set_val('traj.phase.final_states:x4', 0.0)
+            p.set_val('traj.phase.final_states:x5', 0.0)
+        else:
+            p.set_val('traj.phase.states:x0', phase.interp('x0', [4.5, 4.5]))
+            p.set_val('traj.phase.states:x1', phase.interp('x1', [0.0, 2 * np.pi / 3]))
+            p.set_val('traj.phase.states:x2', phase.interp('x2', [np.pi / 4, np.pi / 4]))
+            p.set_val('traj.phase.states:x3', phase.interp('x3', [0.0, 0.0]))
+            p.set_val('traj.phase.states:x4', phase.interp('x4', [0.0, 0.0]))
+            p.set_val('traj.phase.states:x5', phase.interp('x5', [0.0, 0.0]))
 
         return p
 
     def test_partials(self):
-        p = self.make_problem(transcription=Radau, optimizer='SLSQP')
+        p = self.make_problem(tx=dm.Radau(num_segments=30, order=3), optimizer='SLSQP')
         p.run_model()
         with printoptions(linewidth=1024, edgeitems=100):
             cpd = p.check_partials(method='fd', compact_print=True, out_stream=None)
 
     @require_pyoptsparse(optimizer='IPOPT')
     def test_robot_arm_radau(self):
-        p = self.make_problem(transcription=Radau, optimizer='IPOPT', numseg=12)
+        p = self.make_problem(tx=dm.Radau(num_segments=12, order=3), optimizer='IPOPT', numseg=12)
         dm.run_problem(p)
 
         t = p.get_val('traj.phase.timeseries.time')
@@ -129,7 +142,59 @@ class TestRobotArm(unittest.TestCase):
 
     @require_pyoptsparse(optimizer='IPOPT')
     def test_robot_arm_gl(self):
-        p = self.make_problem(transcription=GaussLobatto, optimizer='IPOPT', numseg=20)
+        p = self.make_problem(tx=dm.GaussLobatto(num_segments=20, order=3), optimizer='IPOPT')
+        dm.run_problem(p)
+
+        t = p.get_val('traj.phase.timeseries.time')
+        rho = p.get_val('traj.phase.timeseries.x0')
+        theta = p.get_val('traj.phase.timeseries.x1')
+        phi = p.get_val('traj.phase.timeseries.x2')
+        u0 = p.get_val('traj.phase.timeseries.u0')
+        u1 = p.get_val('traj.phase.timeseries.u1')
+        u2 = p.get_val('traj.phase.timeseries.u2')
+
+        if show_plots:
+            fig, axs = plt.subplots(2, 3)
+            axs[0, 0].plot(t, rho)
+            axs[0, 1].plot(t, theta)
+            axs[0, 2].plot(t, phi)
+            axs[1, 0].plot(t, u0)
+            axs[1, 1].plot(t, u1)
+            axs[1, 2].plot(t, u2)
+            plt.show()
+
+        assert_near_equal(t[-1], 9.14138, tolerance=1e-3)
+
+    @require_pyoptsparse(optimizer='IPOPT')
+    def test_robot_arm_birkhoff_lgl(self):
+        tx = dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=1, nodes_per_seg=30, grid_type='lgl'))
+        p = self.make_problem(tx=tx, optimizer='IPOPT')
+        dm.run_problem(p)
+
+        t = p.get_val('traj.phase.timeseries.time')
+        rho = p.get_val('traj.phase.timeseries.x0')
+        theta = p.get_val('traj.phase.timeseries.x1')
+        phi = p.get_val('traj.phase.timeseries.x2')
+        u0 = p.get_val('traj.phase.timeseries.u0')
+        u1 = p.get_val('traj.phase.timeseries.u1')
+        u2 = p.get_val('traj.phase.timeseries.u2')
+
+        if show_plots:
+            fig, axs = plt.subplots(2, 3)
+            axs[0, 0].plot(t, rho)
+            axs[0, 1].plot(t, theta)
+            axs[0, 2].plot(t, phi)
+            axs[1, 0].plot(t, u0)
+            axs[1, 1].plot(t, u1)
+            axs[1, 2].plot(t, u2)
+            plt.show()
+
+        assert_near_equal(t[-1], 9.14138, tolerance=1e-3)
+
+    @require_pyoptsparse(optimizer='IPOPT')
+    def test_robot_arm_birkhoff_cgl(self):
+        tx = dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=1, nodes_per_seg=30, grid_type='cgl'))
+        p = self.make_problem(tx=tx, optimizer='IPOPT')
         dm.run_problem(p)
 
         t = p.get_val('traj.phase.timeseries.time')
diff --git a/dymos/examples/shuttle_reentry/test/test_reentry.py b/dymos/examples/shuttle_reentry/test/test_reentry.py
index 5c097e78..09974640 100644
--- a/dymos/examples/shuttle_reentry/test/test_reentry.py
+++ b/dymos/examples/shuttle_reentry/test/test_reentry.py
@@ -1,12 +1,12 @@
 import unittest
 
 import numpy as np
-
 import openmdao.api as om
+import dymos as dm
+
 from openmdao.utils.assert_utils import assert_near_equal, assert_check_partials
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 
-from dymos import Trajectory, GaussLobatto, Phase, Radau, run_problem
 from dymos.examples.shuttle_reentry.shuttle_ode import ShuttleODE
 
 # Expected results from Betts
@@ -17,8 +17,7 @@ expected_results = {'constrained': {'time': 2198.67, 'theta': 30.6255},
 @use_tempdirs
 class TestReentry(unittest.TestCase):
 
-    def make_problem(self, constrained=True, transcription=GaussLobatto, optimizer='SLSQP',
-                     numseg=50):
+    def make_problem(self, constrained=True, transcription=dm.GaussLobatto, optimizer='SLSQP',):
         p = om.Problem(model=om.Group())
         p.driver = om.pyOptSparseDriver(optimizer=optimizer)
         p.driver.declare_coloring()
@@ -28,13 +27,22 @@ class TestReentry(unittest.TestCase):
             p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'
             p.driver.opt_settings['print_level'] = 0
             p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'
-            p.driver.opt_settings['tol'] = 1.0E-7
+            p.driver.opt_settings['tol'] = 1.0E-2
+            p.driver.opt_settings['constr_viol_tol'] = 1.0E-6
             p.driver.opt_settings['mu_strategy'] = 'monotone'
-
-        traj = p.model.add_subsystem('traj', Trajectory())
+            p.driver.opt_settings['mu_init'] = 0.01
+            p.driver.opt_settings['bound_mult_init_method'] = 'mu-based'
+        elif optimizer == 'SNOPT':
+            p.driver.declare_coloring()
+            p.driver.opt_settings['iSumm'] = 6
+            p.driver.opt_settings['Major iterations limit'] = 100
+            p.driver.opt_settings['Major feasibility tolerance'] = 1.0E-5
+            p.driver.opt_settings['Major optimality tolerance'] = 1.0E-3
+
+        traj = p.model.add_subsystem('traj', dm.Trajectory())
         phase0 = traj.add_phase('phase0',
-                                Phase(ode_class=ShuttleODE,
-                                      transcription=transcription(num_segments=numseg, order=3)))
+                                dm.Phase(ode_class=ShuttleODE,
+                                         transcription=transcription))
 
         phase0.set_time_options(fix_initial=True, units='s', duration_ref=200)
         phase0.add_state('h', fix_initial=True, fix_final=True, units='ft', rate_source='hdot',
@@ -83,18 +91,32 @@ class TestReentry(unittest.TestCase):
         p.set_val('traj.phase0.controls:beta',
                   phase0.interp('beta', [-75 * np.pi / 180, 0 * np.pi / 180]), units='rad')
 
+        if isinstance(transcription, dm.Birkhoff):
+            p.set_val('traj.phase0.initial_states:h', 260000, units='ft')
+            p.set_val('traj.phase0.initial_states:gamma', -1 * np.pi / 180, units='rad')
+            p.set_val('traj.phase0.initial_states:phi', 0.0, units='rad')
+            p.set_val('traj.phase0.initial_states:psi', 90 * np.pi / 180, units='rad')
+            p.set_val('traj.phase0.initial_states:theta', 0.0, units='rad')
+            p.set_val('traj.phase0.initial_states:v', 25600, units='ft/s')
+
+            p.set_val('traj.phase0.final_states:h', 80000, units='ft')
+            p.set_val('traj.phase0.final_states:gamma', -5 * np.pi / 180, units='rad')
+            p.set_val('traj.phase0.final_states:v', 2500, units='ft/s')
+
         return p
 
     @require_pyoptsparse(optimizer='SLSQP')
     def test_partials(self):
-        p = self.make_problem(constrained=True, transcription=Radau, optimizer='SLSQP', numseg=5)
+        tx = dm.Radau(num_segments=5, order=3)
+        p = self.make_problem(constrained=True, transcription=tx, optimizer='SLSQP')
         p.run_model()
         cpd = p.check_partials(method='cs', compact_print=True, out_stream=None)
         assert_check_partials(cpd, atol=1.0E-4, rtol=1.1)
 
     @require_pyoptsparse(optimizer='IPOPT')
     def test_reentry_constrained_radau(self):
-        p = self.make_problem(constrained=True, transcription=Radau, optimizer='IPOPT')
+        tx = dm.Radau(num_segments=50, order=3)
+        p = self.make_problem(constrained=True, transcription=tx, optimizer='IPOPT')
 
         p.run_driver()
         assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1],
@@ -107,7 +129,36 @@ class TestReentry(unittest.TestCase):
 
     @require_pyoptsparse(optimizer='IPOPT')
     def test_reentry_constrained_gauss_lobatto(self):
-        p = self.make_problem(constrained=True, transcription=GaussLobatto, optimizer='IPOPT')
+        tx = dm.GaussLobatto(num_segments=50, order=3)
+        p = self.make_problem(constrained=True, transcription=tx, optimizer='IPOPT')
+
+        p.run_driver()
+        assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1],
+                          expected_results['constrained']['time'],
+                          tolerance=1e-2)
+
+        assert_near_equal(p.get_val('traj.phase0.timeseries.theta', units='deg')[-1],
+                          expected_results['constrained']['theta'],
+                          tolerance=1e-2)
+
+    @require_pyoptsparse(optimizer='IPOPT')
+    def test_reentry_constrained_birkhoff_lgl(self):
+        tx = dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=1, nodes_per_seg=30, grid_type='lgl'))
+        p = self.make_problem(constrained=True, transcription=tx, optimizer='IPOPT')
+
+        p.run_driver()
+        assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1],
+                          expected_results['constrained']['time'],
+                          tolerance=1e-2)
+
+        assert_near_equal(p.get_val('traj.phase0.timeseries.theta', units='deg')[-1],
+                          expected_results['constrained']['theta'],
+                          tolerance=1e-2)
+
+    @require_pyoptsparse(optimizer='IPOPT')
+    def test_reentry_constrained_birkhoff_cgl(self):
+        tx = dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=1, nodes_per_seg=30, grid_type='cgl'))
+        p = self.make_problem(constrained=True, transcription=tx, optimizer='IPOPT')
 
         p.run_driver()
         assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1],
@@ -120,7 +171,8 @@ class TestReentry(unittest.TestCase):
 
     @require_pyoptsparse(optimizer='IPOPT')
     def test_reentry_unconstrained_radau(self):
-        p = self.make_problem(constrained=False, transcription=Radau, optimizer='IPOPT')
+        tx = dm.Radau(num_segments=50, order=3)
+        p = self.make_problem(constrained=False, transcription=tx, optimizer='IPOPT')
 
         p.run_driver()
         assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1],
@@ -133,7 +185,36 @@ class TestReentry(unittest.TestCase):
 
     @require_pyoptsparse(optimizer='IPOPT')
     def test_reentry_unconstrained_gauss_lobatto(self):
-        p = self.make_problem(constrained=False, transcription=GaussLobatto, optimizer='IPOPT')
+        tx = dm.GaussLobatto(num_segments=50, order=3)
+        p = self.make_problem(constrained=False, transcription=tx, optimizer='IPOPT')
+
+        p.run_driver()
+        assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1],
+                          expected_results['unconstrained']['time'],
+                          tolerance=1e-2)
+
+        assert_near_equal(p.get_val('traj.phase0.timeseries.theta', units='deg')[-1],
+                          expected_results['unconstrained']['theta'],
+                          tolerance=1e-2)
+
+    @require_pyoptsparse(optimizer='IPOPT')
+    def test_reentry_unconstrained_birkhoff_lgl(self):
+        tx = dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=1, nodes_per_seg=60, grid_type='lgl'))
+        p = self.make_problem(constrained=False, transcription=tx, optimizer='IPOPT')
+
+        p.run_driver()
+        assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1],
+                          expected_results['unconstrained']['time'],
+                          tolerance=1e-2)
+
+        assert_near_equal(p.get_val('traj.phase0.timeseries.theta', units='deg')[-1],
+                          expected_results['unconstrained']['theta'],
+                          tolerance=1e-2)
+
+    @require_pyoptsparse(optimizer='IPOPT')
+    def test_reentry_unconstrained_birkhoff_cgl(self):
+        tx = dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=1, nodes_per_seg=60, grid_type='cgl'))
+        p = self.make_problem(constrained=False, transcription=tx, optimizer='IPOPT')
 
         p.run_driver()
         assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1],
@@ -156,10 +237,10 @@ class TestReentry(unittest.TestCase):
         p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'
         p.driver.opt_settings['mu_strategy'] = 'monotone'
 
-        traj = p.model.add_subsystem('traj', Trajectory())
+        traj = p.model.add_subsystem('traj', dm.Trajectory())
         phase0 = traj.add_phase('phase0',
-                                Phase(ode_class=ShuttleODE,
-                                      transcription=Radau(num_segments=30, order=3)))
+                                dm.Phase(ode_class=ShuttleODE,
+                                         transcription=dm.Radau(num_segments=30, order=3)))
 
         phase0.set_time_options(fix_initial=True, units='s', duration_ref=200)
         phase0.add_state('h', fix_initial=True, fix_final=True, units='ft', rate_source='hdot',
@@ -211,7 +292,7 @@ class TestReentry(unittest.TestCase):
         p.set_val('traj.phase0.polynomial_controls:beta',
                   phase0.interp('beta', [-20, 0]), units='deg')
 
-        run_problem(p, simulate=True)
+        dm.run_problem(p, simulate=True)
 
         sol = om.CaseReader('dymos_solution.db').get_case('final')
         sim = om.CaseReader('dymos_simulation.db').get_case('final')
diff --git a/dymos/transcriptions/pseudospectral/components/test/test_birkhoff_collocation_defect.py b/dymos/transcriptions/pseudospectral/components/test/test_birkhoff_collocation_defect.py
new file mode 100644
index 00000000..8fab8d4f
--- /dev/null
+++ b/dymos/transcriptions/pseudospectral/components/test/test_birkhoff_collocation_defect.py
@@ -0,0 +1,138 @@
+import unittest
+
+import numpy as np
+from numpy.testing import assert_almost_equal
+
+import openmdao.api as om
+from dymos.utils.testing_utils import assert_check_partials
+
+import dymos as dm
+from dymos.transcriptions.pseudospectral.components.birkhoff_collocation_comp import BirkhoffCollocationComp
+from dymos.transcriptions.grid_data import BirkhoffGrid
+
+# Modify class so we can run it standalone.
+from dymos.utils.misc import CompWrapperConfig
+
+
+class _PhaseStub():
+    """
+    A stand-in for the Phase during config_io for testing.
+    It just supports the classify_var method and returns "state", the only value needed for unittests.
+    """
+    def classify_var(self, name):
+        return None
+
+
+CollocationComp = CompWrapperConfig(BirkhoffCollocationComp, [_PhaseStub()])
+
+
+class TestCollocationComp(unittest.TestCase):
+
+    def make_problem(self, grid_type='lgl'):
+        dm.options['include_check_partials'] = True
+
+        gd = BirkhoffGrid(num_segments=1, segment_ends=np.array([0., 10.]),
+                          nodes_per_seg=21, grid_type=grid_type)
+        n = gd.subset_num_nodes['col']
+        tau = gd.node_stau
+        t = 5 * tau + 5
+
+        self.p = om.Problem(model=om.Group())
+
+        state_options = {'x': {'units': 'm', 'shape': (1,), 'fix_initial': True,
+                               'fix_final': False, 'solve_segments': False,
+                               'input_initial': False, 'rate_source': 'foo'},
+                         'y': {'units': 'm', 'shape': (2, 2), 'fix_initial': False,
+                               'fix_final': True, 'solve_segments': False,
+                               'input_initial': False, 'rate_source': 'bar'}}
+
+        indep_comp = om.IndepVarComp()
+        self.p.model.add_subsystem('indep', indep_comp, promotes_outputs=['*'])
+
+        # Testing the basic ODE xdot = -x, x(0) = 10
+        # Solution is x(t) = 10*exp(-t)
+
+        x_val = 10*np.exp(-t)
+        y_val = np.zeros((n, 2, 2))
+        y_val[:, 0, 0] = x_val
+        y_val[:, 1, 1] = x_val
+
+        indep_comp.add_output(
+            'dt_dstau',
+            val=np.ones(n) * 5, units='s')
+        indep_comp.add_output(
+            'state_value:x',
+            val=x_val, units='m')
+        indep_comp.add_output(
+            'f_value:x',
+            val=-x_val * 5, units='m')
+        indep_comp.add_output(
+            'f_computed:x',
+            val=-x_val, units='m/s')
+        indep_comp.add_output(
+            'state_value:y',
+            val=y_val, units='m')
+        indep_comp.add_output(
+            'f_value:y',
+            val=-y_val * 5, units='m')
+        indep_comp.add_output(
+            'f_computed:y',
+            val=-y_val, units='m/s')
+
+        self.p.model.add_subsystem('defect_comp',
+                                   subsys=CollocationComp(grid_data=gd,
+                                                          state_options=state_options,
+                                                          time_units='s'))
+
+        if grid_type == 'radau-ps':
+            src_indices = om.slicer[:-1]
+        else:
+            src_indices = om.slicer[:]
+
+        self.p.model.connect('state_value:x', 'defect_comp.states:x', src_indices=src_indices)
+        self.p.model.connect('f_value:x', 'defect_comp.state_rates:x', src_indices=src_indices)
+        self.p.model.connect('state_value:y', 'defect_comp.states:y')
+        self.p.model.connect('f_value:y', 'defect_comp.state_rates:y')
+        self.p.model.connect('f_computed:x', 'defect_comp.f_computed:x', src_indices=src_indices)
+        self.p.model.connect('f_computed:y', 'defect_comp.f_computed:y')
+        self.p.model.connect('dt_dstau', 'defect_comp.dt_dstau')
+
+        self.p.setup(force_alloc_complex=True)
+
+        self.p.set_val('defect_comp.initial_states:x', 10.0)
+        self.p.set_val('defect_comp.final_states:x', 10 * np.exp(-10))
+        self.p.set_val('defect_comp.state_segment_ends:x', [[10.0, 10 * np.exp(-10)]])
+
+        self.p.set_val('defect_comp.initial_states:y', np.array([[10.0, 0.0], [0.0, 10.0]]))
+        self.p.set_val('defect_comp.final_states:y', np.array([[10 * np.exp(-10), 0.0], [0.0, 10 * np.exp(-10)]]))
+        self.p.set_val('defect_comp.state_segment_ends:y', np.array([[[10.0, 0.0],
+                                                                      [0.0, 10.0]],
+                                                                     [[10 * np.exp(-10), 0.0],
+                                                                      [0.0, 10 * np.exp(-10)]]]))
+
+        self.p.run_model()
+
+    def tearDown(self):
+        dm.options['include_check_partials'] = False
+
+    def test_results(self):
+        for grid_type in ('lgl', 'cgl'):
+            with self.subTest(msg=grid_type):
+                self.make_problem(grid_type=grid_type)
+                assert_almost_equal(self.p['defect_comp.state_defects:x'], 0.0)
+                assert_almost_equal(self.p['defect_comp.state_rate_defects:x'], 0.0)
+                assert_almost_equal(self.p['defect_comp.final_state_defects:x'], 0.0)
+                assert_almost_equal(self.p['defect_comp.state_defects:y'], 0.0)
+                assert_almost_equal(self.p['defect_comp.state_rate_defects:y'], 0.0)
+                assert_almost_equal(self.p['defect_comp.final_state_defects:y'], 0.0)
+
+    def test_partials(self):
+        for grid_type in ('lgl', 'cgl'):
+            with self.subTest(msg=grid_type):
+                self.make_problem(grid_type=grid_type)
+                cpd = self.p.check_partials(compact_print=True, method='fd', out_stream=None)
+                assert_check_partials(cpd)
+
+
+if __name__ == '__main__':  # pragma: no cover
+    unittest.main()
diff --git a/dymos/transcriptions/pseudospectral/components/test/test_birkhoff_iter_group.py b/dymos/transcriptions/pseudospectral/components/test/test_birkhoff_iter_group.py
new file mode 100644
index 00000000..6a43da2f
--- /dev/null
+++ b/dymos/transcriptions/pseudospectral/components/test/test_birkhoff_iter_group.py
@@ -0,0 +1,183 @@
+import unittest
+
+import numpy as np
+
+import openmdao.api as om
+
+import dymos
+from openmdao.utils.assert_utils import assert_check_partials, assert_near_equal
+
+from dymos.utils.misc import GroupWrapperConfig
+from dymos.transcriptions.pseudospectral.components import BirkhoffIterGroup
+from dymos.phase.options import StateOptionsDictionary, TimeOptionsDictionary
+from dymos.transcriptions.grid_data import BirkhoffGrid
+
+
+class _PhaseStub():
+    """
+    A stand-in for the Phase during config_io for testing.
+    It just supports the classify_var method and returns "state", the only value needed for unittests.
+    """
+    def __init__(self):
+        self.nonlinear_solver = None
+        self.linear_solver = None
+
+    def classify_var(self, name):
+        return 'ode'
+
+
+BirkhoffIterGroup = GroupWrapperConfig(BirkhoffIterGroup, [_PhaseStub()])
+
+
+class SimpleODE(om.ExplicitComponent):
+    """
+    A simple ODE from https://math.okstate.edu/people/yqwang/teaching/math4513_fall11/Notes/rungekutta.pdf
+    """
+    def initialize(self):
+        self.options.declare('num_nodes', types=(int,))
+
+    def setup(self):
+        nn = self.options['num_nodes']
+        self.add_input('x', shape=(nn,), units='s**2')
+        self.add_input('t', shape=(nn,), units='s')
+        self.add_input('p', shape=(nn,), units='s**2')
+
+        self.add_output('x_dot', shape=(nn,), units='s')
+
+        ar = np.arange(nn, dtype=int)
+        self.declare_partials(of='x_dot', wrt='x', rows=ar, cols=ar, val=1.0)
+        self.declare_partials(of='x_dot', wrt='t', rows=ar, cols=ar)
+        self.declare_partials(of='x_dot', wrt='p', rows=ar, cols=ar, val=1.0)
+
+    def compute(self, inputs, outputs):
+        x = inputs['x']
+        t = inputs['t']
+        p = inputs['p']
+        outputs['x_dot'] = x - t**2 + p
+
+    def compute_partials(self, inputs, partials):
+        t = inputs['t']
+        partials['x_dot', 't'] = -2*t
+
+
+class TestBirkhoffIterGroup(unittest.TestCase):
+
+    def test_solve_segments_gl_fwd(self):
+
+        for grid_type in ['lgl', 'cgl']:
+            with self.subTest(msg=grid_type):
+                with dymos.options.temporary(include_check_partials=True):
+
+                    state_options = {'x': StateOptionsDictionary()}
+
+                    state_options['x']['shape'] = (1,)
+                    state_options['x']['units'] = 's**2'
+                    state_options['x']['targets'] = ['x']
+                    state_options['x']['initial_bounds'] = (None, None)
+                    state_options['x']['final_bounds'] = (None, None)
+                    state_options['x']['solve_segments'] = 'forward'
+                    state_options['x']['rate_source'] = 'x_dot'
+
+                    time_options = TimeOptionsDictionary()
+                    grid_data = BirkhoffGrid(num_segments=1, nodes_per_seg=31, grid_type=grid_type)
+                    nn = grid_data.subset_num_nodes['all']
+                    ode_class = SimpleODE
+
+                    p = om.Problem()
+                    p.model.add_subsystem('birkhoff', BirkhoffIterGroup(state_options=state_options,
+                                                                        time_options=time_options,
+                                                                        grid_data=grid_data,
+                                                                        ode_class=ode_class))
+
+                    birkhoff = p.model._get_subsystem('birkhoff')
+
+                    birkhoff.nonlinear_solver = om.NewtonSolver(solve_subsystems=True)
+                    birkhoff.linear_solver = om.DirectSolver()
+
+                    p.setup(force_alloc_complex=True)
+
+                    # Instead of using the TimeComp just transform the node segment taus onto [0, 2]
+                    times = grid_data.node_stau + 1
+
+                    solution = np.reshape(times**2 + 2 * times + 1 - 0.5 * np.exp(times), (nn, 1))
+                    dsolution_dt = np.reshape(2 * times + 2 - 0.5 * np.exp(times), (nn, 1))
+
+                    p.set_val('birkhoff.initial_states:x', 0.5)
+                    p.set_val('birkhoff.ode_all.t', times)
+                    p.set_val('birkhoff.ode_all.p', 1.0)
+
+                    # We don't need to provide guesses for these values in this case
+                    # p.set_val('birkhoff.final_states:x', solution[-1])
+                    # p.set_val('birkhoff.states:x', solution)
+                    # p.set_val('birkhoff.state_rates:x', dsolution_dt)
+
+                    p.final_setup()
+                    p.run_model()
+
+                    assert_near_equal(solution, p.get_val('birkhoff.states:x'), tolerance=1.0E-9)
+                    assert_near_equal(dsolution_dt, p.get_val('birkhoff.state_rates:x'), tolerance=1.0E-9)
+                    assert_near_equal(solution[np.newaxis, 0], p.get_val('birkhoff.initial_states:x'), tolerance=1.0E-9)
+                    assert_near_equal(solution[np.newaxis, -1], p.get_val('birkhoff.final_states:x'), tolerance=1.0E-9)
+
+                    cpd = p.check_partials(method='cs', compact_print=True, out_stream=None)
+                    assert_check_partials(cpd)
+
+    def test_solve_segments_gl_bkwd(self):
+
+        for grid_type in ['lgl', 'cgl']:
+            with self.subTest(msg=grid_type):
+                with dymos.options.temporary(include_check_partials=True):
+
+                    state_options = {'x': StateOptionsDictionary()}
+
+                    state_options['x']['shape'] = (1,)
+                    state_options['x']['units'] = 's**2'
+                    state_options['x']['targets'] = ['x']
+                    state_options['x']['initial_bounds'] = (None, None)
+                    state_options['x']['final_bounds'] = (None, None)
+                    state_options['x']['solve_segments'] = 'backward'
+                    state_options['x']['rate_source'] = 'x_dot'
+
+                    time_options = TimeOptionsDictionary()
+                    grid_data = BirkhoffGrid(num_segments=1, nodes_per_seg=31, grid_type=grid_type)
+                    nn = grid_data.subset_num_nodes['all']
+                    ode_class = SimpleODE
+
+                    p = om.Problem()
+                    p.model.add_subsystem('birkhoff', BirkhoffIterGroup(state_options=state_options,
+                                                                        time_options=time_options,
+                                                                        grid_data=grid_data,
+                                                                        ode_class=ode_class))
+
+                    birkhoff = p.model._get_subsystem('birkhoff')
+
+                    birkhoff.nonlinear_solver = om.NewtonSolver(solve_subsystems=True)
+                    birkhoff.linear_solver = om.DirectSolver()
+
+                    p.setup(force_alloc_complex=True)
+
+                    # Instead of using the TimeComp just transform the node segment taus onto [0, 2]
+                    times = grid_data.node_stau + 1
+
+                    solution = np.reshape(times**2 + 2 * times + 1 - 0.5 * np.exp(times), (nn, 1))
+                    dsolution_dt = np.reshape(2 * times + 2 - 0.5 * np.exp(times), (nn, 1))
+
+                    p.set_val('birkhoff.final_states:x', solution[-1])
+                    p.set_val('birkhoff.ode_all.t', times)
+                    p.set_val('birkhoff.ode_all.p', 1.0)
+
+                    # We don't need to provide guesses for these values in this case
+                    # p.set_val('birkhoff.final_states:x', solution[-1])
+                    # p.set_val('birkhoff.states:x', solution)
+                    # p.set_val('birkhoff.state_rates:x', dsolution_dt)
+
+                    p.final_setup()
+                    p.run_model()
+
+                    assert_near_equal(solution, p.get_val('birkhoff.states:x'), tolerance=1.0E-9)
+                    assert_near_equal(dsolution_dt, p.get_val('birkhoff.state_rates:x'), tolerance=1.0E-9)
+                    assert_near_equal(solution[np.newaxis, 0], p.get_val('birkhoff.initial_states:x'), tolerance=1.0E-9)
+                    assert_near_equal(solution[np.newaxis, -1], p.get_val('birkhoff.final_states:x'), tolerance=1.0E-9)
+
+                    cpd = p.check_partials(method='cs', compact_print=True, out_stream=None)
+                    assert_check_partials(cpd)
diff --git a/dymos/utils/test/test_birkhoff.py b/dymos/utils/test/test_birkhoff.py
new file mode 100644
index 00000000..022f34da
--- /dev/null
+++ b/dymos/utils/test/test_birkhoff.py
@@ -0,0 +1,26 @@
+import unittest
+
+import numpy as np
+from numpy.testing import assert_almost_equal
+
+from dymos.utils.birkhoff import birkhoff_matrix
+from dymos.utils.cgl import cgl, clenshaw_curtis
+from dymos.utils.lgl import lgl
+
+
+class TestBirkhoffMatrix(unittest.TestCase):
+
+    def test_birkhoff_matrix_lgl(self):
+        x, w = lgl(9)
+        B = birkhoff_matrix(x, w, grid_type='lgl')
+        assert_almost_equal(B[-1, :], w)
+
+    def test_birkhoff_matrix_cgl(self):
+        x, w = cgl(9)
+        B = birkhoff_matrix(x, w, grid_type='cgl')
+        _, wB = clenshaw_curtis(9)
+        assert_almost_equal(B[-1, :], wB)
+
+
+if __name__ == '__main__':  # pragma: no cover
+    unittest.main()
diff --git a/dymos/utils/test/test_cgl.py b/dymos/utils/test/test_cgl.py
new file mode 100644
index 00000000..647a68b7
--- /dev/null
+++ b/dymos/utils/test/test_cgl.py
@@ -0,0 +1,44 @@
+import unittest
+
+import numpy as np
+from numpy.testing import assert_almost_equal
+
+from dymos.utils.cgl import cgl
+
+# Known solutions
+x_i = {2: [-1.0, 1.0],
+       3: [-1.0, 0.0, 1.0],
+       4: [-1.0, -0.5, 0.5, 1.0],
+       5: [-1.0, -1/np.sqrt(2), 0.0, 1/np.sqrt(2), 1.0]}
+
+w_i = {2: [np.pi/2, np.pi/2],
+       3: [np.pi/4, np.pi/2, np.pi/4],
+       4: [np.pi/6, np.pi/3, np.pi/3, np.pi/6],
+       5: [np.pi/8, np.pi/4, np.pi/4, np.pi/4, np.pi/8]}
+
+
+class TestCGL(unittest.TestCase):
+
+    def test_nodes_and_weights_2(self):
+        x_2, w_2 = cgl(2)
+        assert_almost_equal(x_2, x_i[2], decimal=6)
+        assert_almost_equal(w_2, w_i[2], decimal=6)
+
+    def test_nodes_and_weights_3(self):
+        x_3, w_3 = cgl(3)
+        assert_almost_equal(x_3, x_i[3], decimal=6)
+        assert_almost_equal(w_3, w_i[3], decimal=6)
+
+    def test_nodes_and_weights_4(self):
+        x_4, w_4 = cgl(4)
+        assert_almost_equal(x_4, x_i[4], decimal=6)
+        assert_almost_equal(w_4, w_i[4], decimal=6)
+
+    def test_nodes_and_weights_5(self):
+        x_5, w_5 = cgl(5)
+        assert_almost_equal(x_5, x_i[5], decimal=6)
+        assert_almost_equal(w_5, w_i[5], decimal=6)
+
+
+if __name__ == '__main__':  # pragma: no cover
+    unittest.main()