Repository: OpenMDAO/dymos

Issue:
Add multiple-shooting capability
### Issue Type

- [ ] Bug
- [x] Enhancement
- [ ] Docs
- [ ] Miscellaneous

### Description

The experimental ExplicitShooting transcription currently performs single shooting.  That is, `phase.states:{name}` contains the initial value of the state at the start of the phase and it is propagated across each segment continuously.

dymos should also support multiple shooting, where `phase.states:{name}` contains the state value at the start of each segment, and each segment is propagated independently.  Constraints on discontinuities between state values at the segment boundaries should be provided to the optimizer.  The first state value in each segment should be a design variable (except for the first segment in the case of `fix_initial=True`).

### Example

None

### Environment

N/A

Hints from issue comments:

Solution patch:
diff --git a/.github/workflows/dymos_docs_workflow.yml b/.github/workflows/dymos_docs_workflow.yml
index 94dede7a4..d8b604654 100644
--- a/.github/workflows/dymos_docs_workflow.yml
+++ b/.github/workflows/dymos_docs_workflow.yml
@@ -25,16 +25,16 @@ jobs:
         include:
           # baseline versions except with pyoptsparse but no SNOPT
           # build docs to verify those that use pyoptsparse do not use SNOPT
-          - NAME: baseline_no_snopt
-            PY: '3.11'
-            NUMPY: '1.26'
-            SCIPY: '1.13'
-            PETSc: '3.19'
-            PYOPTSPARSE: 'v2.11.0'
-            OPENMDAO: 'latest'
-            OPTIONAL: '[docs]'
-            JAX: '0.4.28'
-            PUBLISH_DOCS: 1
+          # - NAME: baseline_no_snopt
+          #   PY: '3.11'
+          #   NUMPY: '1.26'
+          #   SCIPY: '1.13'
+          #   PETSc: '3.19'
+          #   PYOPTSPARSE: 'v2.11.0'
+          #   OPENMDAO: 'latest'
+          #   OPTIONAL: '[docs]'
+          #   JAX: '0.4.28'
+          #   PUBLISH_DOCS: 1
 
           # make sure the latest versions of things don't break the docs
           # sticking with Python 3.12 for now, 3.13 requires NumPy 2.1 which does not work yet with PETSc/pyoptsparse
@@ -49,7 +49,7 @@ jobs:
             OPENMDAO: 'dev'
             OPTIONAL: '[docs]'
             JAX: 'latest'
-            PUBLISH_DOCS: 0
+            PUBLISH_DOCS: 1
 
     steps:
       - name: Display run details
diff --git a/.github/workflows/publish-image.yml b/.github/workflows/publish-image.yml
index 9b025049d..9758dc9d0 100644
--- a/.github/workflows/publish-image.yml
+++ b/.github/workflows/publish-image.yml
@@ -1,4 +1,4 @@
-name: Publish Image
+name: Publish Images
 
 on:
   release:
@@ -7,16 +7,30 @@ on:
   # Run the workflow manually from the Actions tab
   workflow_dispatch:
 
-env:
-  # Will use docker.io for Docker Hub if empty
-  REGISTRY: ghcr.io
-  # github.repository as <account>/<repo>
-  IMAGE_NAME: ${{ github.repository }}
-
 jobs:
-  build:
 
-    runs-on: ubuntu-latest
+  build_image:
+
+    strategy:
+
+      fail-fast: false
+
+      matrix:
+        include:
+          - NAME: ubuntu
+            OS: ubuntu-24.04
+            FILE: Dockerfile
+            IMAGE_NAME: ${{ github.repository }}
+            REGISTRY: ghcr.io
+
+          - NAME: ubuntu-arm
+            OS: ubuntu-24.04-arm
+            FILE: Dockerfile.arm
+            IMAGE_NAME: ${{ github.repository }}-arm
+            REGISTRY: ghcr.io
+
+    runs-on: ${{ matrix.OS }}
+
     permissions:
       contents: read
       packages: write
@@ -37,11 +51,11 @@ jobs:
           git describe --tags > tag.txt
           echo "TAG=$(cat tag.txt)" >> $GITHUB_ENV
 
-      # Login against the registry
-      - name: Log into registry ${{ env.REGISTRY }}
+      # Login to container registry
+      - name: Log into registry ${{ matrix.REGISTRY }}
         uses: docker/login-action@v3
         with:
-          registry: ${{ env.REGISTRY }}
+          registry: ${{ matrix.REGISTRY }}
           username: ${{ github.actor }}
           password: ${{ secrets.GITHUB_TOKEN }}
 
@@ -50,7 +64,7 @@ jobs:
         id: meta
         uses: docker/metadata-action@v5
         with:
-          images: ${{ env.REGISTRY }}/${{ env.IMAGE_NAME }}
+          images: ${{ matrix.REGISTRY }}/${{ matrix.IMAGE_NAME }}
           tags: |
             type=raw,value=${{ env.TAG }}
 
@@ -59,6 +73,6 @@ jobs:
         uses: docker/build-push-action@v6
         with:
           push: true
-          context: ${{ inputs.image_name }}
+          file: ${{ matrix.FILE }}
           tags: ${{ steps.meta.outputs.tags }}
           labels: ${{ steps.meta.outputs.labels }}
diff --git a/Dockerfile b/Dockerfile
index 8e81ea9c7..e08d6fa7d 100644
--- a/Dockerfile
+++ b/Dockerfile
@@ -5,7 +5,7 @@ FROM ubuntu:24.04
 SHELL ["/bin/bash", "-c"]
 
 # Install updates
-RUN apt-get update -y && apt-get -y install wget git g++ gfortran make libblas-dev liblapack-dev
+RUN apt-get update -y && apt-get -y install sudo vim curl wget git g++ gfortran make cmake graphviz libblas-dev liblapack-dev
 
 RUN wget https://dl.google.com/linux/direct/google-chrome-stable_current_amd64.deb ;\
     apt-get install -y ./google-chrome-stable_current_amd64.deb
@@ -15,6 +15,7 @@ ENV USER=omdao
 RUN adduser --shell /bin/bash --disabled-password ${USER}
 RUN adduser ${USER} sudo
 RUN echo '%sudo ALL=(ALL) NOPASSWD:ALL' >> /etc/sudoers
+RUN chmod u+s /usr/bin/sudo
 
 USER ${USER}
 WORKDIR /home/${USER}
@@ -31,16 +32,16 @@ RUN source $HOME/miniforge3/etc/profile.d/conda.sh ;\
     #
     # Create conda environment
     #
-    conda create -n mdaowork python=3.12 'numpy<2' scipy cython swig -y ;\
+    conda create -n mdaowork python=3.12 'numpy<2' scipy cython swig -q -y ;\
     conda activate mdaowork ;\
-    conda install matplotlib graphviz -y ;\
-    conda install mpi4py petsc4py=3.20 -y ;\
+    conda install matplotlib graphviz -q -y ;\
+    conda install mpi4py openmpi petsc4py=3.20 pyoptsparse -q -y ;\
     python -m pip install pyparsing psutil objgraph plotly pyxdsm pydot ;\
     #
     # Install pyoptsparse
     #
     python -m pip install git+https://github.com/openmdao/build_pyoptsparse ;\
-    build_pyoptsparse -v ;\
+    # build_pyoptsparse -v ;\
     #
     # Install OpenMDAO
     #
diff --git a/Dockerfile.arm b/Dockerfile.arm
new file mode 100644
index 000000000..714e8e501
--- /dev/null
+++ b/Dockerfile.arm
@@ -0,0 +1,73 @@
+# Define image with OpenMDAO and Dymos
+
+
+FROM arm64v8/ubuntu:noble
+
+SHELL ["/bin/bash", "-c"]
+
+# Install updates
+RUN apt-get update -y && apt-get -y install sudo vim curl wget git g++ gfortran make cmake graphviz libblas-dev liblapack-dev
+
+# install a browser
+RUN apt-get update -y && apt-get install firefox --fix-missing -y
+
+# Create user
+ENV USER=omdao
+RUN adduser --shell /bin/bash --disabled-password ${USER}
+RUN adduser ${USER} sudo
+RUN echo '%sudo ALL=(ALL) NOPASSWD:ALL' >> /etc/sudoers
+RUN chmod u+s /usr/bin/sudo
+
+USER ${USER}
+WORKDIR /home/${USER}
+
+# Install Miniforge
+RUN wget -q -O Miniforge3.sh "https://github.com/conda-forge/miniforge/releases/latest/download/Miniforge3-$(uname)-$(uname -m).sh" ;\
+    bash Miniforge3.sh -b ;\
+    rm Miniforge3.sh ;\
+    export PATH=$HOME/miniforge3/bin:$PATH ;\
+    conda init bash
+
+# Create conda environment
+RUN source $HOME/miniforge3/etc/profile.d/conda.sh ;\
+    #
+    # Create conda environment
+    #
+    conda create -n mdaowork python=3.12 'numpy<2' scipy cython swig -q -y ;\
+    conda activate mdaowork ;\
+    conda install matplotlib graphviz -q -y ;\
+    conda install mpi4py openmpi petsc4py=3.20 pyoptsparse -q -y ;\
+    python -m pip install pyparsing psutil objgraph plotly pyxdsm pydot ;\
+    #
+    # Install build_pyoptsparse
+    # (this will allow the user additional options for installing pyoptsparse, beyond the conda install above)
+    #
+    python -m pip install git+https://github.com/openmdao/build_pyoptsparse ;\
+    # build_pyoptsparse -v ;\
+    #
+    # Install OpenMDAO
+    #
+    git clone https://github.com/OpenMDAO/OpenMDAO.git ;\
+    python -m pip install -e 'OpenMDAO[all]' ;\
+    #
+    # Install Dymos
+    #
+    git clone https://github.com/OpenMDAO/dymos.git ;\
+    python -m pip install -e 'dymos[all]'
+
+# Modify .bashrc
+RUN echo "## Always activate mdaowork environment on startup ##" >> ~/.bashrc ;\
+    echo "conda activate mdaowork" >> ~/.bashrc ;\
+    echo "" >> ~/.bashrc ;\
+    echo "## OpenMPI settings" >> ~/.bashrc ;\
+    echo "export OMPI_MCA_rmaps_base_oversubscribe=1" >> ~/.bashrc ;\
+    echo "export OMPI_MCA_btl=^openib" >> ~/.bashrc ;\
+    echo "" >> ~/.bashrc ;\
+    echo "## Required for some newer MPI / libfabric instances" >> ~/.bashrc ;\
+    echo "export RDMAV_FORK_SAFE=true" >> ~/.bashrc ;\
+    echo "" >> ~/.bashrc
+
+# Set up a work directory that can be shared with the host operating system
+WORKDIR /home/${USER}/work
+
+ENTRYPOINT ["tail", "-f", "/dev/null"]
diff --git a/docs/dymos_book/api/transcriptions_api.ipynb b/docs/dymos_book/api/transcriptions_api.ipynb
index a64c482fd..36ea707ce 100644
--- a/docs/dymos_book/api/transcriptions_api.ipynb
+++ b/docs/dymos_book/api/transcriptions_api.ipynb
@@ -140,12 +140,26 @@
     "om.show_options_table(tx)"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## PicardShooting Options"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": null,
-   "metadata": {},
+   "metadata": {
+    "tags": [
+     "remove-input"
+    ]
+   },
    "outputs": [],
-   "source": []
+   "source": [
+    "tx = dm.transcriptions.PicardShooting(num_segments=1, nodes_per_seg=5)\n",
+    "om.show_options_table(tx)"
+   ]
   }
  ],
  "metadata": {
@@ -159,7 +173,7 @@
    }
   },
   "kernelspec": {
-   "display_name": "Python 3 (ipykernel)",
+   "display_name": "py312",
    "language": "python",
    "name": "python3"
   },
@@ -173,7 +187,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.11.0"
+   "version": "3.12.8"
   }
  },
  "nbformat": 4,
diff --git a/docs/dymos_book/examples/implicit_duration/brachistochrone_implicit_duration.ipynb b/docs/dymos_book/examples/implicit_duration/brachistochrone_implicit_duration.ipynb
index 8fd7d6857..fa859e25d 100644
--- a/docs/dymos_book/examples/implicit_duration/brachistochrone_implicit_duration.ipynb
+++ b/docs/dymos_book/examples/implicit_duration/brachistochrone_implicit_duration.ipynb
@@ -47,23 +47,19 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# Implicit Duration\n",
+    "# Boundary Balance\n",
     "\n",
     "In some cases, using dymos without the need for an optimizer is desirable.\n",
     "Sometimes in the absense of a control, you may just want to propagate the problem.\n",
     "Solve segments allows us to do this for a specified amount of time, but if the duration of the phase is itself an implicit output, then we need dymos to understand the residual to associate with the phase duration.\n",
     "\n",
-    "The requirement for this procedure to work is that the dynamics have to be obeyed at each iteration.\n",
-    "A time-stepping transcription (`ExplicitShooting`) or a pseudospectral transcription with `solve_segments=True` is required.\n",
+    "This feature requires that some output of the system is directly coupled to one of the phase parameters, such as t_duration.\n",
+    "In a collocation phase where an optimizer satisfies the defect equations, changing the duration of the phase does not directly\n",
+    "impact the final conditions of the system.\n",
+    "The transcription needs to be such that invoking run_model is enough to generate outputs that change with the desired parameters.\n",
+    "The transcription must be one of `ExplicitShooting`, `PicardShooting`, or a pseudospectral transcription with `solve_segments=True`.\n",
     "\n",
-    "One important note: In the pseudospectral transcriptions, using `solve_segments` with option `compressed=False`, results in a \"multiple shooting\" approach, whereby the state discontinuitites between segments are treated as constraints to be nulled out by the optimizer. This approach should therefore only be used with \"single shooting\" approaches (`compressed=True` in the case of the pseudospectral transcriptions) \n",
-    "\n",
-    "## Solving the brachistochrone without an optimizer.\n",
-    "\n",
-    "Consider the brachistochrone problem. Below we implement the brachistochrone equations of motion using OpenMDAO's `JaxExplicitComponent` so that we can use automatic differentiation to provide the partials.\n",
-    "\n",
-    "The implicit solution to the brachistochrone is very robust, but what if we wanted to simplify it?\n",
-    "For instance, we leave the duration of the bead's trajectory as a design variable that we want to minimize, but we also require the enforcement of the final position of the brachistochrone to some (x, y) coordinate pair."
+    "## Solving the brachistochrone without an optimizer."
    ]
   },
   {
@@ -72,14 +68,13 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "import jax.numpy as jnp\n",
+    "import numpy as np\n",
     "import openmdao.api as om\n",
-    "import dymos as dm\n",
     "\n",
     "\n",
-    "class BrachistochroneODE(om.JaxExplicitComponent):\n",
+    "class BrachistochroneODE(om.ExplicitComponent):\n",
     "    \"\"\"\n",
-    "    The brachistochrone EOM assuming \n",
+    "    The brachistochrone EOM assuming\n",
     "    \"\"\"\n",
     "\n",
     "    def initialize(self):\n",
@@ -91,8 +86,7 @@
     "        # static parameters\n",
     "        self.add_input('v', desc='speed of the bead on the wire', shape=(nn,), units='m/s')\n",
     "        self.add_input('theta', desc='angle between wire tangent and the nadir', shape=(nn,), units='rad')\n",
-    "        self.add_input('g', desc='gravitational acceleration', tags=['dymos.static_target'], units='m/s**2',\n",
-    "                       shape=(nn,))\n",
+    "        self.add_input('g', desc='gravitational acceleration', shape=(1,), tags=['dymos.static_target'], units='m/s**2')\n",
     "\n",
     "        self.add_output('xdot', desc='velocity component in x', shape=(nn,), units='m/s',\n",
     "                        tags=['dymos.state_rate_source:x', 'dymos.state_units:m'])\n",
@@ -103,28 +97,17 @@
     "        self.add_output('vdot', desc='acceleration magnitude', shape=(nn,), units='m/s**2',\n",
     "                        tags=['dymos.state_rate_source:v', 'dymos.state_units:m/s'])\n",
     "\n",
-    "        # Setup partials\n",
-    "        ar = jnp.arange(nn, dtype=int)\n",
-    "\n",
-    "        self.declare_partials(of='vdot', wrt='v', rows=ar, cols=ar)\n",
-    "        self.declare_partials(of='vdot', wrt='theta', rows=ar, cols=ar)\n",
-    "        self.declare_partials(of='vdot', wrt='g')\n",
+    "        self.declare_coloring(method='cs')\n",
     "\n",
-    "        self.declare_partials(of='xdot', wrt='v', rows=ar, cols=ar)\n",
-    "        self.declare_partials(of='xdot', wrt='theta', rows=ar, cols=ar)\n",
+    "    def compute(self, inputs, outputs):\n",
+    "        v, theta, g = inputs.values()\n",
     "\n",
-    "        self.declare_partials(of='ydot', wrt='v', rows=ar, cols=ar)\n",
-    "        self.declare_partials(of='ydot', wrt='theta', rows=ar, cols=ar)\n",
+    "        sin_theta = np.sin(theta)\n",
+    "        cos_theta = np.cos(theta)\n",
     "\n",
-    "    def compute_primal(self, v, theta, g):\n",
-    "        sin_theta = jnp.sin(theta)\n",
-    "        cos_theta = jnp.cos(theta)\n",
-    "\n",
-    "        vdot = g * cos_theta\n",
-    "        xdot = v * sin_theta\n",
-    "        ydot = -v * cos_theta\n",
-    "\n",
-    "        return xdot, ydot, vdot"
+    "        outputs['vdot'] = g * cos_theta\n",
+    "        outputs['xdot'] = v * sin_theta\n",
+    "        outputs['ydot'] = -v * cos_theta"
    ]
   },
   {
@@ -135,7 +118,13 @@
     "\n",
     "In the following example, we set `solve_segments='forward'` such that all of the collocation defects associate with the Radau transcription are satisfied by a Newton solver rather than the optimizer.\n",
     "This is the first step to solving the problem without an optimizer.\n",
-    "This ensures that the equations of motion are satisfied to the extent possible assuming the transcription grid."
+    "This ensures that the equations of motion are satisfied to the extent possible assuming the transcription grid.\n",
+    "\n",
+    "In pseudospectral phases, solve_segments is an optional setting that can be applied to all states in the phase via the `solve_segments` option on the pseudospectral transcriptions, or to individual states.\n",
+    "\n",
+    "In phases that use PicardShooting, `solve_segments` must be `'forward'` or `'backward'`. The default direction of the solve is `'forward'`, but individual states may override this.\n",
+    "\n",
+    "You can experiment with using different transcriptions below by choosing which value of `TX` you want to use. Note that because the initial states of the brachistochrone are fixed, `solve_segments='forward'` should be used."
    ]
   },
   {
@@ -144,6 +133,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "# TX = dm.PicardShooting(num_segments=5, nodes_per_seg=4, solve_segments='forward')\n",
+    "TX = dm.Radau(num_segments=10, order=3, solve_segments='forward', compressed=True)\n",
+    "# TX = dm.ExplicitShooting(num_segments=10, order=3)\n",
+    "\n",
     "# Define the OpenMDAO problem\n",
     "p = om.Problem(model=om.Group())\n",
     "\n",
@@ -151,9 +144,9 @@
     "traj = dm.Trajectory()\n",
     "p.model.add_subsystem('traj', subsys=traj)\n",
     "\n",
-    "# Define a Dymos Phase object with Radau Transcription\n",
+    "# Define a Dymos Phase\n",
     "phase = dm.Phase(ode_class=BrachistochroneODE,\n",
-    "                 transcription=dm.Radau(num_segments=5, order=3, solve_segments='forward'))\n",
+    "                 transcription=TX)\n",
     "traj.add_phase(name='phase0', phase=phase)\n",
     "\n",
     "# Set the time options\n",
@@ -181,7 +174,7 @@
     "# Doing so can significant speed up the execution of dymos.\n",
     "p.driver.declare_coloring()\n",
     "\n",
-    "phase.nonlinear_solver = om.NewtonSolver(maxiter=1000, solve_subsystems=True, stall_limit=3, iprint=2)\n",
+    "phase.nonlinear_solver = om.NewtonSolver(maxiter=100, solve_subsystems=True, stall_limit=3, iprint=2)\n",
     "phase.linear_solver = om.DirectSolver()\n",
     "\n",
     "# Setup the problem\n",
@@ -229,9 +222,9 @@
    "source": [
     "## Stopping the propagation at the desired time.\n",
     "\n",
-    "We can utilize the `set_implicit_duration` method on `Phase` to turn `t_duration` into an implicit output and provide it with a residual.\n",
+    "We can utilize the `add_boundary_balance` method on `Phase` to turn `t_duration` into an implicit output and provide it with a residual.\n",
     "\n",
-    "In our case, we set `phase.set_duration_balance('x', 10.0)`\n",
+    "In our case, we set `phase.add_boundary_balance(param='t_duration', name='x', tgt_val=10.0)`\n",
     "\n",
     "This specifies that we want `t_duration` to be associated with the following residual:\n",
     "\n",
@@ -248,6 +241,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "#TX = dm.PicardShooting(num_segments=5, nodes_per_seg=4, solve_segments='forward')\n",
+    "TX = dm.Radau(num_segments=10, order=3, solve_segments='forward', compressed=True)\n",
+    "# TX = dm.ExplicitShooting(num_segments=10, order=3)\n",
+    "\n",
     "# Define the OpenMDAO problem\n",
     "p = om.Problem(model=om.Group())\n",
     "\n",
@@ -255,15 +252,14 @@
     "traj = dm.Trajectory()\n",
     "p.model.add_subsystem('traj', subsys=traj)\n",
     "\n",
-    "# Define a Dymos Phase object with Radau Transcription\n",
+    "# Define a Dymos Phase\n",
     "phase = dm.Phase(ode_class=BrachistochroneODE,\n",
-    "                 transcription=dm.Radau(num_segments=5, order=3, solve_segments='forward', compressed=True))\n",
+    "                 transcription=TX)\n",
     "traj.add_phase(name='phase0', phase=phase)\n",
     "\n",
     "# Set the time options\n",
-    "phase.set_time_options(fix_initial=True,\n",
-    "                       duration_bounds=(0.5, 10.0))\n",
-    "phase.set_duration_balance('x', 10.0)\n",
+    "phase.set_time_options(fix_initial=True)\n",
+    "phase.add_boundary_balance(param='t_duration', name='x', tgt_val=10.0)\n",
     "\n",
     "# Set the state options\n",
     "phase.add_state('x', rate_source='xdot',\n",
@@ -325,7 +321,7 @@
     "\n",
     "print('final time', p.get_val('traj.phase0.timeseries.time')[-1])\n",
     "print('final y', y[-1])\n",
-    "print('final theta', jnp.degrees(theta[-1]))"
+    "print('final theta', np.degrees(theta[-1]))"
    ]
   },
   {
@@ -355,6 +351,10 @@
    "metadata": {},
    "outputs": [],
    "source": [
+    "# TX = dm.PicardShooting(num_segments=5, nodes_per_seg=4, solve_segments='forward')\n",
+    "TX = dm.Radau(num_segments=10, order=3, solve_segments='forward', compressed=True)\n",
+    "# TX = dm.ExplicitShooting(num_segments=10, order=3)\n",
+    "\n",
     "# Define the OpenMDAO problem\n",
     "p = om.Problem(model=om.Group())\n",
     "\n",
@@ -362,14 +362,13 @@
     "traj = dm.Trajectory()\n",
     "p.model.add_subsystem('traj', subsys=traj)\n",
     "\n",
-    "# Define a Dymos Phase object with Radau Transcription\n",
+    "# Define a Dymos Phase\n",
     "phase = dm.Phase(ode_class=BrachistochroneODE,\n",
-    "                 transcription=dm.Radau(num_segments=10, order=3, solve_segments='forward', compressed=True))\n",
+    "                 transcription=TX)\n",
     "traj.add_phase(name='phase0', phase=phase)\n",
     "\n",
     "# Set the time options\n",
     "phase.set_time_options(fix_initial=True)\n",
-    "phase.set_duration_balance('x', 10.0)\n",
     "\n",
     "# Set the state options\n",
     "phase.add_state('x', rate_source='xdot',\n",
@@ -399,11 +398,9 @@
     "phase.nonlinear_solver.linesearch = None\n",
     "phase.linear_solver = om.DirectSolver()\n",
     "\n",
-    "theta_rate_balance = p.model.add_subsystem('theta_rate_balance', om.BalanceComp())\n",
-    "\n",
-    "theta_rate_balance.add_balance('theta_rate', units='rad/s', eq_units='m', lhs_name='y_final', rhs_val=5.0)\n",
-    "p.model.connect('theta_rate_balance.theta_rate', 'traj.phase0.parameters:theta_rate')\n",
-    "p.model.connect('traj.phase0.timeseries.y', 'theta_rate_balance.y_final', src_indices=[-1])\n",
+    "# Define balance residuals to be converged by varying t_duration and theta_rate.\n",
+    "phase.add_boundary_balance('t_duration', name='x', tgt_val=10.0, loc='final', eq_units='m')\n",
+    "phase.add_boundary_balance('theta_rate', name='y', tgt_val=5.0, loc='final', eq_units='m')\n",
     "\n",
     "# Now we need a solver to converge the loop around the entire problem.\n",
     "p.model.nonlinear_solver = om.NewtonSolver(maxiter=100, solve_subsystems=True, stall_limit=3,\n",
@@ -422,7 +419,7 @@
     "phase.set_parameter_val('g', val=9.80665)\n",
     "\n",
     "# Note that we set theta_rate at the balance comp, which is then passed into the phase as a parameter.\n",
-    "p.model.set_val('theta_rate_balance.theta_rate', val=0.5, units='rad/s')\n",
+    "phase.set_parameter_val('theta_rate', 0.5, units='rad/s')\n",
     "\n",
     "# Run the driver to solve the problem\n",
     "p.run_model()\n"
@@ -449,20 +446,20 @@
     "\n",
     "print('final time', p.get_val('traj.phase0.timeseries.time')[-1])\n",
     "print('final y', y[-1])\n",
-    "print('final theta', jnp.degrees(theta[-1]))"
+    "print('final theta', np.degrees(theta[-1]))"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Using a combination of solve_segments, implicit phase duration, and an outer loop to make a phase parameter an implicit output, we've managed to solve a two point boundary value problem without the need of an optimizer. Since we're solving the segments in an inner loop here, each converged phase provides a physically valid traejctory, though it may not satisfy our desired final conditions unless the solvers are converged at all levels.\n",
+    "Using a combination of solve_segments and boundary balances, we've managed to solve the brachistochrone problem without the need of an optimizer.\n",
     "\n",
     "Note that solving problems in this way can be filled with pitfalls.\n",
     "Here we have a formulation of the brachistochrone that is robust to a variety of values for the parameter `theta_rate`.\n",
     "An alternative formulation is to keep the ratio $\\frac{\\sin{\\theta}}{v}$ constant. However, using the inverse sine in that formulation can easily result in domain issues which makes convergence difficult.\n",
     "\n",
-    "Even in this case, setting the `solve_subsystems` option on the outermost Newton solver to `False` seems to be too challenging and the solver fails to converge."
+    "Unlike `solve_segments` options in pseudospectral phases, the PicardShooting method uses multiple levels of solvers. The inner-most solver converges the state history in each segment. Outside of that a \"multiple shooting\" solver ensures state continuity between segments within a phase. Finally, the entire phase uses a NewtonSolver to find the value of any \"boundary balance\" outputs by driving their associated residuals to zero. In some situations"
    ]
   }
  ],
@@ -477,7 +474,7 @@
    }
   },
   "kernelspec": {
-   "display_name": "py312",
+   "display_name": "py311_2",
    "language": "python",
    "name": "python3"
   },
@@ -491,7 +488,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.12.8"
+   "version": "3.11.4"
   }
  },
  "nbformat": 4,
diff --git a/dymos/__init__.py b/dymos/__init__.py
index cbb6bb70c..63aba722a 100644
--- a/dymos/__init__.py
+++ b/dymos/__init__.py
@@ -2,8 +2,10 @@
 
 
 from .phase import Phase, AnalyticPhase
-from .transcriptions import GaussLobatto, Radau, ExplicitShooting, Analytic, Birkhoff
-from .transcriptions.grid_data import GaussLobattoGrid, RadauGrid, UniformGrid, BirkhoffGrid
+from .transcriptions import GaussLobatto, Radau, ExplicitShooting, Analytic, \
+    Birkhoff, PicardShooting
+from .transcriptions.grid_data import GaussLobattoGrid, ChebyshevGaussLobattoGrid, \
+    RadauGrid, UniformGrid, BirkhoffGrid
 from .trajectory.trajectory import Trajectory
 from .run_problem import run_problem
 from .load_case import load_case
diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index d7e28c50a..a582c2c4b 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -21,7 +21,7 @@
     TimeseriesOutputOptionsDictionary, PhaseTimeseriesOptionsDictionary
 
 from ..transcriptions.transcription_base import TranscriptionBase
-from ..transcriptions.grid_data import GaussLobattoGrid, RadauGrid, UniformGrid, BirkhoffGrid
+from ..transcriptions.grid_data import ChebyshevGaussLobattoGrid, GaussLobattoGrid, RadauGrid, UniformGrid, BirkhoffGrid
 from ..transcriptions import ExplicitShooting, GaussLobatto, Radau
 from ..utils.indexing import get_constraint_flat_idxs
 from ..utils.introspection import configure_time_introspection, _configure_constraint_introspection, \
@@ -197,6 +197,8 @@ def initialize(self):
                              desc='Options for each parameter in this phase.')
         self.options.declare('control_options', types=dict, default={},
                              desc='Options for each control in this phase.')
+        self.options.declare('boundary_balance_options', types=dict, default={},
+                             desc='Options specifying boundary conditions to be satisfied with a solver.')
 
     @property
     def time_options(self):
@@ -227,6 +229,10 @@ def parameter_options(self):
     def control_options(self):
         return self.options['control_options']
 
+    @property
+    def boundary_balance_options(self):
+        return self.options['boundary_balance_options']
+
     def add_state(self, name, units=_unspecified, shape=_unspecified,
                   rate_source=_unspecified, targets=_unspecified,
                   val=_unspecified, fix_initial=_unspecified, fix_final=_unspecified,
@@ -2087,6 +2093,49 @@ def set_parameter_val(self, name, val=None, units=None):
             units = self.parameter_options[name]['units']
         self.set_val(f'parameters:{name}', val=val, units=units)
 
+    def add_boundary_balance(self, param, name, tgt_val=0.0, loc='final', index=0, **kwargs):
+        """
+        Turn param into an implicit output in the phase when using shooting methods.
+
+        This capability requires the variable specified by `name` to be dependent upon the given `param` name, and
+        thus won't work with purely implicit transcriptions.
+        The transcription must be one of ExplicitShooting, PicardShooting, Analytic, or Pseudospectral with `solve_segments=True`.
+
+        Param will be output by a BalanceComp at the end of the phase. with a residual defined by `name - val` at either
+        the initial or final point in the phase (specified by 'loc').
+
+        Param is the name of the implicit variable within the phase. Typically one of 't_initial', 't_duration',
+        'initial_states:{state_name}', 'final_states:{state_name}', or 'parameters:{param_name}'.
+
+        Parameters
+        ----------
+        param : str
+            The dymos phase variable whose value is being set to satisfy the residual.
+        name : str
+            The state or time-dependent phase output that provides the residual value for param.
+        tgt_val : float or array-like, optional
+            The target value for `name`. Default is 0.0.
+        loc : str, optional
+            Whether the value given by name is being evaluated at the 'initial' or 'final'
+            point in the phase. Default is 'final'.
+        index : int or slice, optional
+            If the variable given by name is an array at each point in time, C-order index into
+            that variable at the intiial or final time. Default is 0.
+        **kwargs : dict
+            Additional keyword arguments forwarded to the `add_balance` call of an OpenMDAO
+            BalanceComp used to define the residual.
+        """
+        bbos = self.boundary_balance_options
+        if param in bbos:
+            name = bbos[param]['name']
+            tgt_val = bbos[param]['tgt_val']
+            loc = bbos[param]['loc']
+            raise ValueError(f'Phase variable {param} is already an implicit output for a boundary balance.\n'
+                             f'R({param}) = {name}[{loc}] - {tgt_val} = 0')
+
+        bbos[param] = {'param': param, 'name': name, 'tgt_val': tgt_val, 'loc': loc, 'index': index}
+        bbos[param].update(kwargs)
+
     def set_duration_balance(self, name, val=0.0, index=None, units=None, mult_val=None, normalize=False):
         """
         Adds a condition for the duration of the phase. This is satisfied using a nonlinear solver.
@@ -2103,43 +2152,19 @@ def set_duration_balance(self, name, val=0.0, index=None, units=None, mult_val=N
             If variable is an array at each point in time, this indicates which index is to be
             used as the objective, assuming C-ordered flattening.
         units : str, optional
-            The units of the objective function.  If None, use the units associated with the target.
-            If provided, must be compatible with the target units.
+            The units of the duration residual.
         mult_val : float, optional
             Default value for the LHS multiplier.
         normalize : bool, optional
             Specifies whether the resulting residual should be normalized by a quadratic
             function of the RHS.
         """
-        if self.time_options['fix_duration']:
-            raise ValueError('Cannot implicitly solve for phase duration when fix_duration is True')
-        elif self.time_options['input_duration']:
-            raise ValueError('Cannot implicitly solve for phase duration when input_duration is True')
-
-        if isinstance(self.options['transcription'], ExplicitShooting):
-            raise NotImplementedError('Transcription ExplicitShooting does not implement method setup_duration_balance')
-
-        options = {'name': name,
-                   'val': val,
-                   'index': index,
-                   'units': units,
-                   'mult_val': mult_val,
-                   'normalize': normalize}
-
-        if '=' in name:
-            balance_name = name.split('=')[0].strip()
-            self.add_calc_expr(name)
-        else:
-            balance_name = name
-
-        options['balance_name'] = balance_name
-
-        self.time_options['t_duration_balance_options'] = options
-
-        var_type = self.classify_var(name)
-        if var_type == 'ode':
-            if balance_name not in self._timeseries['timeseries']['outputs']:
-                self.add_timeseries_output(name, output_name=balance_name, units=units)
+        issue_warning("phase.set_duration_balance is deprecated. "
+                      "In the futurec use phase.add_boundary_balance('t_duration', ...",
+                      category=om.OMDeprecationWarning)
+        idx = 0 if index is None else index
+        self.add_boundary_balance('t_duration', name=name, tgt_val=val, loc='final', index=idx,
+                                  eq_units=units, mult_val=mult_val, normalize=normalize)
 
     def classify_var(self, var):
         """
@@ -2215,7 +2240,7 @@ def setup(self):
 
         transcription.setup_timeseries_outputs(self)
 
-        transcription.setup_duration_balance(self)
+        transcription.setup_boundary_balance(self)
 
         transcription.setup_defects(self)
         transcription.setup_solvers(self)
@@ -2256,6 +2281,7 @@ def configure(self):
         _configure_constraint_introspection(self)
 
         transcription.configure_boundary_constraints(self)
+        transcription.configure_boundary_balance(self)
 
         transcription.configure_path_constraints(self)
 
@@ -2269,8 +2295,6 @@ def configure(self):
 
         transcription.configure_timeseries_outputs(self)
 
-        transcription.configure_duration_balance(self)
-
         transcription.configure_solvers(self)
 
     def check_time_options(self):
@@ -2541,8 +2565,8 @@ def get_simulation_phase(self, times_per_seg=_unspecified, method=_unspecified,
         elif isinstance(self_tx, Radau):
             grid = RadauGrid(num_segments=num_seg, nodes_per_seg=seg_order + 1, segment_ends=seg_ends,
                              compressed=compressed)
-        elif isinstance(self_tx.grid_data, GaussLobattoGrid) or \
-                isinstance(self_tx.grid_data, RadauGrid) or isinstance(self_tx.grid_data, BirkhoffGrid):
+        elif isinstance(self_tx.grid_data, (GaussLobattoGrid, RadauGrid,
+                                            BirkhoffGrid, ChebyshevGaussLobattoGrid)):
             grid = self_tx.grid_data
         else:
             raise RuntimeError(f'Unexpected grid class for {self_tx.grid_data}. Only phases with GaussLobatto '
diff --git a/dymos/transcriptions/__init__.py b/dymos/transcriptions/__init__.py
index d4e2a1fc1..a971e497e 100644
--- a/dymos/transcriptions/__init__.py
+++ b/dymos/transcriptions/__init__.py
@@ -3,3 +3,4 @@
 from .pseudospectral.gauss_lobatto import GaussLobatto
 from .pseudospectral.radau_pseudospectral import Radau
 from .pseudospectral.birkhoff import Birkhoff
+from .picard_shooting.picard_shooting import PicardShooting
diff --git a/dymos/transcriptions/analytic/analytic.py b/dymos/transcriptions/analytic/analytic.py
index 0ee2b4339..70fd39de5 100644
--- a/dymos/transcriptions/analytic/analytic.py
+++ b/dymos/transcriptions/analytic/analytic.py
@@ -239,28 +239,6 @@ def setup_defects(self, phase):
         """
         pass
 
-    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 for this transcription.
diff --git a/dymos/transcriptions/common/control_comp.py b/dymos/transcriptions/common/control_comp.py
index 4f6065a26..74a1449b4 100644
--- a/dymos/transcriptions/common/control_comp.py
+++ b/dymos/transcriptions/common/control_comp.py
@@ -303,6 +303,73 @@ def _configure_desvars(self):
 
                 # self.set_input_defaults(name=dvname, val=default_val, units=options['units'])
 
+    def _configure_desvars(self):
+        control_options = self.options['control_options']
+        gd = self.options['grid_data']
+
+        for name, options in control_options.items():
+            if options['control_type'] == 'polynomial':
+                num_input_nodes = options['order'] + 1
+                shape = options['shape']
+                # default_val = reshape_val(options['val'], shape, num_input_nodes)
+                if options['opt']:
+
+                    desvar_indices = np.arange(num_input_nodes, dtype=int)
+                    if options['fix_initial']:
+                        desvar_indices = desvar_indices[1:]
+                    if options['fix_final']:
+                        desvar_indices = desvar_indices[:-1]
+
+                    lb = -INF_BOUND if options['lower'] is None else options['lower']
+                    ub = INF_BOUND if options['upper'] is None else options['upper']
+
+                    self.add_design_var(f'controls:{name}',
+                                        lower=lb,
+                                        upper=ub,
+                                        ref=options['ref'],
+                                        ref0=options['ref0'],
+                                        adder=options['adder'],
+                                        scaler=options['scaler'],
+                                        indices=desvar_indices,
+                                        flat_indices=True)
+
+                # self.set_input_defaults(name=f'controls:{name}', val=default_val, units=options['units'])
+            else:
+                num_input_nodes = gd.subset_num_nodes['control_input']
+
+                dvname = f'controls:{name}'
+                shape = options['shape']
+                size = np.prod(shape)
+                if options['opt']:
+                    desvar_indices = get_desvar_indices(size, num_input_nodes,
+                                                        options['fix_initial'], options['fix_final'])
+
+                    if len(desvar_indices) > 0:
+                        coerce_desvar_option = CoerceDesvar(num_input_nodes, desvar_indices,
+                                                            options=options)
+
+                        lb = np.zeros_like(desvar_indices, dtype=float)
+                        lb[:] = -INF_BOUND if coerce_desvar_option('lower') is None else \
+                            coerce_desvar_option('lower')
+
+                        ub = np.zeros_like(desvar_indices, dtype=float)
+                        ub[:] = INF_BOUND if coerce_desvar_option('upper') is None else \
+                            coerce_desvar_option('upper')
+
+                        self.add_design_var(name=dvname,
+                                            lower=lb,
+                                            upper=ub,
+                                            scaler=coerce_desvar_option('scaler'),
+                                            adder=coerce_desvar_option('adder'),
+                                            ref0=coerce_desvar_option('ref0'),
+                                            ref=coerce_desvar_option('ref'),
+                                            indices=desvar_indices,
+                                            flat_indices=True)
+
+                # default_val = reshape_val(options['val'], shape, num_input_nodes)
+
+                # self.set_input_defaults(name=dvname, val=default_val, units=options['units'])
+
     def configure_io(self):
         """
         I/O creation is delayed until configure so we can determine shape and units for the states.
diff --git a/dymos/transcriptions/explicit_shooting/explicit_shooting.py b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
index c6de23f01..d31b2c86f 100644
--- a/dymos/transcriptions/explicit_shooting/explicit_shooting.py
+++ b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
@@ -7,7 +7,7 @@
 
 from .explicit_shooting_continuity_comp import ExplicitShootingContinuityComp
 from ..transcription_base import TranscriptionBase
-from ..grid_data import BirkhoffGrid, GaussLobattoGrid, RadauGrid, UniformGrid
+from ..grid_data import BirkhoffGrid, GaussLobattoGrid, RadauGrid, UniformGrid, ChebyshevGaussLobattoGrid
 from .ode_integration_comp import ODEIntegrationComp
 from ...utils.misc import get_rate_units, CoerceDesvar
 from ...utils.indexing import get_src_indices_by_row
@@ -61,10 +61,12 @@ def initialize(self):
                                   '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, BirkhoffGrid, str), allow_none=True, default=None,
+        self.options.declare('grid', types=(GaussLobattoGrid, ChebyshevGaussLobattoGrid,
+                                            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, BirkhoffGrid), allow_none=True,
+        self.options.declare('output_grid', types=(GaussLobattoGrid, ChebyshevGaussLobattoGrid, 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 '
@@ -550,28 +552,6 @@ def configure_timeseries_outputs(self, phase):
         """
         super().configure_timeseries_outputs(phase)
 
-    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 for this transcription.
diff --git a/dymos/transcriptions/explicit_shooting/ode_integration_comp.py b/dymos/transcriptions/explicit_shooting/ode_integration_comp.py
index bf52150cc..a249fa82e 100644
--- a/dymos/transcriptions/explicit_shooting/ode_integration_comp.py
+++ b/dymos/transcriptions/explicit_shooting/ode_integration_comp.py
@@ -722,8 +722,8 @@ def _propagate(self, inputs, propagate_derivs=None, x_out=None, t_out=None, dx_d
             state_initial_val = inputs[self._state_input_names[state_name]]
             x0[self.state_idxs[state_name]] = state_initial_val.ravel()
 
-        theta[0] = t_initial = inputs['t_initial']
-        theta[1] = t_duration = inputs['t_duration']
+        theta[0] = t_initial = inputs['t_initial'].ravel()[0]
+        theta[1] = t_duration = inputs['t_duration'].ravel()[0]
 
         for param_name in self.parameter_options:
             param_val = inputs[self._param_input_names[param_name]]
diff --git a/dymos/transcriptions/grid_data.py b/dymos/transcriptions/grid_data.py
index eb668c4e9..fbdec6afe 100644
--- a/dymos/transcriptions/grid_data.py
+++ b/dymos/transcriptions/grid_data.py
@@ -38,6 +38,8 @@ def gauss_lobatto_subsets_and_nodes(n, seg_idx, compressed=False):
         'all' Gives all node indices.
     np.array
         The location of all nodes on [-1, 1].
+    np.array
+        The weights of all nodes on [-1, 1].
 
     Notes
     -----
@@ -60,7 +62,61 @@ def gauss_lobatto_subsets_and_nodes(n, seg_idx, compressed=False):
         'solution': np.arange(n, dtype=int),
     }
 
-    return subsets, lgl(n)[0]
+    return subsets, *lgl(n)
+
+
+def chebyshev_gauss_lobatto_subsets_and_nodes(n, seg_idx, compressed=False):
+    """
+    Provides node information and the location of the nodes for n Legendre-Gauss-Lobatto nodes on the range [-1, 1].
+
+    Parameters
+    ----------
+    n : int
+        The total number of nodes in the Gauss-Lobatto segment.  Must be
+        an odd number.
+    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].
+    np.array
+        The weights of all nodes on [-1, 1].
+
+    Notes
+    -----
+    Subset 'state_input' is the same as subset 'state_disc' if `compressed == False` or
+    `first_seg == True`.  The same is true of subsets 'control_input' and 'control_disc'.
+    """
+    if n < 2:
+        raise ValueError('The number of nodes must be larger than 1.')
+
+    subsets = {
+        'state_disc': np.arange(n, dtype=int),
+        'state_input': np.arange(n, dtype=int) if not compressed or seg_idx == 0
+        else np.arange(1, n, 2, dtype=int),
+        'control_disc': np.arange(n, dtype=int),
+        'control_input': np.arange(n, dtype=int) if not compressed or seg_idx == 0
+        else np.arange(1, n, dtype=int),
+        'segment_ends': np.array([0, n-1], dtype=int),
+        'col': np.arange(n, dtype=int),
+        'all': np.arange(n, dtype=int),
+        'solution': np.arange(n, dtype=int),
+    }
+
+    return subsets, *cgl(n)
 
 
 def radau_pseudospectral_subsets_and_nodes(n, seg_idx, compressed=False):
@@ -89,6 +145,8 @@ def radau_pseudospectral_subsets_and_nodes(n, seg_idx, compressed=False):
         'all' gives all node indices.
     np.array
         The location of all nodes on [-1, 1].
+    np.array
+        The weights of all nodes on [-1, 1].
 
     Notes
     -----
@@ -108,7 +166,7 @@ def radau_pseudospectral_subsets_and_nodes(n, seg_idx, compressed=False):
         'solution': np.arange(n + 1, dtype=int),
     }
 
-    return subsets, lgr(n, include_endpoint=True)[0]
+    return subsets, *lgr(n, include_endpoint=True)
 
 
 def birkhoff_subsets_and_nodes(n, grid, seg_idx, compressed=False):
@@ -139,6 +197,8 @@ def birkhoff_subsets_and_nodes(n, grid, seg_idx, compressed=False):
         'all' gives all node indices.
     np.array
         The location of all nodes on [-1, 1].
+    np.array
+        The weights of all nodes on [-1, 1].
 
     Notes
     -----
@@ -159,18 +219,18 @@ def birkhoff_subsets_and_nodes(n, grid, seg_idx, compressed=False):
     }
 
     if grid == 'lgl':
-        nodes = lgl(n)[0]
+        nodes, weights = lgl(n)
         subsets['all'] = np.arange(n, dtype=int)
     elif grid == 'lgr':
-        nodes = lgr(n, include_endpoint=False)[0]
+        nodes, weights = lgr(n, include_endpoint=False)
         subsets['all'] = np.arange(n, dtype=int)
     elif grid == 'cgl':
-        nodes = cgl(n)[0]
+        nodes, weights = cgl(n)
         subsets['all'] = np.arange(n, dtype=int)
     else:
         raise ValueError(f'Unrecognized grid. Acceptable values are one of {acceptable_grids}')
 
-    return subsets, nodes
+    return subsets, nodes, weights
 
 
 def uniform_subsets_and_nodes(n, *args, **kwargs):
@@ -216,7 +276,9 @@ def uniform_subsets_and_nodes(n, *args, **kwargs):
         'all': np.arange(n + 1, dtype=int),
         'solution': np.arange(n + 1, dtype=int),
     }
-    return subsets, np.linspace(-1, 1, n + 1)
+    weights = np.ones(n + 1)
+    weights[1:-1] = 2.0
+    return subsets, np.linspace(-1, 1, n + 1), weights
 
 
 def make_subset_map(from_subset_idxs, to_subset_idxs):
@@ -290,6 +352,8 @@ class GridData(object):
         The number of steps to take in each segment of the phase, for explicit phases.
     num_nodes : int
         The total number of nodes in the phase
+    node_weight : ndarray
+        The quadrature weight for each node on the interval [-1, 1].
     node_stau : ndarray
         The locations of each node in non-dimensional segment time (segment tau space).
     node_ptau : ndarray
@@ -343,6 +407,8 @@ def __init__(self, num_segments, transcription, transcription_order=None,
 
         self.num_steps_per_segment = 0
 
+        self.node_weight = np.empty(0,)
+
         self.node_stau = np.empty(0,)
 
         self.node_ptau = np.empty(0,)
@@ -364,10 +430,12 @@ def __init__(self, num_segments, transcription, transcription_order=None,
         self.input_maps = {'state_input_to_disc': np.empty(0, dtype=int),
                            'dynamic_control_input_to_disc': np.empty(0, dtype=int)}
 
-        if transcription.lower() in ['radau', 'radau-ps']:
+        if transcription.lower() in ['radau', 'radau-ps', 'lgr']:
             self.transcription = 'radau-ps'
         elif transcription.lower() in ['gausslobatto', 'gauss-lobatto', 'lgl']:
             self.transcription = 'gauss-lobatto'
+        elif transcription.lower() in ['chebyshev-gauss-lobatto', 'cgl']:
+            self.transcription = 'chebyshev-gauss-lobatto'
         elif transcription.lower() in ['birkhoff']:
             self.transcription = 'birkhoff'
         elif transcription.lower() in ['uniform']:
@@ -378,6 +446,8 @@ def __init__(self, num_segments, transcription, transcription_order=None,
         # Define get_subsets and node points based on the transcription scheme
         if self.transcription == 'gauss-lobatto':
             get_subsets_and_nodes = gauss_lobatto_subsets_and_nodes
+        if self.transcription == 'chebyshev-gauss-lobatto':
+            get_subsets_and_nodes = chebyshev_gauss_lobatto_subsets_and_nodes
         elif self.transcription == 'radau-ps':
             get_subsets_and_nodes = radau_pseudospectral_subsets_and_nodes
         elif self.transcription == 'uniform':
@@ -403,9 +473,9 @@ def __init__(self, num_segments, transcription, transcription_order=None,
         self.segment_indices[0, 0] = 0
         ind0 = 0  # index of the first node in the segment
         for iseg in range(num_segments):
-            subsets_i, nodes_i = get_subsets_and_nodes(self.transcription_order[iseg],
-                                                       seg_idx=iseg,
-                                                       compressed=compressed)
+            subsets_i, nodes_i, weights_i = get_subsets_and_nodes(self.transcription_order[iseg],
+                                                                  seg_idx=iseg,
+                                                                  compressed=compressed)
 
             if iseg == 0:
                 subset_ind0 = {name: 0 for name in subsets_i}
@@ -419,6 +489,9 @@ def __init__(self, num_segments, transcription, transcription_order=None,
             # Append our nodes in segment tau space
             self.node_stau = np.concatenate((self.node_stau, nodes_i))
 
+            # Append our node weights
+            self.node_weight = np.concatenate((self.node_weight, weights_i))
+
             # Append our nodes in phase tau space
             v0 = segment_ends[iseg]
             v1 = segment_ends[iseg + 1]
@@ -491,6 +564,9 @@ def __eq__(self, other):
         else:
             return False
 
+    def __repr__(self):
+        return f'{self.__class__.__name__}(num_seg={self.num_segments}, order={self.transcription_order}) at <{id(self)}>'
+
     def is_aligned_with(self, other, tol=1.0E-12):
         """
         Check that the segment distribution in GridData object `other` matches that of this GridData object.
@@ -673,11 +749,38 @@ class GaussLobattoGrid(GridData):
         to the appropriate indices.
     """
     def __init__(self, num_segments, nodes_per_seg, segment_ends=None, compressed=False):
+        self.grid_type = 'lgl'
         super().__init__(num_segments=num_segments, transcription='gauss-lobatto',
                          transcription_order=np.asarray(nodes_per_seg, dtype=int),
                          segment_ends=segment_ends, compressed=compressed)
 
 
+class ChebyshevGaussLobattoGrid(GridData):
+    """
+    A GridData object that provides the node information for a Gauss-Lobatto distribution.
+
+    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.
+    """
+    def __init__(self, num_segments, nodes_per_seg, segment_ends=None, compressed=False):
+        self.grid_type = 'cgl'
+        super().__init__(num_segments=num_segments, transcription='chebyshev-gauss-lobatto',
+                         transcription_order=np.asarray(nodes_per_seg, dtype=int),
+                         segment_ends=segment_ends, compressed=compressed)
+
+
 class BirkhoffGrid(GridData):
     """
     A GridData object that provides the node information for the Birkhoff transcription.
@@ -715,6 +818,7 @@ class RadauGrid(GridData):
         to the appropriate indices.
     """
     def __init__(self, num_segments, nodes_per_seg, segment_ends=None, compressed=False):
+        self.grid_type = 'lgr'
         super().__init__(num_segments=num_segments, transcription='radau-ps',
                          transcription_order=np.asarray(nodes_per_seg, dtype=int) - 1,
                          segment_ends=segment_ends, compressed=compressed)
@@ -740,6 +844,7 @@ class UniformGrid(GridData):
         to the appropriate indices.
     """
     def __init__(self, num_segments, nodes_per_seg, segment_ends=None, compressed=False):
+        self.grid_type = 'uniform'
         super().__init__(num_segments=num_segments, transcription='uniform',
                          transcription_order=np.asarray(nodes_per_seg) - 1,
                          segment_ends=segment_ends, compressed=compressed)
diff --git a/dymos/transcriptions/picard_shooting/__init__.py b/dymos/transcriptions/picard_shooting/__init__.py
new file mode 100644
index 000000000..e69de29bb
diff --git a/dymos/transcriptions/picard_shooting/birkhoff_picard_iter_group.py b/dymos/transcriptions/picard_shooting/birkhoff_picard_iter_group.py
new file mode 100644
index 000000000..ac39cc150
--- /dev/null
+++ b/dymos/transcriptions/picard_shooting/birkhoff_picard_iter_group.py
@@ -0,0 +1,113 @@
+import openmdao.api as om
+
+from .birkhoff_picard_update_comp import PicardUpdateComp
+from .states_comp import StatesComp
+
+from ..grid_data import GridData
+from ...utils.ode_utils import _make_ode_system
+
+
+class BirkhoffPicardIterGroup(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 __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        self._implicit_outputs = set()
+
+    def initialize(self):
+        """
+        Declare group options.
+        """
+        self.options.declare('state_options', types=dict,
+                             desc='Dictionary of options for the states.')
+        self.options.declare('time_units', types=str,
+                             desc='Units 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')
+        self.options.declare('calc_exprs', types=dict, default={},
+                             desc='phase calculation expressions.')
+        self.options.declare('parameter_options', types=dict, default={},
+                             desc='phase parameter options')
+
+    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_units = self.options['time_units']
+
+        self.add_subsystem('states_comp',
+                           subsys=StatesComp(grid_data=gd, state_options=state_options))
+
+        ode = _make_ode_system(ode_class=self.options['ode_class'],
+                               num_nodes=nn,
+                               calc_exprs=self.options['calc_exprs'],
+                               ode_init_kwargs=self.options['ode_init_kwargs'],
+                               parameter_options=self.options['parameter_options'])
+
+        self.add_subsystem('ode_all', subsys=ode)
+
+        self.add_subsystem('picard_update_comp',
+                           subsys=PicardUpdateComp(grid_data=gd,
+                                                   state_options=state_options,
+                                                   time_units=time_units))
+
+    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.
+        """
+        states_comp = self._get_subsystem('states_comp')
+        states_comp.configure_io(phase)
+
+        self.promotes('states_comp', inputs=['*'], outputs=['*'])
+
+        picard_update_comp = self._get_subsystem('picard_update_comp')
+        picard_update_comp.configure_io(phase)
+
+        state_options = self.options['state_options']
+
+        for name, options in state_options.items():
+            rate_source = options['rate_source']
+
+            for tgt in options['targets']:
+                self.connect(f'state_val:{name}', f'ode_all.{tgt}')
+
+            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'picard_update_comp.f_computed:{name}')
+
+            promotes = [f'states:{name}']
+            if options['solve_segments'] == 'forward':
+                promotes += [f'final_states:{name}']
+            else:
+                promotes += [f'initial_states:{name}']
+
+            self.promotes('picard_update_comp', any=promotes)
diff --git a/dymos/transcriptions/picard_shooting/birkhoff_picard_update_comp.py b/dymos/transcriptions/picard_shooting/birkhoff_picard_update_comp.py
new file mode 100644
index 000000000..859599cb9
--- /dev/null
+++ b/dymos/transcriptions/picard_shooting/birkhoff_picard_update_comp.py
@@ -0,0 +1,269 @@
+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.birkhoff import birkhoff_matrix
+
+
+class PicardUpdateComp(om.ExplicitComponent):
+    """
+    Class definition for the PicardUpdateComp.
+
+    Given the initial state values (for forward propagation) or the final state values
+    (for backward propagation), compute the next state value for picard iteration
+    using a NonlinearBlockGS solver.
+
+    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['all']
+        num_segs = gd.num_segments
+        time_units = self.options['time_units']
+        state_options = self.options['state_options']
+
+        B_blocks = []
+
+        self._seg_repeats = gd.subset_num_nodes_per_segment['all']
+
+        start_idx = 0
+        for i in range(num_segs):
+            nnps_i = gd.subset_num_nodes_per_segment['all'][i]
+            tau_i = gd.node_stau[start_idx: start_idx + nnps_i]
+            w_i = gd.node_weight[start_idx: start_idx + nnps_i]
+
+            B_i = birkhoff_matrix(tau_i, w_i, grid_type=gd.grid_type)
+            B_blocks.append(B_i)
+            start_idx += nnps_i
+
+        self._B = scipy.sparse.block_diag(B_blocks, format='csr')
+
+        self.add_input('dt_dstau', units=self.options['time_units'], shape=(num_nodes,))
+
+        self.var_names = var_names = {}
+        for state_name, options in state_options.items():
+            var_names[state_name] = {
+                'f_computed': f'f_computed:{state_name}',
+                'x_0': f'seg_initial_states:{state_name}',
+                'x_f': f'seg_final_states:{state_name}',
+                'x_a': f'initial_states:{state_name}',
+                'x_b': f'final_states:{state_name}',
+                'x_hat': f'states:{state_name}'
+            }
+
+        for state_name, options in state_options.items():
+            shape = options['shape']
+            size = np.prod(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)
+
+            if rate_source_type == 'state':
+                var_names['f_computed'] = f'state_val:{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_output(
+                name=var_names['x_hat'],
+                shape=(num_nodes,) + shape,
+                units=units
+            )
+
+            if options['solve_segments'] == 'forward':
+                self.add_input(
+                    name=var_names['x_0'],
+                    shape=(num_segs,) + shape,
+                    desc=f'Initial value of state {state_name} in each segment',
+                    units=units
+                )
+                self.add_output(
+                    name=var_names['x_b'],
+                    shape=(1,) + shape,
+                    desc=f'Final value of state {state_name} in the phase',
+                    units=units
+                )
+                # The first row of these matrices are zero.
+                rs = np.repeat(np.arange(1, num_nodes, dtype=int), num_nodes)
+                cs = np.tile(np.arange(num_nodes, dtype=int), num_nodes - 1)
+                self.declare_partials(of=var_names['x_hat'],
+                                      wrt='dt_dstau',
+                                      rows=rs, cols=cs)
+                self.declare_partials(of=var_names['x_hat'],
+                                      wrt=var_names['f_computed'],
+                                      rows=rs, cols=cs)
+                self.declare_partials(of=var_names['x_b'],
+                                      wrt='dt_dstau')
+                self.declare_partials(of=var_names['x_b'],
+                                      wrt=var_names['f_computed'])
+
+                rs = np.arange(num_nodes, dtype=int)
+                cs = np.repeat(np.arange(num_segs, dtype=int), self._seg_repeats)
+                self.declare_partials(of=var_names['x_hat'],
+                                      wrt=var_names['x_0'],
+                                      rows=rs, cols=cs, val=1.0)
+                template = sp.lil_array((1, num_segs), dtype=int)
+                template[0, -1] = 1
+                template = sp.kron(template.tocsr(), sp.eye(size, dtype=int))
+                rs, cs = template.nonzero()
+                self.declare_partials(of=var_names['x_b'],
+                                      wrt=var_names['x_0'],
+                                      rows=rs, cols=cs, val=1.0)
+
+            elif options['solve_segments'] == 'backward':
+                self.add_input(name=var_names['x_f'],
+                               shape=(num_segs,) + shape,
+                               desc=f'Final value of state {state_name} in each segment',
+                               units=units)
+                self.add_output(name=var_names['x_a'],
+                                shape=(1,) + shape,
+                                desc=f'Initial value of state {state_name} in the phase',
+                                units=units)
+                rs = np.repeat(np.arange(num_nodes - 1, dtype=int), num_nodes)
+                cs = np.tile(np.arange(num_nodes, dtype=int), num_nodes - 1)
+                self.declare_partials(of=var_names['x_hat'], wrt='dt_dstau',
+                                      rows=rs, cols=cs)
+                self.declare_partials(of=var_names['x_hat'], wrt=var_names['f_computed'],
+                                      rows=rs, cols=cs)
+                self.declare_partials(of=var_names['x_a'],
+                                      wrt='dt_dstau')
+                self.declare_partials(of=var_names['x_a'],
+                                      wrt=var_names['f_computed'])
+                rs = np.arange(num_nodes, dtype=int)
+                cs = np.repeat(np.arange(num_segs, dtype=int), self._seg_repeats)
+                self.declare_partials(of=var_names['x_hat'],
+                                      wrt=var_names['x_f'],
+                                      rows=rs, cols=cs, val=1.0)
+                template = sp.lil_array((1, num_segs), dtype=int)
+                template[0, 0] = 1
+                template = sp.kron(template.tocsr(), sp.eye(size, dtype=int))
+                rs, cs = template.nonzero()
+                self.declare_partials(of=var_names['x_a'], wrt=var_names['x_f'],
+                                      rows=rs, cols=cs, 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.
+        """
+        dt_dstau = np.atleast_2d(inputs['dt_dstau']).T
+        nn = self.options['grid_data'].num_nodes
+
+        for state_name, options in self.options['state_options'].items():
+            var_names = self.var_names[state_name]
+
+            # Multiplication by B results in the integral in segment tau space,
+            # so we need to convert f from being dx/dt to dx/dstau.
+            f = np.einsum('i...,i...->i...', inputs[var_names['f_computed']], dt_dstau)
+            f_flat = f.reshape(nn, -1)
+
+            if options['solve_segments'] == 'forward':
+                x_0 = inputs[var_names['x_0']]
+                x_0_repeated = np.repeat(x_0, self._seg_repeats, axis=0)
+                outputs[var_names['x_hat']] = x_0_repeated + (self._B @ f_flat).reshape(f.shape)
+                outputs[var_names['x_b']][...] = outputs[var_names['x_hat']][-1, ...]
+
+            elif options['solve_segments'] == 'backward':
+                x_f = inputs[var_names['x_f']]
+                x_f_repeated = np.repeat(x_f, self._seg_repeats, axis=0)
+                outputs[var_names['x_hat']] = x_f_repeated - (self._B[::-1, ...] @ f_flat[::-1, ...]).reshape(f.shape)
+                outputs[var_names['x_a']][...] = outputs[var_names['x_hat']][0, ...]
+
+    def compute_partials(self, inputs, partials, discrete_inputs=None):
+        """
+        Compute component partials.
+
+        Parameters
+        ----------
+        inputs : Vector
+            Unscaled, dimensional input variables read via inputs[key].
+        partials : dict
+            Partial values to be returend in partials[of, wrt] format.
+        discrete_inputs : dict or None
+            If not None, dict containing discrete input values.
+        """
+        gd = self.options['grid_data']
+        num_nodes = gd.subset_num_nodes['all']
+
+        dt_dstau = np.atleast_2d(inputs['dt_dstau']).T
+
+        for state_name, options in self.options['state_options'].items():
+            var_names = self.var_names[state_name]
+
+            x_a_name = var_names['x_a']
+            x_b_name = var_names['x_b']
+            x_name = var_names['x_hat']
+            f_name = var_names['f_computed']
+
+            # Multiplication by B results in the integral in tau space,
+            # so we need to convert f from being dx/dt to dx/dtau.
+            f_t = inputs[var_names['f_computed']]
+
+            if options['solve_segments'] == 'forward':
+                partials[x_name, f_name] = (self._B.multiply(dt_dstau.T)).todense()[1:, ...].ravel()
+                partials[x_b_name, f_name] = partials[x_name, f_name][-num_nodes:]
+                partials[x_name, 'dt_dstau'] = (self._B.multiply(f_t.T)).todense()[1:, ...].ravel()
+                partials[x_b_name, 'dt_dstau'] = partials[x_name, 'dt_dstau'][-num_nodes:]
+
+            elif options['solve_segments'] == 'backward':
+                B_flip = self._B[::-1, ...]
+                dt_dstau_flip = dt_dstau[::-1, ...]
+
+                partials[x_name, f_name] = -(B_flip.multiply(dt_dstau_flip.T)).todense()[:-1, ::-1].ravel()
+                partials[x_a_name, f_name] = partials[x_name, f_name][:num_nodes]
+                partials[x_name, 'dt_dstau'] = -(B_flip.multiply(f_t[::-1, ...].T)).todense()[:-1, ::-1].ravel()
+                partials[x_a_name, 'dt_dstau'] = partials[x_name, 'dt_dstau'][:num_nodes]
diff --git a/dymos/transcriptions/picard_shooting/multiple_shooting_iter_group.py b/dymos/transcriptions/picard_shooting/multiple_shooting_iter_group.py
new file mode 100644
index 000000000..e1c907dc1
--- /dev/null
+++ b/dymos/transcriptions/picard_shooting/multiple_shooting_iter_group.py
@@ -0,0 +1,107 @@
+import openmdao.api as om
+
+from ..grid_data import GridData
+from .birkhoff_picard_iter_group import BirkhoffPicardIterGroup
+from .multiple_shooting_update_comp import MultipleShootingUpdateComp
+
+
+class MultipleShootingIterGroup(om.Group):
+    """
+    Class definition for the MultipleShootingIterGroup.
+
+    This group allows for iteration of the initial or final state variable values in each segment
+    (in forward shooting or backward shooting, respectively). This group feeds back new guesses
+    for these "starting" state values based on the results of integration the causal segment
+
+    Parameters
+    ----------
+    **kwargs : dict
+        Dictionary of optional arguments.
+    """
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+
+    def initialize(self):
+        """
+        Declare group options.
+        """
+        self.options.declare('state_options', types=dict,
+                             desc='Dictionary of options for the states.')
+        self.options.declare('time_units', types=str,
+                             desc='Units 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')
+        self.options.declare('ode_nonlinear_solver', default=om.NonlinearBlockGS(maxiter=100, use_aitken=True, iprint=0),
+                             desc='Nonlinear solver used to resolve Picard iteration.', recordable=False)
+        self.options.declare('ode_linear_solver', default=om.DirectSolver(),
+                             desc='Linear solver used to linearize the picard iteration subsystem.', recordable=False)
+        self.options.declare('ms_nonlinear_solver', default=om.NonlinearBlockGS(maxiter=100, use_aitken=True, iprint=0),
+                             desc='Nonlinear solver used to resolve Picard iteration.', recordable=False)
+        self.options.declare('ms_linear_solver', default=om.DirectSolver(),
+                             desc='Linear solver used to linearize the picard iteration subsystem.', recordable=False)
+        self.options.declare('calc_exprs', types=dict, default={},
+                             desc='phase calculation expressions.')
+        self.options.declare('parameter_options', types=dict, default={},
+                             desc='phase parameter options')
+
+    def setup(self):
+        """
+        Define the structure of the control group.
+        """
+        gd = self.options['grid_data']
+        state_options = self.options['state_options']
+        time_units = self.options['time_units']
+        ode_class = self.options['ode_class']
+        ode_init_kwargs = self.options['ode_init_kwargs']
+        calc_exprs = self.options['calc_exprs']
+        parameter_options = self.options['parameter_options']
+
+        seg_prop_group = self.add_subsystem('segment_prop_group',
+                                            subsys=BirkhoffPicardIterGroup(grid_data=gd,
+                                                                           state_options=state_options,
+                                                                           time_units=time_units,
+                                                                           ode_class=ode_class,
+                                                                           ode_init_kwargs=ode_init_kwargs,
+                                                                           calc_exprs=calc_exprs,
+                                                                           parameter_options=parameter_options),
+                                            promotes_inputs=['*'], promotes_outputs=['*'])
+
+        seg_prop_group.nonlinear_solver = self.options['ode_nonlinear_solver']
+        seg_prop_group.linear_solver = self.options['ode_linear_solver']
+
+        self.nonlinear_solver = self.options['ms_nonlinear_solver']
+        self.linear_solver = self.options['ms_linear_solver']
+
+        self.add_subsystem('ms_update_comp',
+                           MultipleShootingUpdateComp(grid_data=gd,
+                                                      state_options=state_options),
+                           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.
+        """
+        state_options = self.options['state_options']
+
+        segment_prop_group = self._get_subsystem('segment_prop_group')
+        segment_prop_group.configure_io(phase)
+
+        ms_update_comp = self._get_subsystem('ms_update_comp')
+        ms_update_comp.configure_io(phase)
+
+        for state_name, options in state_options.items():
+            if options['solve_segments'] == 'forward':
+                self.connect(f'seg_initial_states:{state_name}',
+                             f'picard_update_comp.seg_initial_states:{state_name}')
+            else:
+                self.connect(f'seg_final_states:{state_name}',
+                             f'picard_update_comp.seg_final_states:{state_name}')
diff --git a/dymos/transcriptions/picard_shooting/multiple_shooting_update_comp.py b/dymos/transcriptions/picard_shooting/multiple_shooting_update_comp.py
new file mode 100644
index 000000000..b3f5e2cf2
--- /dev/null
+++ b/dymos/transcriptions/picard_shooting/multiple_shooting_update_comp.py
@@ -0,0 +1,216 @@
+import numpy as np
+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
+
+
+class MultipleShootingUpdateComp(om.ExplicitComponent):
+    """
+    Class definition for the PicardUpdateComp.
+
+    Given the initial state values (for forward propagation) or the final state values
+    (for backward propagation), compute the next state value for picard iteration
+    using a NonlinearBlockGS solver.
+
+    Parameters
+    ----------
+    **kwargs : dict
+        Dictionary of optional arguments.
+
+    Attributes
+    ----------
+    _no_check_partials : bool
+        If True, this component will be excluded when partials are checked.
+    _M_fwd : array-like
+        A selection matrix that pulls appropriate values from x into x_0
+        when in forward propagation mode.
+    _M_bkwd : array-like
+        A selection matrix that pulls appropriate values from x into x_f
+        when in backwards propagation mode.
+    """
+
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        self._no_check_partials = not dymos_options['include_check_partials']
+        self._M_fwd = None
+        self._M_bkwd = None
+
+    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['all']
+        num_segs = gd.num_segments
+        time_units = self.options['time_units']
+        state_options = self.options['state_options']
+
+        # Construct the forward and backward mapping matrices
+        seg_start_nodes = gd.subset_node_indices['segment_ends'][::2]
+        seg_end_nodes = gd.subset_node_indices['segment_ends'][1::2]
+
+        M_fwd = sp.lil_array((num_segs, num_nodes), dtype=int)
+        for iseg in range(1, gd.num_segments):
+            # The ith row of M_fwd contains in the column pertaining to first node in the (i-1)th segment
+            M_fwd[iseg, seg_end_nodes[iseg-1]] = 1
+
+        M_bkwd = sp.lil_array((num_segs, num_nodes), dtype=int)
+        for iseg in range(gd.num_segments-1):
+            # The ith row of M_bkwd contains in the column pertaining to first node in the (i+1)th segment
+            M_bkwd[iseg, seg_start_nodes[iseg+1]] = 1
+
+        self._M_fwd = {}
+        self._M_bkwd = {}
+
+        self._var_names = var_names = {}
+        for state_name, options in state_options.items():
+            var_names[state_name] = {
+                'x_a': f'initial_states:{state_name}',
+                'x_b': f'final_states:{state_name}',
+                'x_0': f'seg_initial_states:{state_name}',
+                'x_f': f'seg_final_states:{state_name}',
+                'x': f'states:{state_name}'
+            }
+
+        for state_name, options in state_options.items():
+            shape = options['shape']
+            units = options['units']
+            size = np.prod(shape)
+
+            rate_units = get_rate_units(units, time_units)
+            var_names = self._var_names[state_name]
+
+            self.add_input(
+                name=var_names['x'],
+                shape=(num_nodes,) + shape,
+                desc='State value at the nodes.',
+                units=rate_units)
+
+            if options['solve_segments'] == 'forward':
+                self._M_fwd[state_name] = sp.kron(M_fwd, sp.eye(size, dtype=int))
+
+                self.add_input(
+                    name=var_names['x_a'],
+                    shape=(1,) + shape,
+                    desc=f'Desired initial value of state {state_name} in phase.',
+                    units=units
+                )
+                self.add_output(
+                    name=var_names['x_0'],
+                    shape=(num_segs,) + shape,
+                    val=0.0,
+                    desc=f'Given initial value of state {state_name} in each segment.',
+                    units=units
+                )
+
+                ar_size_x_nn0 = np.arange(size, dtype=int)
+                self.declare_partials(of=var_names['x_0'],
+                                      wrt=var_names['x_a'],
+                                      rows=ar_size_x_nn0,
+                                      cols=np.zeros_like(ar_size_x_nn0),
+                                      val=1.0)
+
+                if num_segs > 1:
+                    rs, cs = self._M_fwd[state_name].nonzero()
+                    self.declare_partials(of=var_names['x_0'],
+                                          wrt=var_names['x'],
+                                          rows=rs, cols=cs,
+                                          val=self._M_fwd[state_name].data.ravel())
+
+            elif options['solve_segments'] == 'backward':
+                self._M_bkwd[state_name] = sp.kron(M_bkwd, sp.eye(size, dtype=int))
+
+                self.add_input(
+                    name=var_names['x_b'],
+                    shape=(1,) + shape,
+                    desc=f'Desired initial value of state {state_name} in phase.',
+                    units=units
+                )
+
+                self.add_output(
+                    name=var_names['x_f'],
+                    shape=(num_segs,) + shape,
+                    val=0.0,
+                    desc=f'Given final value of state {state_name} in each segment.',
+                    units=units
+                )
+
+                ar_size_x_nnf = (num_segs - 1) * size + np.arange(size, dtype=int)
+                self.declare_partials(of=var_names['x_f'],
+                                      wrt=var_names['x_b'],
+                                      rows=ar_size_x_nnf,
+                                      cols=np.zeros_like(ar_size_x_nnf),
+                                      val=1.0)
+
+                if num_segs > 1:
+                    rs, cs = self._M_bkwd[state_name].nonzero()
+                    self.declare_partials(of=var_names['x_f'],
+                                          wrt=var_names['x'],
+                                          rows=rs, cols=cs,
+                                          val=self._M_bkwd[state_name].data.ravel())
+
+    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.
+        """
+        gd = self.options['grid_data']
+        nn = gd.num_nodes
+
+        for state_name, options in self.options['state_options'].items():
+            var_names = self._var_names[state_name]
+            x = inputs[var_names['x']]
+            x_flat = x.reshape(nn, -1)
+
+            if options['solve_segments'] == 'forward':
+                M_fwd = self._M_fwd[state_name]
+                x_a = inputs[var_names['x_a']]
+                outputs[var_names['x_0']] = M_fwd @ x_flat
+                outputs[var_names['x_0']][0, ...] = x_a[0, ...]
+                # print(x_flat)
+                # print(x_a)
+                # print(outputs[var_names['x_0']])
+
+            elif options['solve_segments'] == 'backward':
+                M_bkwd = self._M_bkwd[state_name]
+                x_b = inputs[var_names['x_b']]
+                outputs[var_names['x_f']] = M_bkwd @ x_flat
+                outputs[var_names['x_f']][-1, ...] = x_b
+
+            else:
+                raise ValueError(f'{self.msginfo}: Invalid direction of integration: {options["solve_segments"]}')
diff --git a/dymos/transcriptions/picard_shooting/picard_shooting.py b/dymos/transcriptions/picard_shooting/picard_shooting.py
new file mode 100644
index 000000000..158c32be1
--- /dev/null
+++ b/dymos/transcriptions/picard_shooting/picard_shooting.py
@@ -0,0 +1,871 @@
+import numpy as np
+
+import openmdao.api as om
+
+from ..transcription_base import TranscriptionBase
+from ..common import TimeComp, TimeseriesOutputComp, ControlInterpComp
+from .multiple_shooting_iter_group import MultipleShootingIterGroup
+
+from ..grid_data import GaussLobattoGrid, ChebyshevGaussLobattoGrid
+from dymos.utils.introspection import get_promoted_vars, get_source_metadata, get_rate_units
+from dymos.utils.indexing import get_constraint_flat_idxs, get_src_indices_by_row
+from dymos.utils.misc import _format_phase_constraint_alias
+
+
+class PicardShooting(TranscriptionBase):
+    """
+        Picard Shooting transcription.
+
+        Parameters
+        ----------
+        **kwargs : dict
+            Dictionary of optional arguments.
+
+        References
+        ----------
+        TBD
+        """
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+        self._rhs_source = 'ode_iter_group.segment_prop_group.ode_all'
+
+    def initialize(self):
+        """
+        Declare transcription options.
+        """
+        self.options.declare('num_segments', types=int, default=1,
+                             desc='The number of segments in the grid.')
+
+        self.options.declare('nodes_per_seg', types=int, default=10,
+                             desc='The number of segments in the grid.')
+
+        self.options.declare('grid_type', values=('cgl', 'lgl'), default='cgl',
+                             desc='Specifies which type of grid is to be used. '
+                                  'Options are Chebyshev-Gauss-Lobatto ("cgl") '
+                                  'and Legendre-Gauss-Lobatto ("lgl")')
+
+        self.options.declare(name='solve_segments', default='forward',
+                             values=('forward', 'backward'),
+                             desc='The default solve direction for states in this phase.'
+                             'This value may be overridden by setting the solve_segments option of states')
+
+        self.options.declare('ode_nonlinear_solver', default=om.NonlinearBlockGS(maxiter=100, use_aitken=True, iprint=0),
+                             desc='Nonlinear solver used to resolve Picard iteration on each segment.',
+                             recordable=False)
+
+        self.options.declare('ode_linear_solver', default=om.DirectSolver(),
+                             desc='Linear solver used to linearize the Picard iteration subsystem on each segment.',
+                             recordable=False)
+
+        self.options.declare('ms_nonlinear_solver', default=om.NonlinearBlockGS(maxiter=100, use_aitken=True, iprint=0),
+                             desc='Nonlinear solver used to resolve Picard iteration differences between segments.',
+                             recordable=False)
+
+        self.options.declare('ms_linear_solver', default=om.DirectSolver(),
+                             desc='Linear solver used to linearize the Picard iteration subsystem between segments.',
+                             recordable=False)
+
+    def init_grid(self):
+        """
+        Setup the GridData object for the Transcription.
+        """
+        if self.options['grid_type'] == 'cgl':
+            grid_cls = ChebyshevGaussLobattoGrid
+        elif self.options['grid_type'] == 'lgl':
+            grid_cls = GaussLobattoGrid
+        self.grid_data = grid_cls(num_segments=self.options['num_segments'],
+                                  nodes_per_seg=self.options['nodes_per_seg'])
+
+    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().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=[('t_initial', 't_initial_val'), ('t_duration', 't_duration_val')],
+                            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().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)
+
+        for name, targets in [('t_initial', options['t_initial_targets']),
+                              ('t_duration', options['t_duration_targets']),
+                              ('t_final', options['t_final_targets'])]:
+            for t in targets:
+                shape = ode_inputs[t]['shape']
+
+                if shape == (1,):
+                    src_idxs = None
+                    flat_src_idxs = None
+                    src_shape = None
+                else:
+                    src_idxs = np.zeros(self.grid_data.subset_num_nodes['all'], 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)
+            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.
+        """
+        self.any_solved_segs = False
+        self.any_connected_opt_segs = False
+        for options in phase.state_options.values():
+            # Transcription solve_segments overrides state solve_segments if its not set
+            if options['solve_segments'] in (None, False):
+                options['solve_segments'] = self.options['solve_segments']
+
+            if options['solve_segments']:
+                self.any_solved_segs = True
+            elif options['input_initial']:
+                self.any_connected_opt_segs = True
+
+    def setup_controls(self, phase):
+        """
+        Setup the control group.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        phase._check_control_options()
+
+        if phase.control_options:
+            control_comp = ControlInterpComp(control_options=phase.control_options,
+                                             time_units=phase.time_options['units'],
+                                             grid_data=self.grid_data)
+
+            phase.add_subsystem('control_comp',
+                                subsys=control_comp)
+
+            phase.connect('t_duration_val', 'control_comp.t_duration')
+
+    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_comp.configure_io()
+            phase.promotes('control_comp',
+                           any=['*controls:*', '*control_values:*', '*control_rates:*'])
+
+            phase.connect('dt_dstau', 'control_comp.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']])
+
+            if options['rate_targets']:
+                phase.connect(f'control_rates:{name}_rate',
+                              [f'ode_all.{t}' for t in options['rate_targets']])
+
+            if options['rate2_targets']:
+                phase.connect(f'control_rates:{name}_rate2',
+                              [f'ode_all.{t}' for t in options['rate2_targets']])
+
+    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_nonlinear_solver = self.options['ode_nonlinear_solver']
+        ode_linear_solver = self.options['ode_linear_solver']
+        ms_nonlinear_solver = self.options['ms_nonlinear_solver']
+        ms_linear_solver = self.options['ms_linear_solver']
+
+        ode_init_kwargs = phase.options['ode_init_kwargs']
+        calc_exprs = phase._calc_exprs
+        parameter_options = phase.parameter_options
+
+        phase.add_subsystem('ode_iter_group',
+                            subsys=MultipleShootingIterGroup(grid_data=grid_data,
+                                                             state_options=phase.state_options,
+                                                             time_units=phase.time_options['units'],
+                                                             ode_class=ODEClass,
+                                                             ode_init_kwargs=ode_init_kwargs,
+                                                             ode_nonlinear_solver=ode_nonlinear_solver,
+                                                             ode_linear_solver=ode_linear_solver,
+                                                             ms_nonlinear_solver=ms_nonlinear_solver,
+                                                             ms_linear_solver=ms_linear_solver,
+                                                             calc_exprs=calc_exprs,
+                                                             parameter_options=parameter_options),
+                            promotes_inputs=['*'], promotes_outputs=['*'])
+
+    def configure_ode(self, phase):
+        """
+        Create connections to the introspected states.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        phase._get_subsystem('ode_iter_group').configure_io(phase)
+        phase.connect('dt_dstau', 'picard_update_comp.dt_dstau')
+
+    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):
+        """
+        Connect the collocation constraints.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        """
+        num_nodes = self.grid_data.subset_num_nodes['all']
+        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_src_path.startswith('parameter_vals:'):
+                src_idxs = om.slicer[np.zeros(num_nodes, dtype=int), ...]
+            else:
+                src_idxs = None
+
+            if rate_source_type not in ('state', 'ode'):
+                phase.connect(rate_src_path, f'picard_update_comp.f_computed:{name}', src_indices=src_idxs)
+
+    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, requires_solvers=None):
+        """
+        Configure the solvers.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase object to which this transcription instance applies.
+        requires_solvers : dict[str: bool]
+            A dictionary mapping a string descriptor of a reason why a solver is required,
+            and whether a solver is required.
+        """
+        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(timeseries_name)
+            ts_inputs_to_promote = []
+            for input_name, src, src_idxs in timeseries_comp._configure_io(timeseries_options):
+                # If the src was added, promote it if it was a state,
+                # or connect it otherwise.
+                if src.startswith('states:'):
+                    state_name = src.split(':')[-1]
+                    ts_inputs_to_promote.append((input_name, f'states:{state_name}'))
+                else:
+                    phase.connect(src_name=src,
+                                  tgt_name=f'{timeseries_name}.{input_name}',
+                                  src_indices=src_idxs)
+            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():
+            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'])
+            phase.add_subsystem(name, subsys=timeseries_comp)
+
+            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()}
+
+        # Determine the path to the variable which we will be constraining
+        var = 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'] = (num_nodes - 1) * size + flat_idxs
+            else:
+                # Path
+                path_idxs = []
+                for i in range(num_nodes):
+                    path_idxs.extend(size * i + flat_idxs)
+
+                constraint_kwargs['indices'] = path_idxs
+
+        con_path = constraint_kwargs.pop('constraint_path')
+        con_name = constraint_kwargs.pop('constraint_name')
+
+        constraint_kwargs['alias'] = _format_phase_constraint_alias(phase, con_name,
+                                                                    constraint_type,
+                                                                    options['indices'])
+        constraint_kwargs.pop('name')
+        constraint_kwargs.pop('shape')
+        constraint_kwargs['flat_indices'] = True
+
+        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']
+            linear = False
+            constraint_path = f'timeseries.{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.control_options[var]['shape']
+            units = phase.control_options[var]['units']
+            linear = True
+            constraint_path = f'control_values:{var}'
+        elif var_type == 'input_polynomial_control':
+            shape = phase.control_options[var]['shape']
+            units = phase.control_options[var]['units']
+            linear = False
+            constraint_path = f'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.control_options[control_var]['shape']
+            control_units = phase.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'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'timeseries.{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.
+        """
+        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'control_values:{var}'
+        elif var_type == 'input_polynomial_control':
+            rate_path = f'control_values:{var}'
+        elif var_type == 'polynomial_control_rate':
+            control_name = var[:-5]
+            rate_path = f'control_rates:{control_name}_rate'
+        elif var_type == 'polynomial_control_rate2':
+            control_name = var[:-6]
+            rate_path = f'control_rates:{control_name}_rate2'
+        elif var_type == 'parameter':
+            rate_path = f'parameter_vals:{var}'
+        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.isin(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'control_values:{var}'
+            src_units = phase.control_options[var]['units']
+            src_shape = phase.control_options[var]['shape']
+        elif var_type == 'polynomial_control_rate':
+            control_name = var[:-5]
+            path = f'control_rates:{control_name}_rate'
+            control = phase.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'control_rates:{control_name}_rate2'
+            control = phase.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)
+                else:
+                    src_idxs_raw = np.zeros(self.grid_data.subset_num_nodes['all'], dtype=int)
+                    src_idxs = get_src_indices_by_row(src_idxs_raw, options['shape'])
+                    if options['shape'] == (1,):
+                        src_idxs = src_idxs.ravel()
+
+                connection_info.append((f'ode_all.{tgt}', (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
+
+    def _phase_set_state_val(self, phase, name, vals, time_vals, interpolation_kind):
+        """
+        Method to interpolate the provided input and return the variables that need to be set
+        along with their appropriate value.
+
+        Parameters
+        ----------
+        phase : dymos.Phase
+            The phase to which this transcription applies.
+        name : str
+            The name of the phase variable to be set.
+        vals : ndarray or Sequence or float
+            Array of control/state/parameter values.
+        time_vals : ndarray or Sequence or None
+            Array of integration variable values.
+        interpolation_kind : str
+            Specifies the kind of interpolation, as per the scipy.interpolate package.
+            One of ('linear', 'nearest', 'zero', 'slinear', 'quadratic', 'cubic'
+            where 'zero', 'slinear', 'quadratic' and 'cubic' refer to a spline
+            interpolation of zeroth, first, second or third order) or as an
+            integer specifying the order of the spline interpolator to use.
+            Default is 'linear'.
+
+        Returns
+        -------
+        input_data : dict
+            Dict containing the values that need to be set in the phase
+
+        """
+        seg_end_idxs = self.grid_data.subset_node_indices['segment_ends']
+        seg_initial_idxs = seg_end_idxs[::2]
+        seg_final_idxs = seg_end_idxs[1::2]
+
+        input_data = {}
+        if np.isscalar(vals):
+            input_data[f'states:{name}'] = vals
+            input_data[f'initial_states:{name}'] = vals
+            input_data[f'final_states:{name}'] = vals
+
+            if phase.state_options[name]['solve_segments'] == 'forward':
+                input_data[f'picard_update_comp.seg_initial_states:{name}'] = vals
+            else:
+                input_data[f'picard_update_comp.seg_final_states:{name}'] = vals
+        else:
+            interp_vals = phase.interp(name, vals, time_vals,
+                                       nodes='all',
+                                       kind=interpolation_kind)
+            input_data[f'states:{name}'] = interp_vals
+            input_data[f'initial_states:{name}'] = interp_vals[0, ...]
+            input_data[f'final_states:{name}'] = interp_vals[-1, ...]
+
+            if phase.state_options[name]['solve_segments'] == 'forward':
+                input_data[f'picard_update_comp.seg_initial_states:{name}'] = interp_vals[seg_initial_idxs, ...]
+            else:
+                input_data[f'picard_update_comp.seg_final_states:{name}'] = interp_vals[seg_final_idxs, ...]
+
+        return input_data
diff --git a/dymos/transcriptions/picard_shooting/states_comp.py b/dymos/transcriptions/picard_shooting/states_comp.py
new file mode 100644
index 000000000..de7b5583a
--- /dev/null
+++ b/dymos/transcriptions/picard_shooting/states_comp.py
@@ -0,0 +1,109 @@
+"""Definition of the States passthru component."""
+import numpy as np
+
+from openmdao.core.explicitcomponent import ExplicitComponent
+from dymos.transcriptions.grid_data import GridData
+from dymos._options import options as dymos_options
+
+
+class StatesComp(ExplicitComponent):
+    """
+    This component provides inputs for states and outputs their values during Picard iteration.
+
+    This component serves to accept variables 'states:{state_name}' and echo them out as
+    'state_val:{state_name}'.  This is necessary when using NonlinearBlockGS to converge
+    the states.
+
+    Parameters
+    ----------
+    **kwargs : dict
+        Arguments to be passed to the component initialization method.
+
+    Attributes
+    ----------
+    _vars : dict
+        Container mapping name of variables to be muxed with additional data.
+    _input_names : dict
+        Container mapping name of variables to be muxed with associated inputs.
+    """
+
+    def __init__(self, **kwargs):
+        """
+        Instantiate MuxComp and populate private members.
+        """
+        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')
+
+    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['all']
+        state_options = self.options['state_options']
+
+        self.var_names = var_names = {}
+        for state_name, options in state_options.items():
+            var_names[state_name] = {
+                'current_state': f'state_val:{state_name}',
+                'next_state': f'states:{state_name}'
+            }
+
+        ar = np.arange(num_nodes, dtype=int)
+        for state_name, options in state_options.items():
+            shape = options['shape']
+            units = options['units']
+
+            var_names = self.var_names[state_name]
+
+            self.add_input(
+                name=var_names['next_state'],
+                shape=(num_nodes,) + shape,
+                desc=f'Value of state {state_name} available as an input.',
+                units=units)
+
+            self.add_output(
+                name=var_names['current_state'],
+                shape=(num_nodes,) + shape,
+                desc=f'Value of the state {state_name} to be used in the ODE evaluation.',
+                units=units
+            )
+
+            self.declare_partials(of=var_names['current_state'],
+                                  wrt=var_names['next_state'],
+                                  rows=ar, cols=ar, val=1.0)
+
+    def compute(self, inputs, outputs, discrete_inputs=None, discrete_outputs=None):
+        """
+        Mux the inputs into the appropriate outputs.
+
+        Parameters
+        ----------
+        inputs : Vector
+            Unscaled, dimensional input variables read via inputs[key].
+        outputs : Vector
+            Unscaled, dimensional output variables read via outputs[key].
+        discrete_inputs : Vector
+            Discrete input variables read via discrete_inputs[key].
+        discrete_outputs : Vector
+            Discrete output variables read via discrete_outputs[key].
+        """
+        outputs.set_val(inputs.asarray())
diff --git a/dymos/transcriptions/pseudospectral/birkhoff.py b/dymos/transcriptions/pseudospectral/birkhoff.py
index 5617028c7..a87250b35 100644
--- a/dymos/transcriptions/pseudospectral/birkhoff.py
+++ b/dymos/transcriptions/pseudospectral/birkhoff.py
@@ -7,7 +7,7 @@
 from .components import BirkhoffIterGroup, BirkhoffBoundaryGroup
 
 from ..grid_data import BirkhoffGrid
-from dymos.utils.misc import get_rate_units
+from dymos.utils.misc import get_rate_units, _format_phase_constraint_alias
 from dymos.utils.introspection import get_promoted_vars, get_source_metadata
 from dymos.utils.indexing import get_constraint_flat_idxs, get_src_indices_by_row
 
@@ -273,28 +273,6 @@ def configure_defects(self, phase):
             if rate_source_type not in ('state', 'ode'):
                 phase.connect(rate_src_path, f'f_computed:{name}', src_indices=src_idxs)
 
-    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.
@@ -398,7 +376,6 @@ def _get_constraint_kwargs(self, constraint_type, options, phase):
         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 = options['name']
@@ -437,15 +414,13 @@ def _get_constraint_kwargs(self, constraint_type, options, phase):
 
                 constraint_kwargs['indices'] = path_idxs
 
-        alias_map = {'path': 'path_constraint',
-                     'initial': 'initial_boundary_constraint',
-                     'final': 'final_boundary_constraint'}
-
-        str_idxs = '' if options['indices'] is None else f'{options["indices"]}'
+        con_path = constraint_kwargs.pop('constraint_path')
+        con_name = constraint_kwargs.pop('constraint_name')
 
-        constraint_kwargs['alias'] = f'{phase.pathname}->{alias_map[constraint_type]}->{con_name}{str_idxs}'
+        constraint_kwargs['alias'] = _format_phase_constraint_alias(phase, con_name,
+                                                                    constraint_type,
+                                                                    options['indices'])
         constraint_kwargs.pop('name')
-        con_path = constraint_kwargs.pop('constraint_path')
         constraint_kwargs.pop('shape')
         constraint_kwargs['flat_indices'] = True
 
diff --git a/dymos/transcriptions/pseudospectral/pseudospectral_base.py b/dymos/transcriptions/pseudospectral/pseudospectral_base.py
index 2eb94630a..90b1ace5f 100644
--- a/dymos/transcriptions/pseudospectral/pseudospectral_base.py
+++ b/dymos/transcriptions/pseudospectral/pseudospectral_base.py
@@ -6,7 +6,7 @@
 from .components import StateIndependentsComp, StateInterpComp, CollocationComp
 from ..common.timeseries_output_comp import TimeseriesOutputComp
 from ...utils.misc import CoerceDesvar, get_rate_units, reshape_val
-from ...utils.introspection import get_promoted_vars, get_source_metadata, configure_duration_balance_introspection
+from ...utils.introspection import get_promoted_vars, get_source_metadata
 from ...utils.constants import INF_BOUND
 from ...utils.indexing import get_src_indices_by_row
 
@@ -438,55 +438,6 @@ def configure_defects(self, phase):
                               f'continuity_comp.control_rates:{name}_rate2',
                               src_indices=src_idxs, flat_src_indices=True)
 
-    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.
-        """
-
-        if self._implicit_duration:
-            duration_balance_comp = om.BalanceComp()
-            phase.add_subsystem('t_duration_balance_comp', duration_balance_comp)
-
-    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.
-        """
-
-        if self._implicit_duration:
-            duration_balance_comp = phase._get_subsystem('t_duration_balance_comp')
-            configure_duration_balance_introspection(phase)
-            options = phase.time_options['t_duration_balance_options']
-            lower, upper = phase.time_options['duration_bounds']
-            duration_val = phase.time_options['duration_val']
-
-            if options['mult_val'] is None:
-                use_mult = False
-                mult_val = 1.0
-            else:
-                use_mult = True
-                mult_val = options['mult_val']
-
-            src_idx = [-1] if options['index'] is None else [[-1]]+options['index']
-
-            duration_balance_comp.add_balance('t_duration', val=duration_val, lower=lower, upper=upper,
-                                              eq_units=options['units'], units=phase.time_options['units'],
-                                              rhs_val=options['val'], use_mult=use_mult, mult_val=mult_val)
-
-            phase.connect('t_duration_balance_comp.t_duration', 't_duration')
-
-            phase.connect(options['var_path'], 't_duration_balance_comp.lhs:t_duration',
-                          src_indices=src_idx)
-
     def setup_solvers(self, phase):
         """
         Setup the solvers.
diff --git a/dymos/transcriptions/transcription_base.py b/dymos/transcriptions/transcription_base.py
index 6dc1a2017..672751b00 100644
--- a/dymos/transcriptions/transcription_base.py
+++ b/dymos/transcriptions/transcription_base.py
@@ -8,7 +8,8 @@
 from ..utils.constants import INF_BOUND
 from ..utils.indexing import get_constraint_flat_idxs
 from ..utils.introspection import configure_states_introspection, get_promoted_vars, \
-    configure_states_discovery
+    configure_states_discovery, _configure_boundary_balance_introspection
+from ..utils.misc import _unspecified, _format_phase_constraint_alias
 
 
 class TranscriptionBase(object):
@@ -22,7 +23,7 @@ class TranscriptionBase(object):
     """
     def __init__(self, **kwargs):
 
-        self._implicit_duration = False
+        self._implicit_params = False
         self.grid_data = None
 
         self.options = om.OptionsDictionary()
@@ -79,8 +80,8 @@ def setup_time(self, phase):
         # Warn about invalid options
         phase.check_time_options()
 
-        if phase.time_options['t_duration_balance_options']:
-            self._implicit_duration = True
+        if phase.boundary_balance_options:
+            self._implicit_params = True
 
         phase.add_subsystem('param_comp', subsys=ParameterComp(time_options=time_options),
                             promotes_inputs=['*'], promotes_outputs=['*'])
@@ -305,30 +306,75 @@ def setup_ode(self, phase):
         """
         raise NotImplementedError(f'Transcription {self.__class__.__name__} does not implement method setup_ode.')
 
-    def setup_duration_balance(self, phase):
+    def setup_boundary_balance(self, phase):
         """
-        Setup the implicit computation of the phase duration.
+        Setup the implicit computation of the phase boundary balance.
 
         Parameters
         ----------
         phase : dymos.Phase
             The phase object to which this transcription instance applies.
         """
+        if phase.boundary_balance_options:
+            boundary_balance_comp = om.BalanceComp()
+            phase.add_subsystem('boundary_balance_comp', boundary_balance_comp, promotes_outputs=['*'])
 
-        raise NotImplementedError(f'Transcription {self.__class__.__name__} does not implement'
-                                  f' method setup_duration_balance.')
-
-    def configure_duration_balance(self, phase):
+    def configure_boundary_balance(self, phase):
         """
-        Configure the implicit computation of the phase duration.
+        Configure the implicit computation of the phase boundary balance.
 
         Parameters
         ----------
         phase : dymos.Phase
             The phase object to which this transcription instance applies.
         """
-        raise NotImplementedError(f'Transcription {self.__class__.__name__} does not implement'
-                                  f' method setup_duration_balance.')
+        param_balance_comp = phase._get_subsystem('boundary_balance_comp')
+
+        _configure_boundary_balance_introspection(phase)
+
+        for param, options in phase.boundary_balance_options.items():
+            name = options['name']
+            tgt_val = options['tgt_val']
+            loc = options['loc']
+            index = [options['index']] if np.isscalar(options['index']) else options['index']
+
+            # Get the indices to connect based on loc.
+            if loc == 'final':
+                src_idxs = om.slicer[(-1, *index)]
+            elif loc == 'initial':
+                src_idxs = om.slicer[(0, *index)]
+            else:
+                raise ValueError(f'{phase.msginfo}: Value of `loc` for boundary balance `{param}` '
+                                 'must be one of `initial` or `final`, but got `{loc}` instead.')
+
+            # Create the arguments for the balance comp.
+            bal_kwargs = {key: options for key, options in options.items()}
+            try:
+                output_name = bal_kwargs.pop('output_name')
+            except KeyError:
+                output_name = name.split('.')[-1]
+            bal_kwargs.pop('param')
+            bal_kwargs.pop('name')
+            bal_kwargs.pop('tgt_val')
+            bal_kwargs.pop('loc')
+            bal_kwargs.pop('index')
+            bal_kwargs['rhs_val'] = tgt_val
+            bal_kwargs['lhs_name'] = output_name
+
+            prom_param_name = f'parameters:{param}' if param in phase.parameter_options else param
+
+            # Now configure the balance.
+            param_balance_comp.add_balance(name=prom_param_name, **bal_kwargs)
+
+            var_type = phase.classify_var(name)
+            if var_type == 'ode':
+                if name not in phase._timeseries['timeseries']['outputs']:
+                    phase.add_timeseries_output(name, output_name=output_name, units=bal_kwargs.get('eq_units', _unspecified))
+            if var_type == 'state' and name.startswith('initial_states:') or name.startswith('final_states:'):
+                phase.promotes('boundary_balance_comp', inputs=[output_name])
+            else:
+                phase.connect(f'timeseries.{output_name}', f'boundary_balance_comp.{output_name}',
+                              src_indices=src_idxs)
 
     def setup_solvers(self, phase):
         """
@@ -356,7 +402,7 @@ def configure_solvers(self, phase, requires_solvers=None):
         if not phase.options['auto_solvers']:
             return
 
-        req_solvers = {'implicit duration': self._implicit_duration}
+        req_solvers = {'implicit parameters': self._implicit_params}
 
         if requires_solvers is not None:
             req_solvers.update(requires_solvers)
@@ -462,7 +508,6 @@ def _get_constraint_kwargs(self, constraint_type, options, phase):
         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 = options['name']
@@ -509,15 +554,13 @@ def _get_constraint_kwargs(self, constraint_type, options, phase):
 
                 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"]}'
+        con_path = constraint_kwargs.pop('constraint_path')
+        con_name = constraint_kwargs.pop('constraint_name')
 
-        constraint_kwargs['alias'] = f'{phase.pathname}->{alias_map[constraint_type]}->{con_name}{str_idxs}'
+        constraint_kwargs['alias'] = _format_phase_constraint_alias(phase, con_name,
+                                                                    constraint_type,
+                                                                    options['indices'])
         constraint_kwargs.pop('name')
-        con_path = constraint_kwargs.pop('constraint_path')
         constraint_kwargs.pop('shape')
         constraint_kwargs['flat_indices'] = True
 
@@ -628,6 +671,9 @@ def _get_ode(self, phase):
             The OpenMDAO system which serves as the ODE for the given Phase.
 
         """
+        ode = phase._get_subsystem(self._rhs_source)
+        if ode is None:
+            raise AttributeError(f'Phase ODE subsystem {self._rhs_source} not found.')
         return phase._get_subsystem(self._rhs_source)
 
     def get_parameter_connections(self, name, phase):
diff --git a/dymos/utils/introspection.py b/dymos/utils/introspection.py
index 45d9de958..7006bc3eb 100644
--- a/dymos/utils/introspection.py
+++ b/dymos/utils/introspection.py
@@ -6,7 +6,6 @@
 from openmdao.utils.units import simplify_unit
 from openmdao.utils.general_utils import ensure_compatible
 from dymos.utils.misc import _unspecified, _none_or_unspecified
-from .._options import options as dymos_options
 from ..phase.options import StateOptionsDictionary, TimeseriesOutputOptionsDictionary
 from .misc import get_rate_units
 
@@ -40,7 +39,8 @@ def classify_var(var, time_options, state_options, parameter_options, control_op
     -------
     str
         The classification of the given variable, which is one of
-        't', 't_phase', 'state', 'input_control', 'indep_control', 'control_rate',
+        't', 't_phase', 'state',
+        'input_control', 'indep_control', 'control_rate',
         'control_rate2', 'parameter', or 'ode'.
     """
     time_name = time_options['name']
@@ -52,6 +52,10 @@ def classify_var(var, time_options, state_options, parameter_options, control_op
         om.issue_warning('time_phase is deprecated. Please use `t_phase` to obtain the change in the integration '
                          'variable within the current phase.', category=om.OMDeprecationWarning)
         return 't_phase'
+    elif var.startswith('initial_states:'):
+        return 'state'
+    elif var.startswith('final_states:'):
+        return 'state'
     elif var in state_options:
         return 'state'
     elif var in control_options:
@@ -836,7 +840,7 @@ def filter_outputs(patterns, sys):
     return results
 
 
-def configure_duration_balance_introspection(phase):
+def _configure_boundary_balance_introspection(phase):
     """
     Modify duration balance options in-place using introspection of the phase and its ODE.
 
@@ -847,81 +851,189 @@ def configure_duration_balance_introspection(phase):
     """
     time_units = phase.time_options['units']
 
-    options = phase.time_options['t_duration_balance_options']
-
-    # Determine the path to the variable which we will be constraining
-    var = options['name']
-    var_type = phase.classify_var(var)
-
-    if var != options['balance_name'] is not None and var_type != 'ode':
-        om.issue_warning(f"Option 'balance_name' on duration residual {var} is only "
-                         f"valid for ODE outputs. The option is being ignored.", om.UnusedOptionWarning)
-
-    if var_type == 't' or var_type == 't_phase':
-        raise ValueError(f'Cannot use a variable of type {var_type} in the duration balance')
-
-    elif var_type == 'state':
-        prefix = 'states:' if dymos_options['use_timeseries_prefix'] else ''
-        state_shape = phase.state_options[var]['shape']
-        state_units = phase.state_options[var]['units']
-        options['shape'] = state_shape
-        options['units'] = state_units if options['units'] is None else options['units']
-        options['var_path'] = f'timeseries.{prefix}{var}'
-
-    elif var_type == 'parameter':
-        param_shape = phase.parameter_options[var]['shape']
-        param_units = phase.parameter_options[var]['units']
-        options['shape'] = param_shape
-        options['units'] = param_units if options['units'] is None else options['units']
-        options['var_path'] = f'parameter_vals:{var}'
-
-    elif var_type in ['indep_control', 'input_control']:
-        prefix = 'controls:' if dymos_options['use_timeseries_prefix'] else ''
-        control_shape = phase.control_options[var]['shape']
-        control_units = phase.control_options[var]['units']
-
-        options['shape'] = control_shape
-        options['units'] = control_units if options['units'] is None else options['units']
-        options['var_path'] = f'timeseries.{prefix}{var}'
-
-    elif var_type == 'control_rate':
-        prefix = 'control_rates:' if dymos_options['use_timeseries_prefix'] else ''
-        control_name = var[:-5]
-        control_shape = phase.control_options[control_name]['shape']
-        control_units = phase.control_options[control_name]['units']
-        options['shape'] = control_shape
-        options['units'] = get_rate_units(control_units, time_units, deriv=1) \
-            if options['units'] is None else options['units']
-        options['var_path'] = f'timeseries.{prefix}{var}'
-
-    elif var_type == 'control_rate2':
-        prefix = 'control_rates:' if dymos_options['use_timeseries_prefix'] else ''
-        control_name = var[:-6]
-        control_shape = phase.control_options[control_name]['shape']
-        control_units = phase.control_options[control_name]['units']
-        options['shape'] = control_shape
-        options['units'] = get_rate_units(control_units, time_units, deriv=2) \
-            if options['units'] is None else options['units']
-        options['var_path'] = f'timeseries.{prefix}{var}'
-
-    elif var_type == 'timeseries_exec_comp_output':
-        options['shape'] = (1,)
-        options['units'] = None
-        options['var_path'] = f'timeseries.timeseries_exec_comp.{var}'
+    for param_name, options in phase.boundary_balance_options.items():
+        resid_name = options['name']
+        resid_units = options.get('eq_units', _unspecified)
+        resid_shape = options.get('shape', None)
+
+        # Determine the path to the variable which we will be constraining
+        resid_type = phase.classify_var(resid_name)
+
+        if param_name == 't_initial':
+            param_units = phase.time_options['units']
+            param_bounds = phase.time_options['initial_bounds']
+        elif param_name == 't_duration':
+            param_units = phase.time_options['units']
+            param_bounds = phase.time_options['duration_bounds']
+        elif param_name in phase.parameter_options:
+            param_units = phase.parameter_options[param_name]['units']
+            param_bounds = tuple(phase.parameter_options[param_name][k] for k in ('lower', 'upper'))
+        elif param_name.startswith('initial_states:'):
+            state_name = ':'.join(param_name.split(':')[1:])
+            param_units = phase.state_options[state_name]['units']
+            param_bounds = phase.state_options[state_name]['initial_bounds']
+        elif param_name.startswith('final_states:'):
+            state_name = ':'.join(param_name.split(':')[1:])
+            param_units = phase.state_options[state_name]['units']
+            param_bounds = phase.state_options[state_name]['final_bounds']
+        else:
+            raise ValueError(f'{phase.msginfo}: For boundary balance, param must be one of t_initial, t_duration, '
+                             'a parameter in the phase, initial_states:{name}, or final_states:{name}')
 
-    else:
-        # Failed to find variable, assume it is in the ODE. This requires introspection.
-        ode = phase.options['transcription']._get_ode(phase)
+        if options.get('units', _unspecified) is _unspecified:
+            options['units'] = param_units
+
+        if options.get('lower', _unspecified) is _unspecified:
+            if param_bounds is None:
+                options['lower'] = None
+            else:
+                options['lower'] = param_bounds[0]
+
+        if options.get('upper', _unspecified) is _unspecified:
+            if param_bounds is None:
+                options['upper'] = None
+            else:
+                options['upper'] = param_bounds[1]
+
+        if resid_units is _unspecified:
+            if resid_type in ['t', 't_phase']:
+                options['eq_units'] = time_units
+            elif resid_type == 'state':
+                if resid_name.startswith('initial_states:') or resid_name.startswith('final_states:'):
+                    state_name = ':'.join(resid_name.split(':')[1:])
+                else:
+                    state_name = resid_name
+                options['eq_units'] = phase.state_options[state_name]['units']
+            elif resid_type == 'parameter':
+                options['eq_units'] = phase.parameter[resid_name]['units']
+            elif resid_type in ['input_control', 'indep_control']:
+                options['eq_units'] = phase.control_options[resid_name]['units']
+            elif resid_type == 'control_rate':
+                control_units = phase.state_options[resid_name]['units']
+                options['eq_units'] = get_rate_units(control_units, time_units, deriv=1)
+            elif resid_type == 'control_rate2':
+                options['eq_units'] = get_rate_units(control_units, time_units, deriv=2)
+            elif resid_type == 'ode':
+                # Failed to find variable, assume it is in the ODE. This requires introspection.
+                ode = phase.options['transcription']._get_ode(phase)
+                meta = get_source_metadata(ode, src=resid_name, user_units=resid_units, user_shape=resid_shape)
+                options['eq_units'] = meta['units']
+            else:
+                raise ValueError(f'{phase.msginfo}: Unable to find boundary balance name {resid_name}')
+
+
+def _configure_constraint_introspection(phase):
+    """
+    Modify constraint options in-place using introspection of the phase and its ODE.
+
+    Parameters
+    ----------
+    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)]:
+
+        time_units = phase.time_options['units']
+        time_name = phase.time_options['name']
+
+        for con in constraints:
+            # Determine the path to the variable which we will be constraining
+            var = con['name']
+            var_type = phase.classify_var(var)
 
-        meta = get_source_metadata(ode, src=var, user_units=options['units'])
+            if var != con['constraint_name'] is not None and var_type != 'ode':
+                om.issue_warning(f"Option 'constraint_name' on {constraint_type} constraint {var} is only "
+                                 f"valid for ODE outputs. The option is being ignored.", om.UnusedOptionWarning)
 
-        options['shape'] = meta['shape']
-        options['units'] = meta['units']
-        options['var_path'] = f'timeseries.{options["balance_name"]}'
+            if var_type == 't':
+                con['shape'] = (1,)
+                con['units'] = time_units if con['units'] is None else con['units']
+                con['constraint_path'] = f'timeseries.{time_name}'
 
-    if options['shape'] != (1,) and options['index'] is None:
-        raise ValueError(f'Specified variable for the duration balance has shape {options["shape"]} and'
-                         f' has no index specified. The balance may only have shape (1,) or a single index')
+            elif var_type == 't_phase':
+                con['shape'] = (1,)
+                con['units'] = time_units if con['units'] is None else con['units']
+                con['constraint_path'] = f'timeseries.{time_name}_phase'
+
+            elif var_type == 'state':
+                prefix = 'states:' if phase.timeseries_options['use_prefix'] else ''
+                state_shape = phase.state_options[var]['shape']
+                state_units = phase.state_options[var]['units']
+                con['shape'] = state_shape
+                con['units'] = state_units if con['units'] is None else con['units']
+                if birkhoff and constraint_type in ('initial', 'final'):
+                    con['constraint_path'] = f'boundary_vals.{var}'
+                else:
+                    con['constraint_path'] = f'timeseries.{prefix}{var}'
+
+            elif var_type == 'parameter':
+                param_shape = phase.parameter_options[var]['shape']
+                param_units = phase.parameter_options[var]['units']
+                con['shape'] = param_shape
+                con['units'] = param_units if con['units'] is None else con['units']
+                con['constraint_path'] = f'parameter_vals:{var}'
+
+            elif var_type in ['indep_control', 'input_control']:
+                prefix = 'controls:' if phase.timeseries_options['use_prefix'] else ''
+                control_shape = phase.control_options[var]['shape']
+                control_units = phase.control_options[var]['units']
+
+                con['shape'] = control_shape
+                con['units'] = control_units if con['units'] is None else con['units']
+                if birkhoff and constraint_type in ('initial', 'final'):
+                    con['constraint_path'] = f'boundary_vals.{var}'
+                else:
+                    con['constraint_path'] = f'timeseries.{prefix}{var}'
+
+            elif var_type == 'control_rate':
+                prefix = 'control_rates:' if phase.timeseries_options['use_prefix'] else ''
+                control_name = var[:-5]
+                control_shape = phase.control_options[control_name]['shape']
+                control_units = phase.control_options[control_name]['units']
+                con['shape'] = control_shape
+                con['units'] = get_rate_units(control_units, time_units, deriv=1) \
+                    if con['units'] is None else con['units']
+                if birkhoff and constraint_type in ('initial', 'final'):
+                    con['constraint_path'] = f'boundary_vals.{var}'
+                else:
+                    con['constraint_path'] = f'timeseries.{prefix}{var}'
+
+            elif var_type == 'control_rate2':
+                prefix = 'control_rates:' if phase.timeseries_options['use_prefix'] else ''
+                control_name = var[:-6]
+                control_shape = phase.control_options[control_name]['shape']
+                control_units = phase.control_options[control_name]['units']
+                con['shape'] = control_shape
+                con['units'] = get_rate_units(control_units, time_units, deriv=2) \
+                    if con['units'] is None else con['units']
+                if birkhoff and constraint_type in ('initial', 'final'):
+                    con['constraint_path'] = f'boundary_vals.{var}'
+                else:
+                    con['constraint_path'] = f'timeseries.{prefix}{var}'
+
+            elif var_type == 'timeseries_exec_comp_output':
+                con['shape'] = (1,)
+                con['units'] = None
+                con['constraint_path'] = f'timeseries.timeseries_exec_comp.{var}'
+
+            else:
+                # Failed to find variable, assume it is in the ODE. This requires introspection.
+                ode = phase.options['transcription']._get_ode(phase)
+
+                meta = get_source_metadata(ode, src=var, user_units=con['units'], user_shape=con['shape'])
+
+                con['shape'] = meta['shape']
+                con['units'] = meta['units']
+
+                if birkhoff and constraint_type in ('initial', 'final'):
+                    con['constraint_path'] = f'boundary_vals.{var}'
+                else:
+                    con['constraint_path'] = f'timeseries.{con["constraint_name"]}'
 
 
 def _get_targets_metadata(ode, name, user_targets=_unspecified):
diff --git a/dymos/utils/misc.py b/dymos/utils/misc.py
index c0e81a1f4..c0e115063 100644
--- a/dymos/utils/misc.py
+++ b/dymos/utils/misc.py
@@ -78,6 +78,36 @@ def reshape_val(val, shape, num_input_nodes):
     return shaped_val
 
 
+def _format_phase_constraint_alias(phase, con_name, con_type, indices=None):
+    """
+    Get an alias for a constraint of the given type on the given path in the given phase.
+
+    Parameters
+    ----------
+    phase : Phase
+        The dymos phase to which the constraint belongs.
+    con_name : str
+        The name or path of the constraint.
+    con_type : str
+        One of 'initial', 'final', or 'path'.
+    indices : tuple or None
+        The indices of the constraint variable to be constrained. These indices
+        ignore the num_nodes dimension in path constraints.
+
+    Returns
+    -------
+    str
+        The alias of the constraint.
+    """
+    str_idxs = '' if indices is None else f'[{indices}]'
+    if f'.phases.{phase.name}' in phase.pathname:
+        phase_path = phase.pathname.replace(f'.phases.{phase.name}', f'.{phase.name}')
+    else:
+        phase_path = phase.pathname
+    # return f'{phase_path}->{con_type}->{con_name}{str_idxs}'
+    return f'{phase_path}.{con_name}[{con_type}]{str_idxs}'
+
+
 class CoerceDesvar(object):
     """
     Check the desvar options for the appropriate shape and resize accordingly with options.
diff --git a/dymos/utils/ode_utils.py b/dymos/utils/ode_utils.py
index e0dd94baf..1d9ebf5af 100644
--- a/dymos/utils/ode_utils.py
+++ b/dymos/utils/ode_utils.py
@@ -255,6 +255,8 @@ def _make_ode_system(ode_class, num_nodes, ode_init_kwargs=None, calc_exprs=None
     calc_exprs : dict
         A dictionary keyed by exec comp expresions whose associated values
         are the metadata associated with the expression.
+    parameter_options : dict
+        The parameter options of the owning phase.
 
     Returns
     -------
diff --git a/dymos/visualization/timeseries/bokeh_timeseries_report.py b/dymos/visualization/timeseries/bokeh_timeseries_report.py
index 9c7999412..ee818efee 100644
--- a/dymos/visualization/timeseries/bokeh_timeseries_report.py
+++ b/dymos/visualization/timeseries/bokeh_timeseries_report.py
@@ -504,12 +504,14 @@ def make_timeseries_report(prob, solution_record_file=None, simulation_record_fi
                                 sources[f'{var_name}[{str_idxs}]'] = s = f'{var_name}_{str_idxs.replace(",", "_")}'
                                 sol_data_column = sol_data[var_name][(slice(None), *idxs)]
                                 sol_data[s] = sol_data_column
-                                sim_data_column = sim_data[var_name][(slice(None), *idxs)]
-                                sim_data[s] = sim_data_column
+                                if sim_data:
+                                    sim_data_column = sim_data[var_name][(slice(None), *idxs)]
+                                    sim_data[s] = sim_data_column
                         else:
                             sources = {var_name: var_name}
                         sol_source = ColumnDataSource(sol_data)
-                        sim_source = ColumnDataSource(sim_data)
+                        if sim_data:
+                            sim_source = ColumnDataSource(sim_data)
 
                         for var_name_with_idxs, _source in sources.items():
                             legend_items = []


Test patch:
diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_boundary_balance.py b/dymos/examples/brachistochrone/test/test_brachistochrone_boundary_balance.py
new file mode 100644
index 000000000..72eb11ceb
--- /dev/null
+++ b/dymos/examples/brachistochrone/test/test_brachistochrone_boundary_balance.py
@@ -0,0 +1,204 @@
+import os
+import unittest
+from numpy.testing import assert_almost_equal
+
+try:
+    import jax
+    import jax.numpy as jnp
+except ImportError as e:
+    jax = None
+    jnp = None
+
+import openmdao.api as om
+import dymos as dm
+
+from openmdao.utils.testing_utils import use_tempdirs
+from dymos.utils.misc import om_version
+
+
+class GuidedBrachistochroneODE(om.JaxExplicitComponent):
+    """
+    The brachistochrone EOM assuming a linear rate of change in the wire angle wrt time.
+    """
+
+    def initialize(self):
+        self.options.declare('num_nodes', types=int)
+
+    def setup(self):
+        nn = self.options['num_nodes']
+
+        # static parameters
+        self.add_input('v', desc='speed of the bead on the wire', shape=(nn,), units='m/s')
+        self.add_input('g', desc='gravitational acceleration', shape=(1,), units='m/s**2')
+        self.add_input('time_phase', desc='phase elapsed time', shape=(nn,), units='s')
+        self.add_input('theta_rate', desc='guidance parameter', shape=(1,), units='rad/s')
+
+        self.add_output('xdot', desc='velocity component in x', shape=(nn,), units='m/s',
+                        tags=['dymos.state_rate_source:x', 'dymos.state_units:m'])
+
+        self.add_output('ydot', desc='velocity component in y', shape=(nn,), units='m/s',
+                        tags=['dymos.state_rate_source:y', 'dymos.state_units:m'])
+
+        self.add_output('vdot', desc='acceleration magnitude', shape=(nn,), units='m/s**2',
+                        tags=['dymos.state_rate_source:v', 'dymos.state_units:m/s'])
+
+        self.add_output('theta', desc='wire angle', shape=(nn,), val=2, units='rad')
+
+        # TODO: Temporary pending bug fix in JaxExplicitComp
+        self.declare_coloring(method='jax')
+
+    # because our compute primal output depends on static variables, in this case
+    # and self.options['num_nodes'], we must define a get_self_statics method. This method must
+    # return a tuple of all static variables. Their order in the tuple doesn't matter.  If your
+    # component happens to have discrete inputs, do NOT return them here. Discrete inputs are passed
+    # into the compute_primal function individually, after the continuous variables.
+    def get_self_statics(self):
+        # return value must be hashable
+        return self.options['num_nodes'],
+
+    def compute_primal(self, v, g, time_phase, theta_rate):
+        theta = time_phase * theta_rate
+
+        sin_theta = jnp.sin(theta)
+        cos_theta = jnp.cos(theta)
+
+        vdot = g * cos_theta
+        xdot = v * sin_theta
+        ydot = -v * cos_theta
+
+        return xdot, ydot, vdot, theta
+
+
+@unittest.skipIf(jax is None, 'Test Requires Jax')
+@unittest.skipIf(om_version()[0] <= (3, 37, 0), 'Requires OpenMDAO version later than 3.37.0')
+@use_tempdirs
+class TestBrachistochroneBoundaryBalance(unittest.TestCase):
+
+    def make_problem(self, tx):
+        p = om.Problem(model=om.Group())
+
+        traj = dm.Trajectory()
+        phase = dm.Phase(ode_class=GuidedBrachistochroneODE, transcription=tx)
+        p.model.add_subsystem('traj0', traj)
+        traj.add_phase('phase0', phase)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(0.5, 10))
+
+        phase.add_state('x', fix_initial=True, fix_final=False)
+        phase.add_state('y', fix_initial=True, fix_final=False)
+        phase.add_state('v', fix_initial=True, fix_final=False)
+
+        phase.add_parameter('g', units='m/s**2')
+        phase.add_parameter('theta_rate', units='deg/s')
+
+        phase.add_calc_expr('v2 = v * v')
+
+        phase.add_timeseries_output('theta', units='deg')
+        phase.add_timeseries_output('v2', units='m**2/s**2')
+
+        phase.add_boundary_balance(param='t_duration', name='x', tgt_val=10.0, loc='final', lower=0.1, upper=5.0)
+        phase.add_boundary_balance(param='theta_rate', name='y', tgt_val=5.0, loc='final', lower=0.1, upper=100.0, res_ref=1.0)
+
+        phase.nonlinear_solver = om.NewtonSolver(solve_subsystems=True, maxiter=100, iprint=0, atol=1.0E-6, rtol=1.0E-6)
+        phase.nonlinear_solver.linesearch = om.BoundsEnforceLS()
+        phase.linear_solver = om.DirectSolver()
+
+        p.setup()
+
+        phase.set_time_val(initial=0.0, duration=1.8016)
+
+        phase.set_state_val('x', [0, 10])
+        phase.set_state_val('y', [10, 5])
+        phase.set_state_val('v', [0, 9.9])
+
+        phase.set_parameter_val('theta_rate',  20.0)
+        phase.set_parameter_val('g', 9.80665)
+
+        p.final_setup()
+
+        return p
+
+    def run_asserts(self, p):
+
+        t_initial = p.get_val('traj0.phase0.timeseries.time')[0]
+        tf = p.get_val('traj0.phase0.timeseries.time')[-1]
+
+        x0 = p.get_val('traj0.phase0.timeseries.x')[0]
+        xf = p.get_val('traj0.phase0.timeseries.x')[-1]
+
+        y0 = p.get_val('traj0.phase0.timeseries.y')[0]
+        yf = p.get_val('traj0.phase0.timeseries.y')[-1]
+
+        v0 = p.get_val('traj0.phase0.timeseries.v')[0]
+        vf = p.get_val('traj0.phase0.timeseries.v')[-1]
+
+        g = p.get_val('traj0.phase0.parameter_vals:g')[0]
+
+        thetaf = p.get_val('traj0.phase0.timeseries.theta')[-1]
+
+        v2_f = p.get_val('traj0.phase0.timeseries.v2')[-1]
+
+        assert_almost_equal(t_initial, 0.0)
+        assert_almost_equal(x0, 0.0)
+        assert_almost_equal(y0, 10.0)
+        assert_almost_equal(v0, 0.0)
+
+        assert_almost_equal(tf, 1.8016, decimal=4)
+        assert_almost_equal(xf, 10.0, decimal=3)
+        assert_almost_equal(yf, 5.0, decimal=3)
+        assert_almost_equal(vf, 9.902, decimal=3)
+        assert_almost_equal(v2_f, vf**2, decimal=3)
+        assert_almost_equal(g, 9.80665, decimal=3)
+
+        assert_almost_equal(thetaf, 100.12, decimal=0)
+
+        outputs = [op[1]['prom_name'] for op in p.model.list_outputs(out_stream=None, prom_name=True)]
+        self.assertNotIn('traj0.phase0.timeseries.theta_rate', outputs)
+
+    def test_picard(self):
+        for grid_type in 'cgl', 'lgl':
+            with self.subTest(f'{grid_type=}'):
+                tx = dm.PicardShooting(num_segments=1, nodes_per_seg=21, grid_type='cgl')
+                p = self.make_problem(tx=tx)
+
+                p.run_model()
+
+                self.run_asserts(p)
+
+    def test_birkhoff(self):
+        for grid_type in 'cgl', 'lgl':
+            with self.subTest(f'{grid_type=}'):
+                tx = dm.Birkhoff(num_nodes=21, grid_type=grid_type, solve_segments='forward')
+                p = self.make_problem(tx=tx)
+
+                p.run_model()
+
+                self.run_asserts(p)
+
+    def test_radau(self):
+        tx = dm.Radau(num_segments=10, order=3, solve_segments='forward')
+        p = self.make_problem(tx=tx)
+
+        p.run_model()
+
+        self.run_asserts(p)
+
+    def test_gausslobatto(self):
+        tx = dm.GaussLobatto(num_segments=10, order=3, solve_segments='forward')
+        p = self.make_problem(tx=tx)
+
+        p.run_model()
+
+        self.run_asserts(p)
+
+    def test_explicit_shooting(self):
+        tx = dm.ExplicitShooting(num_segments=10, order=3)
+        p = self.make_problem(tx=tx)
+
+        p.run_model()
+
+        self.run_asserts(p)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_callable_ode_class.py b/dymos/examples/brachistochrone/test/test_brachistochrone_callable_ode_class.py
index 0dd0ed9e8..a3043cbbb 100644
--- a/dymos/examples/brachistochrone/test/test_brachistochrone_callable_ode_class.py
+++ b/dymos/examples/brachistochrone/test/test_brachistochrone_callable_ode_class.py
@@ -33,6 +33,9 @@ def _make_problem(self, transcription='gauss-lobatto', num_segments=8, transcrip
         elif transcription == 'birkhoff':
             t = dm.Birkhoff(num_nodes=transcription_order+1,
                             grid_type='cgl')
+        elif transcription.startswith('picard'):
+            grid_type = transcription.split('-')[-1]
+            t = dm.PicardShooting(nodes_per_seg=5, num_segments=8, grid_type=grid_type)
 
         def ode(num_nodes):
             return om.ExecComp(['vdot = g * cos(theta)',
@@ -142,6 +145,10 @@ def test_ex_brachistochrone_birkhoff_uncompressed(self):
         p = self._make_problem(transcription='birkhoff', transcription_order=10)
         self.run_asserts(p)
 
+    def test_ex_brachistochrone_picard_shooting(self):
+        p = self._make_problem(transcription='picard-cgl', run_driver=True)
+        self.run_asserts(p)
+
 
 @use_tempdirs
 class TestInvalidCallableODEClass(unittest.TestCase):
diff --git a/dymos/phase/test/test_implicit_duration.py b/dymos/phase/test/test_implicit_duration.py
index 046b13525..046884b38 100644
--- a/dymos/phase/test/test_implicit_duration.py
+++ b/dymos/phase/test/test_implicit_duration.py
@@ -1,3 +1,4 @@
+from re import M
 import openmdao.api as om
 import dymos as dm
 import numpy as np
@@ -69,6 +70,9 @@ def _make_simple_problem(tx, input_duration=False, fix_duration=False, expr_end=
         else:
             phase.set_duration_balance('pe=9.80665*h', val=0.0, units='m**2/s**2')
 
+        phase.nonlinear_solver = om.NewtonSolver(solve_subsystems=True)
+        phase.linear_solver = om.DirectSolver()
+
         phase.set_simulate_options(rtol=1.0E-9, atol=1.0E-9)
 
         p.driver = om.ScipyOptimizeDriver(optimizer='SLSQP')
@@ -79,10 +83,9 @@ def _make_simple_problem(tx, input_duration=False, fix_duration=False, expr_end=
 
         p.setup()
 
-        p.set_val('traj.phase.t_initial', 0)
-        p.set_val('traj.phase.t_duration', 2)
-        p.set_val('traj.phase.states:h', phase.interp('h', [30, 0]))
-        p.set_val('traj.phase.states:v', phase.interp('v', [0, -10]))
+        phase.set_time_val(initial=0, duration=2)
+        phase.set_state_val('h', [30, 0])
+        phase.set_state_val('v', [0, -10])
 
         return p
 
@@ -101,6 +104,9 @@ def _make_matrix_problem(tx, shape_error=False):
         index = None if shape_error else [[0], [1]]
         phase.set_duration_balance('z', val=0.0, index=index)
 
+        phase.nonlinear_solver = om.NewtonSolver(solve_subsystems=True)
+        phase.linear_solver = om.DirectSolver()
+
         phase.set_simulate_options(rtol=1.0E-9, atol=1.0E-9)
 
         phase.nonlinear_solver = om.NewtonSolver(solve_subsystems=True)
@@ -132,7 +138,7 @@ def test_implicit_duration_gl(self):
 
         p = self._make_simple_problem(tx)
 
-        dm.run_problem(p, run_driver=False, simulate=False)
+        dm.run_problem(p, run_driver=False, simulate=True, make_plots=True)
 
         assert_near_equal(p.get_val('traj.phase.timeseries.time')[-1], 2.4735192, tolerance=1E-6)
         assert_near_equal(p.get_val('traj.phase.timeseries.h')[-1], 0.0, tolerance=1E-6)
@@ -140,43 +146,12 @@ def test_implicit_duration_gl(self):
     def test_implicit_duration_shooting(self):
         tx = dm.ExplicitShooting(num_segments=12, order=3)
 
-        expected = 'Transcription ExplicitShooting does not implement method setup_duration_balance'
-
-        with self.assertRaises(NotImplementedError) as e:
-            self._make_simple_problem(tx)
-
-        self.assertEqual(expected, str(e.exception))
-
-    def test_implicit_duration_input_duration(self):
-        tx = dm.Radau(num_segments=12, order=3, solve_segments=False)
-
-        expected = 'Cannot implicitly solve for phase duration when input_duration is True'
-
-        with self.assertRaises(ValueError) as e:
-            self._make_simple_problem(tx, input_duration=True)
-
-        self.assertEqual(expected, str(e.exception))
-
-    def test_implicit_duration_fix_duration(self):
-        tx = dm.Radau(num_segments=12, order=3, solve_segments=False)
-
-        expected = 'Cannot implicitly solve for phase duration when fix_duration is True'
-
-        with self.assertRaises(ValueError) as e:
-            self._make_simple_problem(tx, fix_duration=True)
-
-        self.assertEqual(expected, str(e.exception))
-
-    def test_implicit_duration_matrix_state_exception(self):
-        tx = dm.Radau(num_segments=12, order=3, solve_segments=False)
-
-        expected = 'Specified variable for the duration balance has shape (2, 2) and' \
-                   ' has no index specified. The balance may only have shape (1,) or a single index'
+        p = self._make_simple_problem(tx)
 
-        with self.assertRaises(ValueError) as e:
-            self._make_matrix_problem(tx, shape_error=True)
+        dm.run_problem(p, run_driver=False, simulate=False)
 
-        self.assertEqual(expected, str(e.exception))
+        assert_near_equal(p.get_val('traj.phase.timeseries.time')[-1], 2.4735192, tolerance=1E-6)
+        assert_near_equal(p.get_val('traj.phase.timeseries.h')[-1], 0.0, tolerance=1E-6)
 
     def test_implicit_duration_matrix_state(self):
         tx = dm.Radau(num_segments=12, order=3, solve_segments=False)
@@ -187,7 +162,7 @@ def test_implicit_duration_matrix_state(self):
 
         assert_near_equal(p.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
 
-    @unittest.skip('ODE Expressions currently cannot be provided in balances.')
+    @unittest.skip("set_duration_balance/add_boundary_balance doesn't currently work with expressions.")
     @require_pyoptsparse(optimizer='IPOPT')
     def test_implicit_duration_radau_expr_condition(self):
         tx = dm.Radau(num_segments=12, order=3, solve_segments=False)
@@ -201,7 +176,7 @@ def test_implicit_duration_radau_expr_condition(self):
 
         print((p.get_val('traj.phase.timeseries.v')[-1]) ** 2 / 2)
 
-    @unittest.skip('ODE Expressions currently cannot be provided in balances.')
+    @unittest.skip("set_duration_balance/add_boundary_balance doesn't currently work with expressions.")
     def test_implicit_duration_gl_expr_condition(self):
         tx = dm.GaussLobatto(num_segments=12, order=3, solve_segments=False)
 
diff --git a/dymos/transcriptions/picard_shooting/test/test_birkhoff_picard_iter_group.py b/dymos/transcriptions/picard_shooting/test/test_birkhoff_picard_iter_group.py
new file mode 100644
index 000000000..65e07a2b9
--- /dev/null
+++ b/dymos/transcriptions/picard_shooting/test/test_birkhoff_picard_iter_group.py
@@ -0,0 +1,301 @@
+import time
+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 openmdao.utils.testing_utils import use_tempdirs
+
+from dymos.utils.misc import GroupWrapperConfig, om_version
+from dymos.utils.testing_utils import PhaseStub, SimpleODE
+from dymos.transcriptions.picard_shooting.birkhoff_picard_iter_group import BirkhoffPicardIterGroup
+from dymos.phase.options import StateOptionsDictionary, TimeOptionsDictionary
+from dymos.transcriptions.grid_data import BirkhoffGrid, GaussLobattoGrid
+
+
+BirkhoffPicardIterGroup = GroupWrapperConfig(BirkhoffPicardIterGroup, [PhaseStub()])
+
+
+@use_tempdirs
+@unittest.skipIf(om_version()[0] < (3, 37, 0), 'Requires OpenMDAO version later than 3.37.0')
+class TestBirkhoffPicardIterGroup(unittest.TestCase):
+
+    def test_birkhoff_picard_solve_segments(self):
+        for direction in ['forward', 'backward']:
+            for grid_type in ['lgl', 'cgl']:
+                for nl_solver in ['newton', 'nlbgs']:
+                    with self.subTest(msg=f'{direction=} {grid_type=} {nl_solver=}'):
+                        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'] = direction
+                            state_options['x']['rate_source'] = 'x_dot'
+
+                            time_options = TimeOptionsDictionary()
+                            grid_data = GaussLobattoGrid(nodes_per_seg=11, num_segments=1)
+
+                            nn = grid_data.subset_num_nodes['all']
+                            ode_class = SimpleODE
+
+                            p = om.Problem()
+                            p.model.add_subsystem('birkhoff', BirkhoffPicardIterGroup(state_options=state_options,
+                                                                                      time_units=time_options['units'],
+                                                                                      grid_data=grid_data,
+                                                                                      ode_class=ode_class))
+
+                            birkhoff = p.model._get_subsystem('birkhoff')
+
+                            if nl_solver == 'newton':
+                                birkhoff.nonlinear_solver = om.NewtonSolver(solve_subsystems=True, maxiter=100)
+                            elif nl_solver == 'nlbgs':
+                                birkhoff.nonlinear_solver = om.NonlinearBlockGS(maxiter=201, use_aitken=True)
+                            else:
+                                birkhoff.nonlinear_solver = om.NonlinearRunOnce()
+                            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) * 1.0
+
+                            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.final_setup()
+
+                            if direction == 'forward':
+                                p.set_val('birkhoff.picard_update_comp.seg_initial_states:x', 0.5)
+                            else:
+                                p.set_val('birkhoff.picard_update_comp.seg_final_states:x', solution[-1])
+                            p.set_val('birkhoff.ode_all.t', times)
+                            p.set_val('birkhoff.ode_all.p', 1.0)
+
+                            p.final_setup()
+
+                            t_start = time.perf_counter()
+                            p.run_model()
+                            t_end = time.perf_counter()
+
+                            print(f"Elapsed time: {t_end-t_start:.4f} seconds")
+
+                            assert_near_equal(solution, p.get_val('birkhoff.states:x'), tolerance=1.0E-9)
+                            assert_near_equal(dsolution_dt.ravel(),
+                                              p.get_val('birkhoff.picard_update_comp.f_computed:x').ravel(),
+                                              tolerance=1.0E-9)
+
+                            cpd = p.check_partials(method='cs', compact_print=True,
+                                                   show_only_incorrect=True, out_stream=None)
+                            assert_check_partials(cpd, atol=1.0E-5, rtol=1.0E-5)
+
+    @unittest.skip('This test is a demonstation of the inability of Birkhoff-Picard '
+                   'iteration to solve highly nonlinear systems.')
+    def test_birkhoff_solve_segments_lorenz(self):
+
+        class LorenzAttractorODE(om.JaxExplicitComponent):
+            """
+            The ODE for the Lorenz attractor.
+            """
+
+            def initialize(self):
+                """
+                All Dymos ODE systems are required to have an option "num_nodes",
+                which is the number of points at which the ODE is simultaneously evaluated.
+
+                This will be set by the Phase during setup once the transcription details are known.
+                """
+                self.options.declare('num_nodes', types=(int,))
+
+            def setup(self):
+                """
+                In setup, we add inputs and outputs.
+
+                The first dimension is assumed to pertain to the index of the node.
+
+                An input that's a scalar at each node should have a shape of
+                (num_nodes, 1) or just (num_nodes,).
+
+                For vectors or matrices, it's just the shape of the matrix at each
+                node prepended with num_nodes.
+
+                We provide units for the scalars, but OpenMDAO doesn't do unit conversion on an index-by-index basis,
+                so we just assume that no unit conversion should be done for the S matrix and K vector.
+                """
+                nn = self.options['num_nodes']
+
+                # ODE inputs
+                self.add_input('x', shape=(nn,))
+                self.add_input('y', shape=(nn,))
+                self.add_input('z', shape=(nn,))
+                self.add_input('s', shape=(1,))
+                self.add_input('r', shape=(1,))
+                self.add_input('b', shape=(1,))
+
+                # State rates
+                self.add_output('x_dot', shape=(nn,), units='1/s', tags=['dymos.state_rate_source:x'])
+                self.add_output('y_dot', shape=(nn,), units='1/s', tags=['dymos.state_rate_source:y'])
+                self.add_output('z_dot', shape=(nn,), units='1/s', tags=['dymos.state_rate_source:z'])
+
+            # because our compute primal output depends on static variables, in this case
+            # and self.options['num_nodes'], we must define a get_self_statics method. This method must
+            # return a tuple of all static variables. Their order in the tuple doesn't matter.  If your
+            # component happens to have discrete inputs, do NOT return them here. Discrete inputs are passed
+            # into the compute_primal function individually, after the continuous variables.
+            def get_self_statics(self):
+                # return value must be hashable
+                return self.options['num_nodes'],
+
+            def setup_partials(self):
+                nn = self.options['num_nodes']
+                ar = np.arange(nn, dtype=int)
+                self.declare_partials('x_dot', 'x', rows=ar, cols=ar)
+                self.declare_partials('x_dot', 'y', rows=ar, cols=ar)
+                self.declare_partials('x_dot', 's')
+
+                self.declare_partials('y_dot', 'x', rows=ar, cols=ar)
+                self.declare_partials('y_dot', 'y', rows=ar, cols=ar)
+                self.declare_partials('y_dot', 'z', rows=ar, cols=ar)
+                self.declare_partials('y_dot', 'r')
+
+                self.declare_partials('z_dot', 'x', rows=ar, cols=ar)
+                self.declare_partials('z_dot', 'y', rows=ar, cols=ar)
+                self.declare_partials('z_dot', 'z', rows=ar, cols=ar)
+                self.declare_partials('z_dot', 'b')
+
+            def compute_primal(self, x, y, z, s=10, r=28, b=2.667):
+                """
+                Parameters
+                ----------
+                x : array-like
+                    The x state.
+                y : array-like
+                    The y state.
+                z : array-like
+                    The z state.
+                s : float
+                    The Lorenz s parameter.
+                r : float
+                    The Lorenz r parameter.
+                b : float
+                    The Lorenz b parameter.
+
+                Returns
+                -------
+                x_dot : array-like
+                    The time-derivative of the x state.
+                y_dot : array-like
+                    The time-derivative of the y state.
+                z_dot : array-like
+                    The time-derivative of the z state.
+                """
+                x_dot = s * (y - x)
+                y_dot = r * x - y - x*z
+                z_dot = x * y - b * z
+                return x_dot, y_dot, z_dot
+
+        for direction in ['forward']:
+            for grid_type in ['lgl']:
+                with self.subTest(msg=grid_type):
+                    with dymos.options.temporary(include_check_partials=True):
+
+                        state_options = {'x': StateOptionsDictionary(),
+                                         'y': StateOptionsDictionary(),
+                                         'z': StateOptionsDictionary()}
+
+                        state_options['x']['shape'] = (1,)
+                        state_options['x']['units'] = None
+                        state_options['x']['targets'] = ['x']
+                        state_options['x']['initial_bounds'] = (None, None)
+                        state_options['x']['final_bounds'] = (None, None)
+                        state_options['x']['solve_segments'] = direction
+                        state_options['x']['rate_source'] = 'x_dot'
+
+                        state_options['y']['shape'] = (1,)
+                        state_options['y']['units'] = None
+                        state_options['y']['targets'] = ['y']
+                        state_options['y']['initial_bounds'] = (None, None)
+                        state_options['y']['final_bounds'] = (None, None)
+                        state_options['y']['solve_segments'] = direction
+                        state_options['y']['rate_source'] = 'y_dot'
+
+                        state_options['z']['shape'] = (1,)
+                        state_options['z']['units'] = None
+                        state_options['z']['targets'] = ['z']
+                        state_options['z']['initial_bounds'] = (None, None)
+                        state_options['z']['final_bounds'] = (None, None)
+                        state_options['z']['solve_segments'] = direction
+                        state_options['z']['rate_source'] = 'z_dot'
+
+                        time_options = TimeOptionsDictionary()
+                        grid_data = BirkhoffGrid(num_nodes=300, grid_type=grid_type)
+                        nn = grid_data.subset_num_nodes['all']
+                        ode_class = LorenzAttractorODE
+
+                        p = om.Problem()
+                        p.model.add_subsystem('birkhoff', BirkhoffPicardIterGroup(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, maxiter=100, iprint=0)
+                        birkhoff.nonlinear_solver = om.NonlinearBlockGS(maxiter=500, use_aitken=True, iprint=0)
+                        birkhoff.linear_solver = om.DirectSolver()
+
+                        p.setup(force_alloc_complex=True)
+
+                        # Instead of using the TimeComp just transform the node segment taus onto [0, 100]
+                        times = (grid_data.node_stau + 1) * 0.3
+                        dt_dstau = (times[-1] / 2.0) * np.ones(nn)
+
+                        if direction == 'forward':
+                            p.set_val('birkhoff.initial_states:x', 0.0)
+                            p.set_val('birkhoff.initial_states:y', 1.0)
+                            p.set_val('birkhoff.initial_states:z', 1.05)
+
+                            p.set_val('birkhoff.states:x', 0.0)
+                            p.set_val('birkhoff.states:y', 1.0)
+                            p.set_val('birkhoff.states:z', 1.05)
+                        else:
+                            raise ValueError('not set up for backwards prop yet')
+                            # p.set_val('birkhoff.final_states:x', solution[-1])
+                        p.set_val('birkhoff.dt_dstau', dt_dstau)
+                        p.set_val('birkhoff.ode_all.s', 10.0)
+                        p.set_val('birkhoff.ode_all.r', 28.0)
+                        p.set_val('birkhoff.ode_all.b', 2.667)
+
+                        p.final_setup()
+                        p.run_model()
+
+                    p.model.list_vars()
+
+                    x = p.get_val('birkhoff.states:x')
+                    y = p.get_val('birkhoff.states:y')
+                    z = p.get_val('birkhoff.states:z')
+
+                    import matplotlib.pyplot as plt
+                    ax = plt.figure().add_subplot(projection='3d')
+                    ax.plot(x, y, z, '-')
+                    plt.show()
+
+                    # assert_near_equal(solution, p.get_val('birkhoff.states:x'), tolerance=1.0E-9)
+                    # assert_near_equal(dsolution_dt.ravel(), 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='fd', compact_print=True, out_stream=None)
+                    # assert_check_partials(cpd, atol=1.0E-5, rtol=1.0E-5)
+
+
+if __name__ == '__main__':
+    unittest.main()
diff --git a/dymos/transcriptions/picard_shooting/test/test_multiple_shooting_iter_group.py b/dymos/transcriptions/picard_shooting/test/test_multiple_shooting_iter_group.py
new file mode 100644
index 000000000..95ddfe438
--- /dev/null
+++ b/dymos/transcriptions/picard_shooting/test/test_multiple_shooting_iter_group.py
@@ -0,0 +1,230 @@
+import time
+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 openmdao.utils.testing_utils import use_tempdirs
+
+from dymos.utils.misc import GroupWrapperConfig
+from dymos.utils.testing_utils import PhaseStub, SimpleODE
+from dymos.transcriptions.common.time_comp import TimeComp
+from dymos.transcriptions.picard_shooting.multiple_shooting_iter_group import MultipleShootingIterGroup
+from dymos.phase.options import StateOptionsDictionary, TimeOptionsDictionary
+from dymos.transcriptions.grid_data import GaussLobattoGrid, ChebyshevGaussLobattoGrid
+
+
+TimeComp = GroupWrapperConfig(TimeComp, [])
+MultipleShootingIterGroupWrapped = GroupWrapperConfig(MultipleShootingIterGroup, [PhaseStub()])
+
+
+class LorenzAttractorODE(om.JaxExplicitComponent):
+    """
+    The ODE for the Lorenz attractor.
+    """
+
+    def initialize(self):
+        """
+        All Dymos ODE systems are required to have an option "num_nodes",
+        which is the number of points at which the ODE is simultaneously evaluated.
+
+        This will be set by the Phase during setup once the transcription details are known.
+        """
+        self.options.declare('num_nodes', types=(int,))
+
+    def setup(self):
+        """
+        In setup, we add inputs and outputs.
+
+        The first dimension is assumed to pertain to the index of the node.
+
+        An input that's a scalar at each node should have a shape of
+        (num_nodes, 1) or just (num_nodes,).
+
+        For vectors or matrices, it's just the shape of the matrix at each
+        node prepended with num_nodes.
+
+        We provide units for the scalars, but OpenMDAO doesn't do unit conversion on an index-by-index basis,
+        so we just assume that no unit conversion should be done for the S matrix and K vector.
+        """
+        nn = self.options['num_nodes']
+
+        # ODE inputs
+        self.add_input('x', shape=(nn,))
+        self.add_input('y', shape=(nn,))
+        self.add_input('z', shape=(nn,))
+        self.add_input('s', shape=(1,))
+        self.add_input('r', shape=(1,))
+        self.add_input('b', shape=(1,))
+
+        # State rates
+        self.add_output('x_dot', shape=(nn,), units='1/s', tags=['dymos.state_rate_source:x'])
+        self.add_output('y_dot', shape=(nn,), units='1/s', tags=['dymos.state_rate_source:y'])
+        self.add_output('z_dot', shape=(nn,), units='1/s', tags=['dymos.state_rate_source:z'])
+
+    # because our compute primal output depends on static variables, in this case
+    # and self.options['num_nodes'], we must define a get_self_statics method. This method must
+    # return a tuple of all static variables. Their order in the tuple doesn't matter.  If your
+    # component happens to have discrete inputs, do NOT return them here. Discrete inputs are passed
+    # into the compute_primal function individually, after the continuous variables.
+    def get_self_statics(self):
+        # return value must be hashable
+        return self.options['num_nodes'],
+
+    def setup_partials(self):
+        nn = self.options['num_nodes']
+        ar = np.arange(nn, dtype=int)
+        self.declare_partials('x_dot', 'x', rows=ar, cols=ar)
+        self.declare_partials('x_dot', 'y', rows=ar, cols=ar)
+        self.declare_partials('x_dot', 's')
+
+        self.declare_partials('y_dot', 'x', rows=ar, cols=ar)
+        self.declare_partials('y_dot', 'y', rows=ar, cols=ar)
+        self.declare_partials('y_dot', 'z', rows=ar, cols=ar)
+        self.declare_partials('y_dot', 'r')
+
+        self.declare_partials('z_dot', 'x', rows=ar, cols=ar)
+        self.declare_partials('z_dot', 'y', rows=ar, cols=ar)
+        self.declare_partials('z_dot', 'z', rows=ar, cols=ar)
+        self.declare_partials('z_dot', 'b')
+
+    def compute_primal(self, x, y, z, s=10, r=28, b=2.667):
+        """
+        Parameters
+        ----------
+        x : array-like
+            The x state.
+        y : array-like
+            The y state.
+        z : array-like
+            The z state.
+        s : float
+            The Lorenz s parameter.
+        r : float
+            The Lorenz r parameter.
+        b : float
+            The Lorenz b parameter.
+
+        Returns
+        -------
+        x_dot : array-like
+            The time-derivative of the x state.
+        y_dot : array-like
+            The time-derivative of the y state.
+        z_dot : array-like
+            The time-derivative of the z state.
+        """
+        x_dot = s * (y - x)
+        y_dot = r * x - y - x*z
+        z_dot = x * y - b * z
+        return x_dot, y_dot, z_dot
+
+
+@use_tempdirs
+class TestMultipleShootingIterGroup(unittest.TestCase):
+
+    def test_multiple_shooting_iter_group(self):
+        for direction in ['forward', 'backward']:
+            for grid_type in [GaussLobattoGrid, ChebyshevGaussLobattoGrid]:
+                for num_seg, nodes_per_seg in [(1, 21), (3, 11)]:
+                    grid_data = grid_type(nodes_per_seg=nodes_per_seg, num_segments=num_seg)
+                    for nl_solver in [om.NonlinearBlockGS(use_aitken=True, maxiter=100, iprint=0),
+                                      om.NewtonSolver(solve_subsystems=True, maxiter=100, iprint=0)]:
+                        with self.subTest(msg=f'{direction=} grid={grid_data}'):
+                            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'] = direction
+                                state_options['x']['rate_source'] = 'x_dot'
+
+                                time_options = TimeOptionsDictionary()
+                                ode_class = SimpleODE
+
+                                p = om.Problem()
+
+                                class PhaseStub(om.Group):
+
+                                    def setup(self):
+                                        self.state_options = state_options
+                                        self.time_options = time_options
+
+                                    def classify_var(self, var):
+                                        return 'ode'
+
+                                    def configure(self):
+                                        self._get_subsystem('ms').configure_io(self)
+                                        self._get_subsystem('time').configure_io()
+
+                                p.model = PhaseStub()
+
+                                p.model.add_subsystem('time', TimeComp(num_nodes=grid_data.num_nodes,
+                                                                       node_ptau=grid_data.node_ptau,
+                                                                       node_dptau_dstau=grid_data.node_dptau_dstau,
+                                                                       units='s'))
+
+                                p.model.add_subsystem('ms', MultipleShootingIterGroup(state_options=state_options,
+                                                                                      time_units=time_options['units'],
+                                                                                      grid_data=grid_data,
+                                                                                      ode_class=ode_class,
+                                                                                      ms_nonlinear_solver=nl_solver),
+                                                      promotes=['ode_all*'])
+                                p.model.connect('time.t', 'ode_all.t')
+                                p.model.connect('time.dt_dstau', 'ms.picard_update_comp.dt_dstau')
+
+                                ms = p.model._get_subsystem('ms')
+
+                                ms.nonlinear_solver = nl_solver
+                                ms.linear_solver = om.DirectSolver()
+
+                                p.setup(force_alloc_complex=True)
+
+                                p.set_val('time.t_initial', 0.0)
+                                p.set_val('time.t_duration', 2.0)
+
+                                def solution(t):
+                                    return t**2 + 2 * t + 1 - 0.5 * np.exp(t)
+
+                                def dsolution_dt(t):
+                                    return 2 * t + 2 - 0.5 * np.exp(t)
+
+                                if direction == 'forward':
+                                    p.set_val('ms.seg_initial_states:x', 0.5)
+                                    p.set_val('ms.initial_states:x', 0.5)
+                                else:
+                                    p.set_val('ms.seg_final_states:x', solution(2.0))
+                                    p.set_val('ms.final_states:x', solution(2.0))
+
+                                p.set_val('ode_all.p', 1.0)
+
+                                p.final_setup()
+
+                                t_start = time.perf_counter()
+                                p.run_model()
+                                t_end = time.perf_counter()
+
+                                # print(f"Elapsed time: {t_end-t_start:.4f} seconds")
+
+                                t = p.get_val('time.t')
+                                x = p.get_val('ms.states:x')
+                                x_dot = p.get_val('ode_all.x_dot')
+
+                                assert_near_equal(solution(t), x.ravel(), tolerance=1.0E-9)
+                                assert_near_equal(dsolution_dt(t), x_dot.ravel(), tolerance=1.0E-9)
+                                assert_near_equal(solution(0), p.get_val('ms.initial_states:x').ravel(), tolerance=1.0E-9)
+                                assert_near_equal(solution(t[-1]), p.get_val('ms.final_states:x').ravel(), tolerance=1.0E-9)
+
+                                cpd = p.check_partials(method='fd', compact_print=False, out_stream=None)
+                                assert_check_partials(cpd, atol=1.0E-5, rtol=1.0E-5)
+
+
+if __name__ == '__main__':
+    unittest.main()