Repository: OpenMDAO/dymos

Issue:
Add ability to use non-time-stopping conditions for phases.
### Proposed feature.

In dymos, Phases currently end based on a time (t_duration), regardless of transcription.
In an effort to make dymos more capable of performing analysis without any need to invoke an optimizer, it may be desirable to end a phase on some other condition.
This capability is implemented in solve_ivp, for instance.

For the dymos implementation, the user would specify the name of an output, a target value of that output, and a residual equation would be formed that would serve as the residual function for `t_duration`.

For instance, to terminate a phase when output `alt` is 50 meters, a user might say

```
phase.set_duration_residual('alt', equals=50, units='m')
```

### Which outputs could be used as residuals?
For non-pseudospectral-based phase transcriptions, altitude could be any output that is dependent upon `t_duration`, such as states or ODE outputs.  Controls and parameters would not be valid residuals.

### What transcriptions would support duration residuals?
Duration residuals are not valid for pseudospectral phases **unless** `solve_segments` is used to force the states and ODE outputs to observe the equations of motion without invoking the optimizer. Doing so without using `solve_segments` will likely result in a singular jacobian error.  The Analytic transcription will not support `set_duration_residual` (raising NotImplementedError).  All other transcriptions will allow it but it will be up to the user to ensure that it is posed in a valid way.

### What happens under the hood when using a duration residual?
For the Explicit transcription, we will probably rely on the event functionality in `scipy.integrate.solve_ivp` to terminate the integration when the exit condition is reached.

For Pseudospectral transcriptions `t_duration` will become an implicit output associated with the user's given residual function. A Newton solver at the phase level will converge the residual.

### Example

```python
TBD
```


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

coverage: 92.495% (-0.3%) from 92.815% when pulling **a581f98b868a0a25d34522cb53ef43e0800deedb on kaushikponnapalli:implicit_duration** into **362737e28eb52b13e657c9cea3c7f8c5744d5062 on OpenMDAO:master**.

review-notebook-app[bot]: Check out this pull request on&nbsp; <a href="https://app.reviewnb.com/OpenMDAO/dymos/pull/939"><img align="absmiddle"  alt="ReviewNB" height="28" class="BotMessageButtonImage" src="https://raw.githubusercontent.com/ReviewNB/support/master/images/button_reviewnb.png"/></a> 

 See visual diffs & provide feedback on Jupyter Notebooks. 

---

 <i>Powered by <a href='https://www.reviewnb.com/?utm_source=gh'>ReviewNB</a></i>

Patch:
diff --git a/docs/dymos_book/_toc.yml b/docs/dymos_book/_toc.yml
index 29561f73..9d3147da 100644
--- a/docs/dymos_book/_toc.yml
+++ b/docs/dymos_book/_toc.yml
@@ -41,6 +41,7 @@ parts:
     - file: examples/water_rocket/water_rocket
     - file: examples/cart_pole/cart_pole
     - file: examples/brachistochrone/brachistochrone_tandem_phases
+    - file: examples/cannonball_implicit_duration/cannonball_implicit_duration
 - caption: Feature Reference
   chapters:
   - file: features/phases/phases_list
diff --git a/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb b/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb
new file mode 100644
index 00000000..2eab4113
--- /dev/null
+++ b/docs/dymos_book/examples/cannonball_implicit_duration/cannonball_implicit_duration.ipynb
@@ -0,0 +1,471 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "active-ipynb",
+     "remove-input",
+     "remove-output"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "# This cell is mandatory in all Dymos documentation notebooks.\n",
+    "missing_packages = []\n",
+    "try:\n",
+    "    import openmdao.api as om\n",
+    "except ImportError:\n",
+    "    if 'google.colab' in str(get_ipython()):\n",
+    "        !python -m pip install openmdao[notebooks]\n",
+    "    else:\n",
+    "        missing_packages.append('openmdao')\n",
+    "try:\n",
+    "    import dymos as dm\n",
+    "except ImportError:\n",
+    "    if 'google.colab' in str(get_ipython()):\n",
+    "        !python -m pip install dymos\n",
+    "    else:\n",
+    "        missing_packages.append('dymos')\n",
+    "try:\n",
+    "    import pyoptsparse\n",
+    "except ImportError:\n",
+    "    if 'google.colab' in str(get_ipython()):\n",
+    "        !pip install -q condacolab\n",
+    "        import condacolab\n",
+    "        condacolab.install_miniconda()\n",
+    "        !conda install -c conda-forge pyoptsparse\n",
+    "    else:\n",
+    "        missing_packages.append('pyoptsparse')\n",
+    "if missing_packages:\n",
+    "    raise EnvironmentError('This notebook requires the following packages '\n",
+    "                           'please install them and restart this notebook\\'s runtime: {\",\".join(missing_packages)}')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Cannonball Implicit Duration\n",
+    "\n",
+    "This example demonstrates determining the duration of a phase using a \n",
+    "nonlinear solver to satisfy an endpoint condition. As in the multi-phase\n",
+    "cannonball problem, we will simultaneously find the optimal design for \n",
+    "the cannonball (its radius) and the optimal flight profile (its launch\n",
+    "angle). However, here we will use a single phase that terminates at the\n",
+    "condition\n",
+    "\n",
+    "\\begin{align}\n",
+    "    h(t_f) = 0\n",
+    "\\end{align}\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The cannonball sizing component\n",
+    "\n",
+    "This component computes the area (needed to compute drag) and mass (needed to compute energy) of a cannonball of a given radius and density.\n",
+    "\n",
+    "This component sits upstream of the trajectory model and feeds its outputs to the trajectory as parameters."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "output_scroll"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "from scipy.interpolate import interp1d\n",
+    "import matplotlib.pyplot as plt\n",
+    "\n",
+    "import openmdao.api as om\n",
+    "\n",
+    "import dymos as dm\n",
+    "from dymos.models.atmosphere.atmos_1976 import USatm1976Data\n",
+    "\n",
+    "#############################################\n",
+    "# Component for the design part of the model\n",
+    "#############################################\n",
+    "class CannonballSizeComp(om.ExplicitComponent):\n",
+    "    \"\"\"\n",
+    "    Compute the area and mass of a cannonball with a given radius and density.\n",
+    "\n",
+    "    Notes\n",
+    "    -----\n",
+    "    This component is not vectorized with 'num_nodes' as is the usual way\n",
+    "    with Dymos, but is instead intended to compute a scalar mass and reference\n",
+    "    area from scalar radius and density inputs. This component does not reside\n",
+    "    in the ODE but instead its outputs are connected to the trajectory via\n",
+    "    input design parameters.\n",
+    "    \"\"\"\n",
+    "    def setup(self):\n",
+    "        self.add_input(name='radius', val=1.0, desc='cannonball radius', units='m')\n",
+    "        self.add_input(name='dens', val=7870., desc='cannonball density', units='kg/m**3')\n",
+    "\n",
+    "        self.add_output(name='mass', shape=(1,), desc='cannonball mass', units='kg')\n",
+    "        self.add_output(name='S', shape=(1,), desc='aerodynamic reference area', units='m**2')\n",
+    "\n",
+    "        self.declare_partials(of='mass', wrt='dens')\n",
+    "        self.declare_partials(of='mass', wrt='radius')\n",
+    "\n",
+    "        self.declare_partials(of='S', wrt='radius')\n",
+    "\n",
+    "    def compute(self, inputs, outputs):\n",
+    "        radius = inputs['radius']\n",
+    "        dens = inputs['dens']\n",
+    "\n",
+    "        outputs['mass'] = (4/3.) * dens * np.pi * radius ** 3\n",
+    "        outputs['S'] = np.pi * radius ** 2\n",
+    "\n",
+    "    def compute_partials(self, inputs, partials):\n",
+    "        radius = inputs['radius']\n",
+    "        dens = inputs['dens']\n",
+    "\n",
+    "        partials['mass', 'dens'] = (4/3.) * np.pi * radius ** 3\n",
+    "        partials['mass', 'radius'] = 4. * dens * np.pi * radius ** 2\n",
+    "\n",
+    "        partials['S', 'radius'] = 2 * np.pi * radius\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## The cannonball ODE component\n",
+    "\n",
+    "This component computes the state rates and the kinetic energy of the cannonball.\n",
+    "By calling the `declare_coloring` method wrt all inputs and using method `'cs'`, we're telling OpenMDAO to automatically determine the sparsity pattern of the outputs with respect to the inputs, **and** to automatically compute those outputs using complex-step approximation."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "class CannonballODE(om.ExplicitComponent):\n",
+    "    \"\"\"\n",
+    "    Cannonball ODE assuming flat earth and accounting for air resistance\n",
+    "    \"\"\"\n",
+    "\n",
+    "    def initialize(self):\n",
+    "        self.options.declare('num_nodes', types=int)\n",
+    "\n",
+    "    def setup(self):\n",
+    "        nn = self.options['num_nodes']\n",
+    "\n",
+    "        # static parameters\n",
+    "        self.add_input('m', units='kg')\n",
+    "        self.add_input('S', units='m**2')\n",
+    "        # 0.5 good assumption for a sphere\n",
+    "        self.add_input('CD', 0.5)\n",
+    "\n",
+    "        # time varying inputs\n",
+    "        self.add_input('h', units='m', shape=nn)\n",
+    "        self.add_input('v', units='m/s', shape=nn)\n",
+    "        self.add_input('gam', units='rad', shape=nn)\n",
+    "\n",
+    "        # state rates\n",
+    "        self.add_output('v_dot', shape=nn, units='m/s**2', tags=['dymos.state_rate_source:v'])\n",
+    "        self.add_output('gam_dot', shape=nn, units='rad/s', tags=['dymos.state_rate_source:gam'])\n",
+    "        self.add_output('h_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:h'])\n",
+    "        self.add_output('r_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:r'])\n",
+    "        self.add_output('ke', shape=nn, units='J')\n",
+    "\n",
+    "        # Ask OpenMDAO to compute the partial derivatives using complex-step\n",
+    "        # with a partial coloring algorithm for improved performance, and use\n",
+    "        # a graph coloring algorithm to automatically detect the sparsity pattern.\n",
+    "        self.declare_coloring(wrt='*', method='cs')\n",
+    "\n",
+    "        alt_data = USatm1976Data.alt * om.unit_conversion('ft', 'm')[0]\n",
+    "        rho_data = USatm1976Data.rho * om.unit_conversion('slug/ft**3', 'kg/m**3')[0]\n",
+    "        self.rho_interp = interp1d(np.array(alt_data, dtype=complex),\n",
+    "                                   np.array(rho_data, dtype=complex),\n",
+    "                                   kind='linear')\n",
+    "\n",
+    "    def compute(self, inputs, outputs):\n",
+    "\n",
+    "        gam = inputs['gam']\n",
+    "        v = inputs['v']\n",
+    "        h = inputs['h']\n",
+    "        m = inputs['m']\n",
+    "        S = inputs['S']\n",
+    "        CD = inputs['CD']\n",
+    "\n",
+    "        GRAVITY = 9.80665  # m/s**2\n",
+    "\n",
+    "        # handle complex-step gracefully from the interpolant\n",
+    "        if np.iscomplexobj(h):\n",
+    "            rho = self.rho_interp(inputs['h'])\n",
+    "        else:\n",
+    "            rho = self.rho_interp(inputs['h']).real\n",
+    "\n",
+    "        q = 0.5*rho*inputs['v']**2\n",
+    "        qS = q * S\n",
+    "        D = qS * CD\n",
+    "        cgam = np.cos(gam)\n",
+    "        sgam = np.sin(gam)\n",
+    "        outputs['v_dot'] = - D/m-GRAVITY*sgam\n",
+    "        outputs['gam_dot'] = -(GRAVITY/v)*cgam\n",
+    "        outputs['h_dot'] = v*sgam\n",
+    "        outputs['r_dot'] = v*cgam\n",
+    "        outputs['ke'] = 0.5*m*v**2"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Generating an initial guess\n",
+    "\n",
+    "Since we will be using a nonlinear solver to converge the dynamics and the duration, a simple linear initial guess will not be sufficient. The initial guess we will use is generated from the kinematics equations that assume no atmospheric drag. We compute the state values for some duration and initial flight path angle\n"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def initial_guess(t_dur, gam0=np.pi/3):\n",
+    "    t = np.linspace(0, t_dur, num=100)\n",
+    "    g = -9.80665\n",
+    "    v0 = -0.5*g*t_dur/np.sin(gam0)\n",
+    "    r = v0*np.cos(gam0)*t\n",
+    "    h = v0*np.sin(gam0)*t + 0.5*g*t**2\n",
+    "    v = np.sqrt(v0*np.cos(gam0)**2 + (v0*np.sin(gam0) + g*t)**2)\n",
+    "    gam = np.arctan2(v0*np.sin(gam0) + g*t, v0*np.cos(gam0)**2)\n",
+    "\n",
+    "    guess = {'t': t, 'r': r, 'h': h, 'v': v, 'gam': gam}\n",
+    "\n",
+    "    return guess\n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Building and running the problem\n",
+    "\n",
+    "The following code defines the components for the physical\n",
+    "cannonball calculations and ODE problem, sets up trajectory using a\n",
+    "single phase, and specifies the stopping condition for the phase.\n",
+    "The initial flight path angle is free, since 45 degrees is not\n",
+    "necessarily optimal once air resistance is taken into account."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "p = om.Problem(model=om.Group())\n",
+    "\n",
+    "p.driver = om.pyOptSparseDriver()\n",
+    "p.driver.options['optimizer'] = 'IPOPT'\n",
+    "p.driver.declare_coloring()\n",
+    "\n",
+    "p.driver.opt_settings['derivative_test'] = 'first-order'\n",
+    "p.driver.opt_settings['mu_strategy'] = 'monotone'\n",
+    "p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'\n",
+    "p.driver.opt_settings['bound_mult_init_method'] = 'mu-based'\n",
+    "p.driver.opt_settings['mu_init'] = 0.01\n",
+    "p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'\n",
+    "\n",
+    "p.set_solver_print(level=0, depth=99)\n",
+    "\n",
+    "p.model.add_subsystem('size_comp', CannonballSizeComp(),\n",
+    "                      promotes_inputs=['radius', 'dens'])\n",
+    "p.model.set_input_defaults('dens', val=7.87, units='g/cm**3')\n",
+    "p.model.add_design_var('radius', lower=0.01, upper=0.10,\n",
+    "                       ref0=0.01, ref=0.10, units='m')\n",
+    "\n",
+    "traj = p.model.add_subsystem('traj', dm.Trajectory())\n",
+    "\n",
+    "transcription = dm.Radau(num_segments=25, order=3, compressed=False)\n",
+    "phase = dm.Phase(ode_class=CannonballODE, transcription=transcription)\n",
+    "\n",
+    "phase = traj.add_phase('phase', phase)\n",
+    "\n",
+    "# Duration is bounded to be greater than 1 to ensure the nonlinear solve\n",
+    "# does not converge to the initial point.\n",
+    "phase.set_time_options(fix_initial=True, duration_bounds=(1, 100), units='s')\n",
+    "\n",
+    "# All initial states except flight path angle are fixed\n",
+    "# The output of the ODE which provides the rate source for each state\n",
+    "# is obtained from the tags used on those outputs in the ODE.\n",
+    "# The units of the states are automatically inferred by multiplying the units\n",
+    "# of those rates by the time units.\n",
+    "phase.add_state('r', fix_initial=True, solve_segments='forward')\n",
+    "phase.add_state('h', fix_initial=True, solve_segments='forward')\n",
+    "phase.add_state('gam', fix_initial=False, solve_segments='forward')\n",
+    "phase.add_state('v', fix_initial=False, solve_segments='forward')\n",
+    "\n",
+    "phase.add_parameter('S', units='m**2', static_target=True)\n",
+    "phase.add_parameter('m', units='kg', static_target=True)\n",
+    "phase.add_parameter('CD', val=0.5, opt=False, static_target=True)\n",
+    "\n",
+    "phase.add_boundary_constraint('ke', loc='initial',\n",
+    "                              upper=400000, lower=0, ref=100000)\n",
+    "\n",
+    "# A duration balance is added setting altitude to zero.\n",
+    "# A nonlinear solver is used to find the duration of required to satisfy the condition.\n",
+    "phase.set_duration_balance('h', val=0.0)\n",
+    "\n",
+    "phase.add_objective('r', loc='final', scaler=-1.0)\n",
+    "\n",
+    "p.model.connect('size_comp.mass', 'traj.phase.parameters:m')\n",
+    "p.model.connect('size_comp.S', 'traj.phase.parameters:S')\n",
+    "\n",
+    "# Finish Problem Setup\n",
+    "p.setup()\n",
+    "\n",
+    "#############################################\n",
+    "# Set constants and initial guesses\n",
+    "#############################################\n",
+    "p.set_val('radius', 0.05, units='m')\n",
+    "p.set_val('dens', 7.87, units='g/cm**3')\n",
+    "\n",
+    "p.set_val('traj.phase.parameters:CD', 0.5)\n",
+    "\n",
+    "guess = initial_guess(t_dur=30.0)\n",
+    "\n",
+    "p.set_val('traj.phase.t_initial', 0.0)\n",
+    "p.set_val('traj.phase.t_duration', 30.0)\n",
+    "\n",
+    "p.set_val('traj.phase.states:r', phase.interp('r', ys=guess['r'], xs=guess['t']))\n",
+    "p.set_val('traj.phase.states:h', phase.interp('h', ys=guess['h'], xs=guess['t']))\n",
+    "p.set_val('traj.phase.states:v', phase.interp('v', ys=guess['v'], xs=guess['t']))\n",
+    "p.set_val('traj.phase.states:gam', phase.interp('gam', ys=guess['gam'], xs=guess['t']), units='rad')\n",
+    "\n",
+    "#####################################################\n",
+    "# Run the optimization and final explicit simulation\n",
+    "#####################################################\n",
+    "dm.run_problem(p)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Plotting the results"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "exp_out = traj.simulate()\n",
+    "\n",
+    "rad = p.get_val('radius', units='m')[0]\n",
+    "print(f'optimal radius: {rad} m ')\n",
+    "mass = p.get_val('size_comp.mass', units='kg')[0]\n",
+    "print(f'cannonball mass: {mass} kg ')\n",
+    "angle = p.get_val('traj.phase.timeseries.states:gam', units='deg')[0, 0]\n",
+    "print(f'launch angle: {angle} deg')\n",
+    "max_range = p.get_val('traj.phase.timeseries.states:r')[-1, 0]\n",
+    "print(f'maximum range: {max_range} m')\n",
+    "\n",
+    "fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(10, 6))\n",
+    "\n",
+    "time_imp = p.get_val('traj.phase.timeseries.time')\n",
+    "\n",
+    "time_exp = exp_out.get_val('traj.phase.timeseries.time')\n",
+    "\n",
+    "r_imp = p.get_val('traj.phase.timeseries.states:r')\n",
+    "\n",
+    "r_exp = exp_out.get_val('traj.phase.timeseries.states:r')\n",
+    "\n",
+    "h_imp = p.get_val('traj.phase.timeseries.states:h')\n",
+    "\n",
+    "h_exp = exp_out.get_val('traj.phase.timeseries.states:h')\n",
+    "\n",
+    "axes.plot(r_imp, h_imp, 'bo')\n",
+    "\n",
+    "axes.plot(r_exp, h_exp, 'r--')\n",
+    "\n",
+    "axes.set_xlabel('range (m)')\n",
+    "axes.set_ylabel('altitude (m)')\n",
+    "\n",
+    "fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(10, 6))\n",
+    "states = ['r', 'h', 'v', 'gam']\n",
+    "for i, state in enumerate(states):\n",
+    "    x_imp = p.get_val(f'traj.phase.timeseries.states:{state}')\n",
+    "\n",
+    "    x_exp = exp_out.get_val(f'traj.phase.timeseries.states:{state}')\n",
+    "\n",
+    "    axes[i].set_ylabel(state)\n",
+    "\n",
+    "    axes[i].plot(time_imp, x_imp, 'bo')\n",
+    "    axes[i].plot(time_exp, x_exp, 'r--')\n",
+    "\n",
+    "plt.show()"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "tags": [
+     "remove-input",
+     "remove-output"
+    ]
+   },
+   "outputs": [],
+   "source": [
+    "from openmdao.utils.assert_utils import assert_near_equal\n",
+    "\n",
+    "assert_near_equal(p.get_val('traj.phase.states:r')[-1],\n",
+    "                  3183.25, tolerance=1.0E-2)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "celltoolbar": "Tags",
+  "jupytext": {
+   "cell_metadata_filter": "-all",
+   "notebook_metadata_filter": "-all",
+   "text_representation": {
+    "extension": ".md",
+    "format_name": "markdown"
+   }
+  },
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.11.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 4
+}
diff --git a/dymos/phase/options.py b/dymos/phase/options.py
index 2fc316b8..4d53e0c9 100644
--- a/dymos/phase/options.py
+++ b/dymos/phase/options.py
@@ -571,6 +571,9 @@ class TimeOptionsDictionary(om.OptionsDictionary):
         self.declare(name='t_duration_targets', allow_none=True, default=[],
                      desc='targets in the ODE to which the total duration of the phase is connected')
 
+        self.declare(name='t_duration_balance_options', default={},
+                     desc='options dictionary for the duration residual')
+
 
 class GridRefinementOptionsDictionary(om.OptionsDictionary):
     """
diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index 1e9a9c76..2b93461f 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -710,7 +710,7 @@ class Phase(om.Group):
                                ref=_unspecified, targets=_unspecified, rate_targets=_unspecified,
                                rate2_targets=_unspecified, shape=_unspecified):
         """
-        Adds an polynomial control variable to be tied to a parameter in the ODE.
+        Adds a polynomial control variable to be tied to a parameter in the ODE.
 
         Polynomial controls are defined by values at the Legendre-Gauss-Lobatto nodes of a
         single polynomial, defined on [-1, 1] in phase tau space.
@@ -1707,6 +1707,67 @@ class Phase(om.Group):
         if name is not _unspecified:
             self.time_options['name'] = name
 
+    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.
+
+        Parameters
+        ----------
+        name : str
+            Name of the variable.  This should be a state or control variable, the path to an output
+            from the top level of the RHS, or an expression to be evaluated. If an expression,
+            it must be provided in the form of an equation with a left- and right-hand side.
+        val :  float
+            The value that the residual must equal at the end  point of the phase.
+        index : int, optional
+            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.
+        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}
+
+        expr_operators = ['(', '+', '-', '/', '*', '&', '%', '@']
+        if '=' in name:
+            is_expr = True
+        elif '=' not in name and any(opr in name for opr in expr_operators):
+            raise ValueError(f'The expression provided `{name}` has invalid format. '
+                             'Expression may be a single variable or an equation '
+                             'of the form `constraint_name = func(vars)`')
+        else:
+            is_expr = False
+
+        balance_name = name.split('=')[0].strip() if is_expr else name
+
+        options['is_expr'] = is_expr
+        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)
+
     def classify_var(self, var):
         """
         Classifies a variable of the given name or path.
@@ -1786,6 +1847,9 @@ class Phase(om.Group):
         transcription.setup_ode(self)
 
         transcription.setup_timeseries_outputs(self)
+
+        transcription.setup_duration_balance(self)
+
         transcription.setup_defects(self)
         transcription.setup_solvers(self)
 
@@ -1845,6 +1909,8 @@ class Phase(om.Group):
 
         transcription.configure_timeseries_outputs(self)
 
+        transcription.configure_duration_balance(self)
+
         transcription.configure_solvers(self)
 
     def check_time_options(self):
diff --git a/dymos/transcriptions/analytic/analytic.py b/dymos/transcriptions/analytic/analytic.py
index 8b5f5524..76493a79 100644
--- a/dymos/transcriptions/analytic/analytic.py
+++ b/dymos/transcriptions/analytic/analytic.py
@@ -247,6 +247,28 @@ class Analytic(TranscriptionBase):
         """
         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.
@@ -260,7 +282,7 @@ class Analytic(TranscriptionBase):
 
     def configure_solvers(self, phase):
         """
-        Setup the solvers for this transcription.
+        Configure the solvers for this transcription.
 
         Parameters
         ----------
diff --git a/dymos/transcriptions/explicit_shooting/explicit_shooting.py b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
index 8fdba24c..a22fcab3 100644
--- a/dymos/transcriptions/explicit_shooting/explicit_shooting.py
+++ b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
@@ -643,6 +643,28 @@ class ExplicitShooting(TranscriptionBase):
         """
         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/pseudospectral/components/state_independents.py b/dymos/transcriptions/pseudospectral/components/state_independents.py
index 45a8ffb3..2c3e3460 100644
--- a/dymos/transcriptions/pseudospectral/components/state_independents.py
+++ b/dymos/transcriptions/pseudospectral/components/state_independents.py
@@ -12,7 +12,7 @@ class StateIndependentsComp(om.ImplicitComponent):
     """
     Class definition for the StateIndependentsComp.
 
-    A simple component that replaces the state indepvarcomps whenver the solver needs to solve for
+    A simple component that replaces the state IndepVarComps whenever the solver needs to solve for
     the state or whenever the initial state is connected to an external source.
 
     Parameters
diff --git a/dymos/transcriptions/pseudospectral/components/time_duration_comp.py b/dymos/transcriptions/pseudospectral/components/time_duration_comp.py
new file mode 100644
index 00000000..50041c58
--- /dev/null
+++ b/dymos/transcriptions/pseudospectral/components/time_duration_comp.py
@@ -0,0 +1,211 @@
+import numpy as np
+import openmdao.api as om
+
+from dymos.transcriptions.grid_data import GridData
+from dymos._options import options as dymos_options
+
+
+class TimeDurationComp(om.ImplicitComponent):
+    """
+    Class definition for the TimeDurationComp.
+
+    A simple component that computes the duration of the phase when the solver needs to solve for
+    the stopping condition.
+
+    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')
+
+    def configure_io(self, state_idx_map=None):
+        """
+        I/O creation is delayed until configure so we can determine shape and units.
+
+        Parameters
+        ----------
+        state_idx_map : dict
+            The state_idx_map generated by the transcription during its configuration.  This dict
+            is keyed by state, which in turn holds a dict of the solve indices in the input state
+            vector (keyed as 'solved') and those indices managed by the optimizer (keyed as 'indep').
+            In some unittest, this map is not provided at the time of setup and is instead called
+            later.
+        """
+        if state_idx_map is None:
+            return
+        self.state_idx_map = state_idx_map
+        state_options = self.options['state_options']
+        grid_data = self.options['grid_data']
+        num_col_nodes = grid_data.subset_num_nodes['col']
+
+        num_state_input_nodes = grid_data.subset_num_nodes['state_input']
+
+        self.var_names = {}
+        for state_name in state_options:
+            self.var_names[state_name] = {
+                'defect': f'defects:{state_name}',
+            }
+
+        for state_name, options in state_options.items():
+
+            shape = options['shape']
+            units = options['units']
+            solved = options['solve_segments']
+            default_val = options['val']
+            if np.isscalar(default_val) or np.asarray(default_val).shape == 1:
+                default_val = float(default_val)
+            elif np.asarray(default_val).shape == shape:
+                default_val = np.repeat(default_val[np.newaxis, ...], num_state_input_nodes, axis=0)
+            var_names = self.var_names[state_name]
+
+            if solved and (options['lower'] or options['upper']):
+                om.issue_warning(f'State {state_name} has bounds but they are not enforced when '
+                                 f'using `solve_segments.` Apply a path constraint to {state_name} '
+                                 f'to enforce bounds.', om.UnusedOptionWarning)
+
+            # only need the implicit variable if this state is solved.
+            # Note: we don't add scaling and bounds here. This may be revisited.
+            self.add_output(name=f'states:{state_name}',
+                            shape=(num_state_input_nodes,) + shape,
+                            val=default_val,
+                            units=units)
+
+            # Input for continuity, which can come from an external source.
+            if options['input_initial']:
+                input_name = f'initial_states:{state_name}'
+                self.add_input(name=input_name, shape=(1, ) + shape, units=units)
+
+            # compute an output constraint value since the optimizer needs it
+            if solved:
+                self.add_input(
+                    name=var_names['defect'],
+                    shape=(num_col_nodes, ) + shape,
+                    desc=f'Constraint value for interior defects of state {state_name}',
+                    units=units)
+
+        # Setup partials
+        for state_name, options in state_options.items():
+            shape = options['shape']
+            size = np.prod(shape)
+            solved = options['solve_segments']
+            state_var_name = f'states:{state_name}'
+
+            if solved:  # only need this deriv if its solved
+                solve_idx = np.array(state_idx_map[state_name]['solver'])
+                indep_idx = np.array(state_idx_map[state_name]['indep'])
+
+                num_indep_nodes = indep_idx.shape[0]
+                num_solve_nodes = solve_idx.shape[0]
+
+                base_idx = np.tile(np.arange(size), num_indep_nodes).reshape(num_indep_nodes, size)
+                row = (indep_idx[:, np.newaxis]*size + base_idx).flatten()
+
+                # anything that looks like an indep
+                self.declare_partials(of=state_var_name, wrt=state_var_name,
+                                      rows=row, cols=row, val=-1.0)
+
+                if options['input_initial']:
+                    wrt = f'initial_states:{state_name}'
+                    row_col = np.arange(np.prod(shape))
+                    self.declare_partials(of=state_var_name, wrt=wrt, rows=row_col, cols=row_col,
+                                          val=1.0)
+
+                col = np.arange(num_solve_nodes * size)
+                base_idx = np.tile(np.arange(size), num_solve_nodes).reshape(num_solve_nodes, size)
+                row = (solve_idx[:, np.newaxis]*size + base_idx).flatten()
+
+                var_names = self.var_names[state_name]
+                self.declare_partials(of=state_var_name,
+                                      wrt=var_names['defect'],
+                                      rows=row, cols=col, val=1.0)
+
+            else:
+                row_col = np.arange(num_state_input_nodes*np.prod(shape), dtype=int)
+                self.declare_partials(of=state_var_name, wrt=state_var_name,
+                                      rows=row_col, cols=row_col, val=-1.0)
+
+                if options['input_initial']:
+                    wrt = f'initial_states:{state_name}'
+                    row_col = np.arange(np.prod(shape))
+                    self.declare_partials(of=state_var_name, wrt=wrt, rows=row_col, cols=row_col,
+                                          val=1.0)
+
+    def apply_nonlinear(self, inputs, outputs, residuals):
+        """
+        Compute residuals given inputs and outputs.
+
+        The model is assumed to be in an unscaled state.
+
+        Parameters
+        ----------
+        inputs : Vector
+            Unscaled, dimensional input variables read via inputs[key].
+        outputs : Vector
+            Unscaled, dimensional output variables read via outputs[key].
+        residuals : Vector
+            Unscaled, dimensional residuals written to via residuals[key].
+        """
+        state_options = self.options['state_options']
+
+        for state_name, options in state_options.items():
+            state_var_name = f'states:{state_name}'
+
+            solve_idx = self.state_idx_map[state_name]['solver']
+            indep_idx = self.state_idx_map[state_name]['indep']
+
+            var_names = self.var_names[state_name]
+
+            if options['solve_segments']:
+                defect = inputs[var_names['defect']]
+                residuals[state_var_name][solve_idx, ...] = defect
+
+            # really is: <idep_val> - \outputs[state_name][indep_idx] but OpenMDAO
+            # implementation details mean we just set it to 0
+            # but derivatives are still based on (<idep_val> - \outputs[state_name][indep_idx]),
+            # so you get -1 wrt state var
+            # NOTE: check_partials will report wrong derivs for the indep vars,
+            #       but don't believe it!
+            residuals[state_var_name][indep_idx, ...] = 0.0
+
+            if options['input_initial']:
+                ic_state_name = f'initial_states:{state_name}'
+
+                residuals[state_var_name][0, ...] = \
+                    inputs[ic_state_name][0, ...] - \
+                    outputs[state_var_name][0, ...]
+
+    def solve_nonlinear(self, inputs, outputs):
+        """
+        Compute outputs given inputs.
+
+        The model is assumed to be in an unscaled state.
+
+        Parameters
+        ----------
+        inputs : Vector
+            Unscaled, dimensional input variables read via inputs[key].
+        outputs : Vector
+            Unscaled, dimensional output variables read via outputs[key].
+        """
+        state_options = self.options['state_options']
+
+        for state_name, options in state_options.items():
+            if options['input_initial']:
+                output_name = f'states:{state_name}'
+                input_name = f'initial_states:{state_name}'
+                outputs[output_name][0, ...] = inputs[input_name]
diff --git a/dymos/transcriptions/pseudospectral/pseudospectral_base.py b/dymos/transcriptions/pseudospectral/pseudospectral_base.py
index ba3db3bf..86ed875f 100644
--- a/dymos/transcriptions/pseudospectral/pseudospectral_base.py
+++ b/dymos/transcriptions/pseudospectral/pseudospectral_base.py
@@ -6,7 +6,7 @@ from ..common import TimeComp, TimeseriesOutputGroup
 from .components import StateIndependentsComp, StateInterpComp, CollocationComp, \
     PseudospectralTimeseriesOutputComp
 from ...utils.misc import CoerceDesvar, get_rate_units, reshape_val
-from ...utils.introspection import get_promoted_vars, get_source_metadata
+from ...utils.introspection import get_promoted_vars, get_source_metadata, configure_duration_balance_introspection
 from ...utils.constants import INF_BOUND
 from ...utils.indexing import get_src_indices_by_row
 
@@ -436,6 +436,53 @@ class PseudospectralBase(TranscriptionBase):
                               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']
+
+            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=options['val'],
+                                              eq_units=options['units'], units=phase.time_options['units'],
+                                              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.
@@ -456,7 +503,7 @@ class PseudospectralBase(TranscriptionBase):
         phase : dymos.Phase
             The phase object to which this transcription instance applies.
         """
-        if self.any_solved_segs:
+        if self.any_solved_segs or self._implicit_duration:
             # Only override the solvers if the user hasn't set them to something else.
             if isinstance(phase.nonlinear_solver, om.NonlinearRunOnce):
                 newton = phase.nonlinear_solver = om.NewtonSolver()
@@ -467,7 +514,7 @@ class PseudospectralBase(TranscriptionBase):
                 newton.linesearch = om.BoundsEnforceLS()
 
         # even though you don't need a nl_solver for connections, you still ln_solver since its implicit
-        if self.any_solved_segs or self.any_connected_opt_segs:
+        if self.any_solved_segs or self.any_connected_opt_segs or self._implicit_duration:
             if isinstance(phase.linear_solver, om.LinearRunOnce):
                 phase.linear_solver = om.DirectSolver()
 
diff --git a/dymos/transcriptions/transcription_base.py b/dymos/transcriptions/transcription_base.py
index 360b090d..68123c66 100644
--- a/dymos/transcriptions/transcription_base.py
+++ b/dymos/transcriptions/transcription_base.py
@@ -24,6 +24,7 @@ class TranscriptionBase(object):
     """
     def __init__(self, **kwargs):
 
+        self._implicit_duration = False
         self.grid_data = None
 
         self.options = om.OptionsDictionary()
@@ -80,9 +81,16 @@ class TranscriptionBase(object):
         # Warn about invalid options
         phase.check_time_options()
 
-        if not time_options['input_initial'] or not time_options['input_duration']:
+        if phase.time_options['t_duration_balance_options']:
+            self._implicit_duration = True
+
+        if not time_options['input_initial']:
             phase.add_subsystem('time_extents', om.IndepVarComp(),
                                 promotes_outputs=['*'])
+        else:
+            if not time_options['input_duration'] and not self._implicit_duration:
+                phase.add_subsystem('time_extents', om.IndepVarComp(),
+                                    promotes_outputs=['*'])
 
         for ts_name, ts_options in phase._timeseries.items():
             if t_name not in ts_options['outputs']:
@@ -120,7 +128,7 @@ class TranscriptionBase(object):
         if not time_options['input_initial']:
             indeps.append('t_initial')
 
-        if not time_options['input_duration']:
+        if not time_options['input_duration'] and not self._implicit_duration:
             indeps.append('t_duration')
 
         for var in indeps:
@@ -379,6 +387,31 @@ class TranscriptionBase(object):
         """
         raise NotImplementedError(f'Transcription {self.__class__.__name__} does not implement method setup_ode.')
 
+    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.
+        """
+
+        raise NotImplementedError(f'Transcription {self.__class__.__name__} does not implement'
+                                  f' method setup_duration_balance.')
+
+    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.
+        """
+        raise NotImplementedError(f'Transcription {self.__class__.__name__} does not implement'
+                                  f' method setup_duration_balance.')
+
     def setup_solvers(self, phase):
         """
         Setup the solvers for this transcription.
diff --git a/dymos/utils/introspection.py b/dymos/utils/introspection.py
index 1bc1776a..0120e660 100644
--- a/dymos/utils/introspection.py
+++ b/dymos/utils/introspection.py
@@ -1061,6 +1061,122 @@ def get_target_metadata(ode, name, user_targets=_unspecified, user_units=_unspec
     return shape, units, static_target
 
 
+def configure_duration_balance_introspection(phase):
+    """
+    Modify duration balance 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.
+    """
+    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['balance_name'] if options['is_expr'] else 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 in ['indep_polynomial_control', 'input_polynomial_control']:
+        prefix = 'polynomial_controls:' if dymos_options['use_timeseries_prefix'] else ''
+        control_shape = phase.polynomial_control_options[var]['shape']
+        control_units = phase.polynomial_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 == 'polynomial_control_rate':
+        prefix = 'polynomial_control_rates:' if dymos_options['use_timeseries_prefix'] else ''
+        control_name = var[:-5]
+        control_shape = phase.polynomial_control_options[control_name]['shape']
+        control_units = phase.polynomial_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 == 'polynomial_control_rate2':
+        prefix = 'polynomial_control_rates:' if dymos_options['use_timeseries_prefix'] else ''
+        control_name = var[:-6]
+        control_shape = phase.polynomial_control_options[control_name]['shape']
+        control_units = phase.polynomial_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}'
+
+    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=options['units'])
+
+        options['shape'] = meta['shape']
+        options['units'] = meta['units']
+        options['var_path'] = f'timeseries.{options["balance_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')
+
+
 def _get_targets_metadata(ode, name, user_targets=_unspecified, user_units=_unspecified,
                           user_shape=_unspecified, control_rate=False, user_static_target=_unspecified):
     """


Test patch:
diff --git a/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py b/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py
new file mode 100644
index 00000000..f2adecaf
--- /dev/null
+++ b/dymos/examples/cannonball/doc/test_doc_cannonball_implicit_duration.py
@@ -0,0 +1,294 @@
+import unittest
+
+import matplotlib.pyplot as plt
+import numpy as np
+
+from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
+
+
+plt.switch_backend('Agg')
+
+
+def initial_guess(t_dur, gam0=np.pi/3):
+    t = np.linspace(0, t_dur, num=100)
+    g = -9.80665
+    v0 = -0.5*g*t_dur/np.sin(gam0)
+    r = v0*np.cos(gam0)*t
+    h = v0*np.sin(gam0)*t + 0.5*g*t**2
+    v = np.sqrt(v0*np.cos(gam0)**2 + (v0*np.sin(gam0) + g*t)**2)
+    gam = np.arctan2(v0*np.sin(gam0) + g*t, v0*np.cos(gam0)**2)
+
+    guess = {'t': t, 'r': r, 'h': h, 'v': v, 'gam': gam}
+
+    return guess
+
+
+@use_tempdirs
+class TestTwoPhaseCannonballForDocs(unittest.TestCase):
+
+    @require_pyoptsparse(optimizer='IPOPT')
+    def test_two_phase_cannonball_for_docs(self):
+        from scipy.interpolate import interp1d
+
+        import openmdao.api as om
+        from openmdao.utils.assert_utils import assert_near_equal
+
+        import dymos as dm
+        from dymos.models.atmosphere.atmos_1976 import USatm1976Data
+
+        #############################################
+        # Component for the design part of the model
+        #############################################
+        class CannonballSizeComp(om.ExplicitComponent):
+            """
+            Compute the area and mass of a cannonball with a given radius and density.
+
+            Notes
+            -----
+            This component is not vectorized with 'num_nodes' as is the usual way
+            with Dymos, but is instead intended to compute a scalar mass and reference
+            area from scalar radius and density inputs. This component does not reside
+            in the ODE but instead its outputs are connected to the trajectory via
+            input design parameters.
+            """
+            def setup(self):
+                self.add_input(name='radius', val=1.0, desc='cannonball radius', units='m')
+                self.add_input(name='dens', val=7870., desc='cannonball density', units='kg/m**3')
+
+                self.add_output(name='mass', shape=(1,), desc='cannonball mass', units='kg')
+                self.add_output(name='S', shape=(1,), desc='aerodynamic reference area', units='m**2')
+
+                self.declare_partials(of='mass', wrt='dens')
+                self.declare_partials(of='mass', wrt='radius')
+
+                self.declare_partials(of='S', wrt='radius')
+
+            def compute(self, inputs, outputs):
+                radius = inputs['radius']
+                dens = inputs['dens']
+
+                outputs['mass'] = (4/3.) * dens * np.pi * radius ** 3
+                outputs['S'] = np.pi * radius ** 2
+
+            def compute_partials(self, inputs, partials):
+                radius = inputs['radius']
+                dens = inputs['dens']
+
+                partials['mass', 'dens'] = (4/3.) * np.pi * radius ** 3
+                partials['mass', 'radius'] = 4. * dens * np.pi * radius ** 2
+
+                partials['S', 'radius'] = 2 * np.pi * radius
+
+        #############################################
+        # Build the ODE class
+        #############################################
+        class CannonballODE(om.ExplicitComponent):
+            """
+            Cannonball ODE assuming flat earth and accounting for air resistance
+            """
+
+            def initialize(self):
+                self.options.declare('num_nodes', types=int)
+
+            def setup(self):
+                nn = self.options['num_nodes']
+
+                # static parameters
+                self.add_input('m', units='kg')
+                self.add_input('S', units='m**2')
+                # 0.5 good assumption for a sphere
+                self.add_input('CD', 0.5)
+
+                # time varying inputs
+                self.add_input('h', units='m', shape=nn)
+                self.add_input('v', units='m/s', shape=nn)
+                self.add_input('gam', units='rad', shape=nn)
+
+                # state rates
+                self.add_output('v_dot', shape=nn, units='m/s**2', tags=['dymos.state_rate_source:v'])
+                self.add_output('gam_dot', shape=nn, units='rad/s', tags=['dymos.state_rate_source:gam'])
+                self.add_output('h_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:h'])
+                self.add_output('r_dot', shape=nn, units='m/s', tags=['dymos.state_rate_source:r'])
+                self.add_output('ke', shape=nn, units='J')
+
+                # Ask OpenMDAO to compute the partial derivatives using complex-step
+                # with a partial coloring algorithm for improved performance, and use
+                # a graph coloring algorithm to automatically detect the sparsity pattern.
+                self.declare_coloring(wrt='*', method='cs')
+
+                alt_data = USatm1976Data.alt * om.unit_conversion('ft', 'm')[0]
+                rho_data = USatm1976Data.rho * om.unit_conversion('slug/ft**3', 'kg/m**3')[0]
+                self.rho_interp = interp1d(np.array(alt_data, dtype=complex),
+                                           np.array(rho_data, dtype=complex),
+                                           kind='linear')
+
+            def compute(self, inputs, outputs):
+
+                gam = inputs['gam']
+                v = inputs['v']
+                h = inputs['h']
+                m = inputs['m']
+                S = inputs['S']
+                CD = inputs['CD']
+
+                GRAVITY = 9.80665  # m/s**2
+
+                # handle complex-step gracefully from the interpolant
+                if np.iscomplexobj(h):
+                    rho = self.rho_interp(inputs['h'])
+                else:
+                    rho = self.rho_interp(inputs['h']).real
+
+                q = 0.5*rho*inputs['v']**2
+                qS = q * S
+                D = qS * CD
+                cgam = np.cos(gam)
+                sgam = np.sin(gam)
+                outputs['v_dot'] = - D/m-GRAVITY*sgam
+                outputs['gam_dot'] = -(GRAVITY/v)*cgam
+                outputs['h_dot'] = v*sgam
+                outputs['r_dot'] = v*cgam
+                outputs['ke'] = 0.5*m*v**2
+
+        #############################################
+        # Setup the Dymos problem
+        #############################################
+
+        p = om.Problem(model=om.Group())
+
+        p.driver = om.pyOptSparseDriver()
+        p.driver.options['optimizer'] = 'IPOPT'
+        p.driver.declare_coloring()
+
+        p.driver.opt_settings['derivative_test'] = 'first-order'
+        p.driver.opt_settings['mu_strategy'] = 'monotone'
+        p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'
+        p.driver.opt_settings['bound_mult_init_method'] = 'mu-based'
+        p.driver.opt_settings['mu_init'] = 0.01
+        p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'
+
+        p.set_solver_print(level=0, depth=99)
+
+        p.model.add_subsystem('size_comp', CannonballSizeComp(),
+                              promotes_inputs=['radius', 'dens'])
+        p.model.set_input_defaults('dens', val=7.87, units='g/cm**3')
+        p.model.add_design_var('radius', lower=0.01, upper=0.10,
+                               ref0=0.01, ref=0.10, units='m')
+
+        traj = p.model.add_subsystem('traj', dm.Trajectory())
+
+        transcription = dm.Radau(num_segments=25, order=3, compressed=False)
+        phase = dm.Phase(ode_class=CannonballODE, transcription=transcription)
+
+        phase = traj.add_phase('phase', phase)
+
+        # All initial states except flight path angle are fixed
+        # The output of the ODE which provides the rate source for each state
+        # is obtained from the tags used on those outputs in the ODE.
+        # The units of the states are automatically inferred by multiplying the units
+        # of those rates by the time units.
+        phase.set_time_options(fix_initial=True, duration_bounds=(1, 100), units='s')
+        phase.add_state('r', fix_initial=True, solve_segments='forward')
+        phase.add_state('h', fix_initial=True, solve_segments='forward')
+        phase.add_state('gam', fix_initial=False, solve_segments='forward')
+        phase.add_state('v', fix_initial=False, solve_segments='forward')
+
+        phase.add_parameter('S', units='m**2', static_target=True)
+        phase.add_parameter('m', units='kg', static_target=True)
+        phase.add_parameter('CD', val=0.5, opt=False, static_target=True)
+
+        phase.add_boundary_constraint('ke', loc='initial',
+                                      upper=400000, lower=0, ref=100000)
+
+        # A duration balance is added setting altitude to zero.
+        # A nonlinear solver is used to find the duration of required to satisfy the condition.
+        # The duration was bounded to be greater than 1 to ensure the solver did not
+        # converge to the initial point.
+        phase.set_duration_balance('h', val=0.0)
+
+        phase.add_objective('r', loc='final', scaler=-1.0)
+
+        p.model.connect('size_comp.mass', 'traj.phase.parameters:m')
+        p.model.connect('size_comp.S', 'traj.phase.parameters:S')
+
+        # Finish Problem Setup
+        p.setup()
+
+        #############################################
+        # Set constants and initial guesses
+        #############################################
+        p.set_val('radius', 0.05, units='m')
+        p.set_val('dens', 7.87, units='g/cm**3')
+
+        p.set_val('traj.phase.parameters:CD', 0.5)
+
+        guess = initial_guess(t_dur=30.0)
+
+        p.set_val('traj.phase.t_initial', 0.0)
+        p.set_val('traj.phase.t_duration', 30.0)
+
+        p.set_val('traj.phase.states:r', phase.interp('r', ys=guess['r'], xs=guess['t']))
+        p.set_val('traj.phase.states:h', phase.interp('h', ys=guess['h'], xs=guess['t']))
+        p.set_val('traj.phase.states:v', phase.interp('v', ys=guess['v'], xs=guess['t']))
+        p.set_val('traj.phase.states:gam', phase.interp('gam', ys=guess['gam'], xs=guess['t']), units='rad')
+
+        #####################################################
+        # Run the optimization and final explicit simulation
+        #####################################################
+        dm.run_problem(p)
+
+        assert_near_equal(p.get_val('traj.phase.states:r')[-1],
+                          3183.25, tolerance=1.0)
+
+        exp_out = traj.simulate()
+
+        #############################################
+        # Plot the results
+        #############################################
+        rad = p.get_val('radius', units='m')[0]
+        print(f'optimal radius: {rad} m ')
+        mass = p.get_val('size_comp.mass', units='kg')[0]
+        print(f'cannonball mass: {mass} kg ')
+        angle = p.get_val('traj.phase.timeseries.states:gam', units='deg')[0, 0]
+        print(f'launch angle: {angle} deg')
+        max_range = p.get_val('traj.phase.timeseries.states:r')[-1, 0]
+        print(f'maximum range: {max_range} m')
+
+        fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(10, 6))
+
+        time_imp = p.get_val('traj.phase.timeseries.time')
+
+        time_exp = exp_out.get_val('traj.phase.timeseries.time')
+
+        r_imp = p.get_val('traj.phase.timeseries.states:r')
+
+        r_exp = exp_out.get_val('traj.phase.timeseries.states:r')
+
+        h_imp = p.get_val('traj.phase.timeseries.states:h')
+
+        h_exp = exp_out.get_val('traj.phase.timeseries.states:h')
+
+        axes.plot(r_imp, h_imp, 'bo')
+
+        axes.plot(r_exp, h_exp, 'b--')
+
+        axes.set_xlabel('range (m)')
+        axes.set_ylabel('altitude (m)')
+
+        fig, axes = plt.subplots(nrows=4, ncols=1, figsize=(10, 6))
+        states = ['r', 'h', 'v', 'gam']
+        for i, state in enumerate(states):
+            x_imp = p.get_val(f'traj.phase.timeseries.states:{state}')
+
+            x_exp = exp_out.get_val(f'traj.phase.timeseries.states:{state}')
+
+            axes[i].set_ylabel(state)
+
+            axes[i].plot(time_imp, x_imp, 'bo')
+            axes[i].plot(time_exp, x_exp, 'b--')
+
+        plt.show()
+
+
+if __name__ == '__main__':  # pragma: no cover
+    unittest.main()
diff --git a/dymos/phase/test/test_implicit_duration.py b/dymos/phase/test/test_implicit_duration.py
new file mode 100644
index 00000000..53c98833
--- /dev/null
+++ b/dymos/phase/test/test_implicit_duration.py
@@ -0,0 +1,205 @@
+import openmdao.api as om
+import dymos as dm
+import numpy as np
+
+import unittest
+from openmdao.utils.assert_utils import assert_near_equal
+from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
+
+
+class ODEComp(om.ExplicitComponent):
+
+    def initialize(self):
+        self.options.declare('num_nodes', types=int)
+
+    def setup(self):
+        nn = self.options['num_nodes']
+        self.add_input('h', shape=nn, units='m')
+        self.add_input('v', shape=nn, units='m/s')
+        self.add_output('hdot', shape=nn, units='m/s')
+        self.add_output('vdot', shape=nn, units='m/s**2')
+
+        self.declare_partials(of='hdot', wrt='v', rows=np.arange(nn), cols=np.arange(nn), val=1.0)
+
+    def compute(self, inputs, outputs):
+        outputs['hdot'] = inputs['v']
+        outputs['vdot'] = -9.80665
+
+
+class MatrixStateCannonball(om.ExplicitComponent):
+
+    def initialize(self):
+        self.options.declare('num_nodes', types=int)
+
+    def setup(self):
+        nn = self.options['num_nodes']
+        # z is the state vector, a nn x 2 x 2 in the form of [[x, y], [vx, vy]]
+        self.add_input('z', shape=(nn, 2, 2), units=None)
+        self.add_output('zdot', shape=(nn, 2, 2), units=None)
+
+        self.declare_partials(of='zdot', wrt='z', method='cs')
+        self.declare_coloring(wrt=['z'], method='cs', num_full_jacs=5, tol=1.0E-12)
+
+    def compute(self, inputs, outputs):
+        outputs['zdot'][:, 0, 0] = inputs['z'][:, 1, 0]
+        outputs['zdot'][:, 0, 1] = inputs['z'][:, 1, 1]
+        outputs['zdot'][:, 1, 0] = 0.0
+        outputs['zdot'][:, 1, 1] = -9.81
+
+
+@use_tempdirs
+class TestImplicitDuration(unittest.TestCase):
+
+    @staticmethod
+    def _make_simple_problem(tx, input_duration=False, fix_duration=False, expr_end=False):
+        p = om.Problem()
+
+        traj = dm.Trajectory()
+        phase = dm.Phase(ode_class=ODEComp, transcription=tx)
+        traj.add_phase('phase', phase)
+        p.model.add_subsystem('traj', traj)
+
+        phase.set_time_options(fix_initial=True, units='s',
+                               input_duration=input_duration, fix_duration=fix_duration)
+        phase.add_state('h', rate_source='hdot', fix_initial=True, units='m', solve_segments='forward')
+        phase.add_state('v', rate_source='vdot', fix_initial=True, units='m/s', solve_segments='forward')
+
+        if not expr_end:
+            phase.set_duration_balance('h', val=0.0, units='m')
+        else:
+            phase.set_duration_balance('pe=9.80665*h', val=0.0, units='m**2/s**2')
+
+        phase.set_simulate_options(rtol=1.0E-9, atol=1.0E-9)
+
+        p.driver = om.ScipyOptimizeDriver(optimizer='SLSQP')
+        p.driver.declare_coloring(tol=1.0E-12)
+
+        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]))
+
+        return p
+
+    @staticmethod
+    def _make_matrix_problem(tx, shape_error=False):
+        p = om.Problem()
+
+        traj = dm.Trajectory()
+        phase = dm.Phase(ode_class=MatrixStateCannonball, transcription=tx)
+        traj.add_phase('phase', phase)
+        p.model.add_subsystem('traj', traj)
+
+        phase.set_time_options(fix_initial=True, units='s')
+        phase.add_state('z', rate_source='zdot', fix_initial=True, solve_segments='forward')
+
+        index = None if shape_error else [[0], [1]]
+        phase.set_duration_balance('z', val=0.0, index=index)
+
+        phase.set_simulate_options(rtol=1.0E-9, atol=1.0E-9)
+
+        p.driver = om.ScipyOptimizeDriver(optimizer='SLSQP')
+        p.driver.declare_coloring(tol=1.0E-12)
+
+        p.setup()
+
+        p.set_val('traj.phase.t_initial', 0)
+        p.set_val('traj.phase.t_duration', 7)
+        p.set_val('traj.phase.states:z', phase.interp('z', [[[0, 0], [10, 10]], [[10, 0], [10, -10]]]))
+
+        return p
+
+    def test_implicit_duration_radau(self):
+        tx = dm.Radau(num_segments=12, order=3, solve_segments=False)
+
+        p = self._make_simple_problem(tx)
+
+        dm.run_problem(p, run_driver=False, simulate=False)
+
+        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.states:h')[-1], 0.0, tolerance=1E-6)
+
+    def test_implicit_duration_gl(self):
+        tx = dm.GaussLobatto(num_segments=12, order=3, solve_segments=False)
+
+        p = self._make_simple_problem(tx)
+
+        dm.run_problem(p, run_driver=False, simulate=False)
+
+        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.states:h')[-1], 0.0, tolerance=1E-6)
+
+    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'
+
+        with self.assertRaises(ValueError) as e:
+            self._make_matrix_problem(tx, shape_error=True)
+
+        self.assertEqual(expected, str(e.exception))
+
+    def test_implicit_duration_matrix_state(self):
+        tx = dm.Radau(num_segments=12, order=3, solve_segments=False)
+
+        p = self._make_matrix_problem(tx, shape_error=False)
+
+        dm.run_problem(p, run_driver=False, simulate=False)
+
+        assert_near_equal(p.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
+
+    @require_pyoptsparse(optimizer='IPOPT')
+    def test_implicit_duration_radau_expr_condition(self):
+        tx = dm.Radau(num_segments=12, order=3, solve_segments=False)
+
+        p = self._make_simple_problem(tx, expr_end=True)
+
+        dm.run_problem(p, run_driver=False, simulate=False)
+
+        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.states:h')[-1], 0.0, tolerance=1E-6)
+
+        print((p.get_val('traj.phase.timeseries.states:v')[-1]) ** 2 / 2)
+
+    def test_implicit_duration_gl_expr_condition(self):
+        tx = dm.GaussLobatto(num_segments=12, order=3, solve_segments=False)
+
+        p = self._make_simple_problem(tx, expr_end=True)
+
+        dm.run_problem(p, run_driver=False, simulate=False)
+
+        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.states:h')[-1], 0.0, tolerance=1E-6)