From 3f895497646a29f6db2084c54286f297b7ebc450 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Fri, 7 Oct 2022 15:24:46 -0500
Subject: [PATCH 01/34] added test problems

---
 dymos/phase/test/test_timeseries.py | 254 +++++++++++++++++++++++++++-
 1 file changed, 252 insertions(+), 2 deletions(-)

diff --git a/dymos/phase/test/test_timeseries.py b/dymos/phase/test/test_timeseries.py
index 0d87db5c4..8c2b87448 100644
--- a/dymos/phase/test/test_timeseries.py
+++ b/dymos/phase/test/test_timeseries.py
@@ -15,7 +15,7 @@
 from dymos.examples.min_time_climb.aero import AeroGroup
 from dymos.examples.min_time_climb.prop import PropGroup
 from dymos.models.eom import FlightPathEOM2D
-
+from dymos.utils.testing_utils import assert_timeseries_near_equal
 
 @use_tempdirs
 class TestTimeseriesOutput(unittest.TestCase):
@@ -241,7 +241,7 @@ def setup(self):
 
 @require_pyoptsparse(optimizer='SLSQP')
 def min_time_climb(num_seg=3, transcription_class=dm.Radau, transcription_order=3,
-                   force_alloc_complex=False):
+                   force_alloc_complex=False, timeseries_expr=False):
 
     p = om.Problem(model=om.Group())
 
@@ -304,6 +304,8 @@ def min_time_climb(num_seg=3, transcription_class=dm.Radau, transcription_order=
 
     # Add all ODE outputs to the timeseries
     phase.add_timeseries_output('*')
+    if timeseries_expr:
+        phase.add_timeseries_output('R = atmos.pres / (atmos.temp * atmos.rho)')
 
     p.model.linear_solver = om.DirectSolver()
 
@@ -343,5 +345,253 @@ def test_duplicate_timeseries_glob_name_gl(self):
         self.assertEqual(str(e.exception), msg)
 
 
+class TestTimeseriesExprBrachistochrone(unittest.TestCase):
+
+    @staticmethod
+    def make_problem_brachistochrone(transcription, polynomial_control=False):
+        p = om.Problem(model=om.Group())
+
+        p.driver = om.ScipyOptimizeDriver()
+        p.driver.declare_coloring()
+
+        phase = dm.Phase(ode_class=BrachistochroneODE,
+                         transcription=transcription)
+
+        p.model.add_subsystem('phase0', phase)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+
+        phase.set_state_options('x', fix_initial=True)
+        phase.set_state_options('y', fix_initial=True)
+        phase.set_state_options('v', fix_initial=True)
+
+        if polynomial_control:
+            phase.add_polynomial_control('theta', order=1, units='deg', lower=0.01, upper=179.9)
+            control_name = 'polynomial_controls:theta'
+        else:
+            phase.add_control('theta', continuity=True, rate_continuity=True, opt=True,
+                              units='deg', lower=0.01, upper=179.9, ref=1, ref0=0)
+            control_name = 'controls:theta'
+        phase.add_boundary_constraint('x', loc='final', equals=10.0)
+        phase.add_boundary_constraint('y', loc='final', equals=5.0)
+
+        phase.add_parameter('g', opt=True, units='m/s**2', val=9.80665, include_timeseries=True)
+
+        phase.add_objective('time_phase', loc='final', scaler=10)
+        phase.add_timeseries_output('z=states:x*states:y', units='m**2')
+        phase.add_timeseries_output(f'f=3*parameter_vals:g*cos({control_name})**2', units='deg**2')
+
+        p.model.options['assembled_jac_type'] = 'csc'
+        p.model.linear_solver = om.DirectSolver()
+        p.setup(check=True)
+
+        p['phase0.t_initial'] = 0.0
+        p['phase0.t_duration'] = 2.0
+
+        p['phase0.states:x'] = phase.interp('x', [0, 10])
+        p['phase0.states:y'] = phase.interp('y', [10, 5])
+        p['phase0.states:v'] = phase.interp('v', [0, 9.9])
+        p[f'phase0.{control_name}'] = phase.interp('theta', [5, 100])
+        p['phase0.parameters:g'] = 9.80665
+        return p
+
+    def test_timeseries_expr_radau(self):
+        tx = dm.Radau(num_segments=5, order=3, compressed=True)
+        p = self.make_problem_brachistochrone(transcription=tx)
+        p.run_driver()
+
+        x = p.get_val('phase0.timeseries.states:x')
+        y = p.get_val('phase0.timeseries.states:y')
+        theta = p.get_val('phase0.timeseries.controls:theta')
+        g = p.get_val('phase0.timeseries.parameters:g')
+
+        z_computed = x * y
+        f_computed = 3 * g * np.cos(theta)**2
+        z_ts = p.get_val('phase0.timeseries.z')
+        f_ts = p.get_val('phase0.timeseries.f')
+        assert_near_equal(z_computed, z_ts, tolerance=1e-12)
+        assert_near_equal(f_computed, f_ts, tolerance=1e-12)
+
+    def test_timeseries_expr_radau_polynomial_control(self):
+        tx = dm.Radau(num_segments=5, order=3, compressed=True)
+        p = self.make_problem_brachistochrone(transcription=tx, polynomial_control=True)
+        p.run_driver()
+
+        x = p.get_val('phase0.timeseries.states:x')
+        y = p.get_val('phase0.timeseries.states:y')
+        theta = p.get_val('phase0.timeseries.polynomial_controls:theta')
+        g = p.get_val('phase0.timeseries.parameters:g')
+
+        z_computed = x * y
+        f_computed = 3 * g * np.cos(theta)**2
+        z_ts = p.get_val('phase0.timeseries.z')
+        f_ts = p.get_val('phase0.timeseries.f')
+        assert_near_equal(z_computed, z_ts, tolerance=1e-12)
+        assert_near_equal(f_computed, f_ts, tolerance=1e-12)
+
+    def test_timeseries_expr_gl(self):
+        tx = dm.GaussLobatto(num_segments=5, order=3, compressed=True)
+        p = self.make_problem_brachistochrone(transcription=tx)
+        p.run_driver()
+        x = p.get_val('phase0.timeseries.states:x')
+        y = p.get_val('phase0.timeseries.states:y')
+        theta = p.get_val('phase0.timeseries.controls:theta')
+        g = p.get_val('phase0.timeseries.parameters:g')
+
+        z_computed = x * y
+        f_computed = 3 * g * np.cos(theta) ** 2
+
+        z_ts = p.get_val('phase0.timeseries.z')
+        f_ts = p.get_val('phase0.timeseries.f')
+        assert_near_equal(z_computed, z_ts, tolerance=1e-12)
+        assert_near_equal(f_computed, f_ts, tolerance=1e-12)
+
+    def test_timeseries_expr_gl_polynomial_control(self):
+        tx = dm.GaussLobatto(num_segments=5, order=3, compressed=True)
+        p = self.make_problem_brachistochrone(transcription=tx, polynomial_control=True)
+        p.run_driver()
+
+        x = p.get_val('phase0.timeseries.states:x')
+        y = p.get_val('phase0.timeseries.states:y')
+        theta = p.get_val('phase0.timeseries.polynomial_controls:theta')
+        g = p.get_val('phase0.timeseries.parameters:g')
+
+        z_computed = x * y
+        f_computed = 3 * g * np.cos(theta)**2
+        z_ts = p.get_val('phase0.timeseries.z')
+        f_ts = p.get_val('phase0.timeseries.f')
+        assert_near_equal(z_computed, z_ts, tolerance=1e-12)
+        assert_near_equal(f_computed, f_ts, tolerance=1e-12)
+
+    def test_timeseries_expr_explicit_shooting(self):
+        tx = dm.ExplicitShooting(num_segments=3, grid='gauss-lobatto',
+                                 method='rk4', order=5,
+                                 num_steps_per_segment=5,
+                                 compressed=True)
+
+        p = self.make_problem_brachistochrone(transcription=tx)
+        p.run_driver()
+        x = p.get_val('phase0.timeseries.states:x')
+        y = p.get_val('phase0.timeseries.states:y')
+        theta = p.get_val('phase0.timeseries.controls:theta')
+        g = p.get_val('phase0.timeseries.parameters:g')
+
+        z_computed = x * y
+        f_computed = 3 * g * np.cos(theta) ** 2
+
+        z_ts = p.get_val('phase0.timeseries.z')
+        f_ts = p.get_val('phase0.timeseries.f')
+        assert_near_equal(z_computed, z_ts, tolerance=1e-12)
+        assert_near_equal(f_computed, f_ts, tolerance=1e-12)
+
+    def test_timeseries_expr_explicit_shooting_polynomial_controls(self):
+        tx = dm.ExplicitShooting(num_segments=3, grid='gauss-lobatto',
+                                 method='rk4', order=5,
+                                 num_steps_per_segment=5,
+                                 compressed=True)
+
+        p = self.make_problem_brachistochrone(transcription=tx, polynomial_control=True)
+        p.run_driver()
+        x = p.get_val('phase0.timeseries.states:x')
+        y = p.get_val('phase0.timeseries.states:y')
+        theta = p.get_val('phase0.timeseries.polynomial_controls:theta')
+        g = p.get_val('phase0.timeseries.parameters:g')
+
+        z_computed = x * y
+        f_computed = 3 * g * np.cos(theta) ** 2
+
+        z_ts = p.get_val('phase0.timeseries.z')
+        f_ts = p.get_val('phase0.timeseries.f')
+        assert_near_equal(z_computed, z_ts, tolerance=1e-12)
+        assert_near_equal(f_computed, f_ts, tolerance=1e-12)
+
+    def test_timeseries_expr_solve_ivp(self):
+        tx = dm.Radau(num_segments=5, order=3, compressed=True)
+
+        p = self.make_problem_brachistochrone(transcription=tx)
+        p.run_driver()
+        x = p.get_val('phase0.timeseries.states:x')
+        y = p.get_val('phase0.timeseries.states:y')
+        theta = p.get_val('phase0.timeseries.controls:theta')
+        g = p.get_val('phase0.timeseries.parameters:g')
+        t = p.get_val('phase0.timeseries.time')
+
+        phase0 = p.model._get_subsystem('phase0')
+        sim = phase0.simulate()
+
+        z_computed = x * y
+        f_computed = 3 * g * np.cos(theta) ** 2
+
+        z_sim = sim.get_val('phase0.timeseries.z')
+        f_sim = sim.get_val('phase0.timeseries.f')
+        t_sim = sim.get_val('phase0.timeseries.time')
+        assert_timeseries_near_equal(t, z_computed, t_sim, z_sim, tolerance=1e-3)
+        assert_timeseries_near_equal(t, f_computed, t_sim, f_sim, tolerance=1e-3)
+
+    def test_timeseries_expr_solve_ivp_polynomial_controls(self):
+        tx = dm.Radau(num_segments=5, order=3, compressed=True)
+
+        p = self.make_problem_brachistochrone(transcription=tx, polynomial_control=True)
+        p.run_driver()
+        x = p.get_val('phase0.timeseries.states:x')
+        y = p.get_val('phase0.timeseries.states:y')
+        theta = p.get_val('phase0.timeseries.polynomial_controls:theta')
+        g = p.get_val('phase0.timeseries.parameters:g')
+        t = p.get_val('phase0.timeseries.time')
+
+        phase0 = p.model._get_subsystem('phase0')
+        sim = phase0.simulate()
+
+        z_computed = x * y
+        f_computed = 3 * g * np.cos(theta) ** 2
+
+        z_sim = sim.get_val('phase0.timeseries.z')
+        f_sim = sim.get_val('phase0.timeseries.f')
+        t_sim = sim.get_val('phase0.timeseries.time')
+        assert_timeseries_near_equal(t, z_computed, t_sim, z_sim, tolerance=1e-3)
+        assert_timeseries_near_equal(t, f_computed, t_sim, f_sim, tolerance=1e-3)
+
+
+class TestTimeseriesMinTimeClimb(unittest.TestCase):
+
+    def test_timeseries_expr_radau(self):
+        p = min_time_climb(transcription_class=dm.Radau, num_seg=12, transcription_order=3, timeseries_expr=True)
+        p.run_model()
+
+        pres = p.get_val('traj.phase0.timeseries.pres')
+        temp = p.get_val('traj.phase0.timeseries.temp')
+        rho = p.get_val('traj.phase0.timeseries.rho')
+
+        R_computed = pres / (temp * rho)
+        R_ts = p.get_val('traj.phase0.timeseries.R')
+        assert_near_equal(R_computed, R_ts, tolerance=1e-12)
+
+    def test_timeseries_expr_gl(self):
+        p = min_time_climb(transcription_class=dm.GaussLobatto, num_seg=12,
+                           transcription_order=3, timeseries_expr=True)
+        p.run_driver()
+
+        pres = p.get_val('traj.phase0.timeseries.pres')
+        temp = p.get_val('traj.phase0.timeseries.temp')
+        rho = p.get_val('traj.phase0.timeseries.rho')
+
+        R_computed = pres / (temp * rho)
+        R_ts = p.get_val('traj.phase0.timeseries.R')
+        assert_near_equal(R_computed, R_ts, tolerance=1e-12)
+
+    def test_timeseries_expr_explicit_shooting(self):
+        p = min_time_climb(transcription_class=dm.ExplicitShooting, num_seg=12,
+                           transcription_order=3, timeseries_expr=True)
+        p.run_model()
+
+        pres = p.get_val('traj.phase0.timeseries.pres')
+        temp = p.get_val('traj.phase0.timeseries.temp')
+        rho = p.get_val('traj.phase0.timeseries.rho')
+
+        R_computed = pres / (temp * rho)
+        R_ts = p.get_val('traj.phase0.timeseries.R')
+        assert_near_equal(R_computed, R_ts, tolerance=1e-12)
+
+
 if __name__ == '__main__':  # pragma: no cover
     unittest.main()

From 7d2b03d1ed95ec641d4ef3b6f2dc5cf4094a3bdb Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Fri, 7 Oct 2022 15:35:15 -0500
Subject: [PATCH 02/34] added option to timeseries options dictionary
 indicating if an output is an expression

---
 dymos/phase/options.py | 9 +++++++--
 1 file changed, 7 insertions(+), 2 deletions(-)

diff --git a/dymos/phase/options.py b/dymos/phase/options.py
index 024d9bb72..ee2220ed6 100644
--- a/dymos/phase/options.py
+++ b/dymos/phase/options.py
@@ -690,13 +690,15 @@ def __init__(self, read_only=False):
                           'with wildcards.')
 
         self.declare(name='shape', default=None, allow_none=True,
-                     desc='The shape of the timeseries output variable. This is generally determined automatically by dymos.')
+                     desc='The shape of the timeseries output variable.'
+                          ' This is generally determined automatically by dymos.')
 
         self.declare(name='units', default=None, allow_none=True,
                      desc='Units to be used for the timeseries output, or None to leave the units unchanged.')
 
         self.declare(name='src', types=str, default=None, allow_none=True,
-                     desc='The phase-relative path of the source of this timeseries output, used when issuing connections.')
+                     desc='The phase-relative path of the source of this timeseries output,'
+                          ' used when issuing connections.')
 
         self.declare(name='src_idxs', types=(Iterable,), default=None, allow_none=True,
                      desc='The indices of the source of this timeseries output to be used when issuing connections.')
@@ -706,3 +708,6 @@ def __init__(self, read_only=False):
 
         self.declare(name='is_rate', default=False, allow_none=False,
                      desc='If True this is a rate.')
+
+        self.declare(name='is_expr', default=False, allow_none=False,
+                     desc='If true the requested timeseries is a mathematical expression')

From e4a7781a24b478e8fa00e98ef2b74063318f052d Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Fri, 7 Oct 2022 15:43:00 -0500
Subject: [PATCH 03/34] add new option when adding timeseries outputs

---
 dymos/phase/phase.py | 13 ++++++++++---
 1 file changed, 10 insertions(+), 3 deletions(-)

diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index 913d10218..80c8c0a3f 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -1271,6 +1271,7 @@ def add_timeseries_rate_output(self, name, output_name=None, units=_unspecified,
         """
         if type(name) is list:
             for i, name_i in enumerate(name):
+                expr = True if '=' in name_i else False
                 if type(units) is dict:  # accept dict for units when using array of name
                     unit = units.get(name_i, None)
                 elif type(units) is list:  # allow matching list for units
@@ -1282,21 +1283,24 @@ def add_timeseries_rate_output(self, name, output_name=None, units=_unspecified,
                                                     units=unit,
                                                     shape=shape,
                                                     timeseries=timeseries,
-                                                    rate=True)
+                                                    rate=True,
+                                                    expr=expr)
 
                 # Handle specific units for wildcard names.
                 if oname is not None and '*' in name_i:
                     self._timeseries[timeseries]['outputs'][oname]['wildcard_units'] = units
 
         else:
+            expr = True if '=' in name else False
             self._add_timeseries_output(name, output_name=output_name,
                                         units=units,
                                         shape=shape,
                                         timeseries=timeseries,
-                                        rate=True)
+                                        rate=True,
+                                        expr=expr)
 
     def _add_timeseries_output(self, name, output_name=None, units=_unspecified, shape=_unspecified,
-                               timeseries='timeseries', rate=False):
+                               timeseries='timeseries', rate=False, expr=False):
         r"""
         Add a single variable or rate to the timeseries outputs of the phase.
 
@@ -1337,6 +1341,8 @@ def _add_timeseries_output(self, name, output_name=None, units=_unspecified, sha
 
         if '*' in name:
             output_name = name
+        elif expr:
+            output_name = name.split()
         elif output_name is None:
             output_name = name.rpartition('.')[-1]
 
@@ -1354,6 +1360,7 @@ def _add_timeseries_output(self, name, output_name=None, units=_unspecified, sha
             ts_output['units'] = units
             ts_output['shape'] = shape
             ts_output['is_rate'] = rate
+            ts_output['is_expr'] = expr
 
             self._timeseries[timeseries]['outputs'][output_name] = ts_output
 

From ac0fb3acacd0867d6a43030dcdecb0fbb8210ae1 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Tue, 11 Oct 2022 15:47:43 -0500
Subject: [PATCH 04/34] added a method to introspect expressions and make the
 connections to its input variables

---
 dymos/phase/phase.py                          | 22 +++--
 dymos/phase/test/test_timeseries.py           | 24 ++---
 .../pseudospectral/pseudospectral_base.py     |  9 ++
 dymos/utils/introspection.py                  | 99 +++++++++++++++++--
 4 files changed, 126 insertions(+), 28 deletions(-)

diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index 80c8c0a3f..527b67ad0 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -66,6 +66,7 @@ def __init__(self, from_phase=None, **kwargs):
         self.parameter_options = {}
         self.refine_options = GridRefinementOptionsDictionary()
         self.simulate_options = SimulateOptionsDictionary()
+        self.timeseries_ec_vars = {}
 
         # Dictionaries of variable options that are set by the user via the API
         # These will be applied over any defaults specified by decorators on the ODE
@@ -1218,6 +1219,7 @@ def add_timeseries_output(self, name, output_name=None, units=_unspecified, shap
         """
         if type(name) is list:
             for i, name_i in enumerate(name):
+                expr = True if '=' in name_i else False
                 if type(units) is dict:  # accept dict for units when using array of name
                     unit = units.get(name_i, None)
                 elif type(units) is list:  # allow matching list for units
@@ -1229,18 +1231,21 @@ def add_timeseries_output(self, name, output_name=None, units=_unspecified, shap
                                                     units=unit,
                                                     shape=shape,
                                                     timeseries=timeseries,
-                                                    rate=False)
+                                                    rate=False,
+                                                    expr=expr)
 
                 # Handle specific units for wildcard names.
                 if oname is not None and '*' in name_i:
                     self._timeseries[timeseries]['outputs'][oname]['wildcard_units'] = units
 
         else:
+            expr = True if '=' in name else False
             self._add_timeseries_output(name, output_name=output_name,
                                         units=units,
                                         shape=shape,
                                         timeseries=timeseries,
-                                        rate=False)
+                                        rate=False,
+                                        expr=expr)
 
     def add_timeseries_rate_output(self, name, output_name=None, units=_unspecified, shape=_unspecified,
                                    timeseries='timeseries'):
@@ -1291,13 +1296,11 @@ def add_timeseries_rate_output(self, name, output_name=None, units=_unspecified,
                     self._timeseries[timeseries]['outputs'][oname]['wildcard_units'] = units
 
         else:
-            expr = True if '=' in name else False
             self._add_timeseries_output(name, output_name=output_name,
                                         units=units,
                                         shape=shape,
                                         timeseries=timeseries,
-                                        rate=True,
-                                        expr=expr)
+                                        rate=True)
 
     def _add_timeseries_output(self, name, output_name=None, units=_unspecified, shape=_unspecified,
                                timeseries='timeseries', rate=False, expr=False):
@@ -1339,10 +1342,10 @@ def _add_timeseries_output(self, name, output_name=None, units=_unspecified, sha
         if timeseries not in self._timeseries:
             raise ValueError(f'Timeseries {timeseries} does not exist in phase {self.pathname}')
 
-        if '*' in name:
+        if expr:
+            output_name = name.split('=')[0]
+        elif '*' in name:
             output_name = name
-        elif expr:
-            output_name = name.split()
         elif output_name is None:
             output_name = name.rpartition('.')[-1]
 
@@ -1700,7 +1703,8 @@ def configure(self):
         try:
             configure_timeseries_output_introspection(self)
         except RuntimeError as val_err:
-            raise RuntimeError(f'Error during configure_timeseries_output_introspection in phase {self.pathname}.') from val_err
+            raise RuntimeError(f'Error during configure_timeseries_output_introspection in phase {self.pathname}.')\
+                from val_err
 
         transcription.configure_timeseries_outputs(self)
 
diff --git a/dymos/phase/test/test_timeseries.py b/dymos/phase/test/test_timeseries.py
index 8c2b87448..215e4092b 100644
--- a/dymos/phase/test/test_timeseries.py
+++ b/dymos/phase/test/test_timeseries.py
@@ -231,8 +231,8 @@ def setup(self):
                            subsys=FlightPathEOM2D(num_nodes=nn),
                            promotes_inputs=['m', 'v', 'gam', 'alpha'])
 
-        foo = self.add_subsystem('foo', om.IndepVarComp())
-        foo.add_output('rho', val=100 * np.ones(nn), units='g/cm**3')
+        # foo = self.add_subsystem('foo', om.IndepVarComp())
+        # foo.add_output('rho', val=100 * np.ones(nn), units='g/cm**3')
 
         self.connect('aero.f_drag', 'flight_dynamics.D')
         self.connect('aero.f_lift', 'flight_dynamics.L')
@@ -378,8 +378,8 @@ def make_problem_brachistochrone(transcription, polynomial_control=False):
         phase.add_parameter('g', opt=True, units='m/s**2', val=9.80665, include_timeseries=True)
 
         phase.add_objective('time_phase', loc='final', scaler=10)
-        phase.add_timeseries_output('z=states:x*states:y', units='m**2')
-        phase.add_timeseries_output(f'f=3*parameter_vals:g*cos({control_name})**2', units='deg**2')
+        phase.add_timeseries_output('z=x*y + x**2', units='m**2')
+        phase.add_timeseries_output('f=3*g*cos(theta)**2', units='deg**2')
 
         p.model.options['assembled_jac_type'] = 'csc'
         p.model.linear_solver = om.DirectSolver()
@@ -405,7 +405,7 @@ def test_timeseries_expr_radau(self):
         theta = p.get_val('phase0.timeseries.controls:theta')
         g = p.get_val('phase0.timeseries.parameters:g')
 
-        z_computed = x * y
+        z_computed = x * y + x**2
         f_computed = 3 * g * np.cos(theta)**2
         z_ts = p.get_val('phase0.timeseries.z')
         f_ts = p.get_val('phase0.timeseries.f')
@@ -422,7 +422,7 @@ def test_timeseries_expr_radau_polynomial_control(self):
         theta = p.get_val('phase0.timeseries.polynomial_controls:theta')
         g = p.get_val('phase0.timeseries.parameters:g')
 
-        z_computed = x * y
+        z_computed = x * y + x**2
         f_computed = 3 * g * np.cos(theta)**2
         z_ts = p.get_val('phase0.timeseries.z')
         f_ts = p.get_val('phase0.timeseries.f')
@@ -438,7 +438,7 @@ def test_timeseries_expr_gl(self):
         theta = p.get_val('phase0.timeseries.controls:theta')
         g = p.get_val('phase0.timeseries.parameters:g')
 
-        z_computed = x * y
+        z_computed = x * y + x**2
         f_computed = 3 * g * np.cos(theta) ** 2
 
         z_ts = p.get_val('phase0.timeseries.z')
@@ -456,7 +456,7 @@ def test_timeseries_expr_gl_polynomial_control(self):
         theta = p.get_val('phase0.timeseries.polynomial_controls:theta')
         g = p.get_val('phase0.timeseries.parameters:g')
 
-        z_computed = x * y
+        z_computed = x * y + x**2
         f_computed = 3 * g * np.cos(theta)**2
         z_ts = p.get_val('phase0.timeseries.z')
         f_ts = p.get_val('phase0.timeseries.f')
@@ -476,7 +476,7 @@ def test_timeseries_expr_explicit_shooting(self):
         theta = p.get_val('phase0.timeseries.controls:theta')
         g = p.get_val('phase0.timeseries.parameters:g')
 
-        z_computed = x * y
+        z_computed = x * y + x**2
         f_computed = 3 * g * np.cos(theta) ** 2
 
         z_ts = p.get_val('phase0.timeseries.z')
@@ -497,7 +497,7 @@ def test_timeseries_expr_explicit_shooting_polynomial_controls(self):
         theta = p.get_val('phase0.timeseries.polynomial_controls:theta')
         g = p.get_val('phase0.timeseries.parameters:g')
 
-        z_computed = x * y
+        z_computed = x * y + x**2
         f_computed = 3 * g * np.cos(theta) ** 2
 
         z_ts = p.get_val('phase0.timeseries.z')
@@ -519,7 +519,7 @@ def test_timeseries_expr_solve_ivp(self):
         phase0 = p.model._get_subsystem('phase0')
         sim = phase0.simulate()
 
-        z_computed = x * y
+        z_computed = x * y + x**2
         f_computed = 3 * g * np.cos(theta) ** 2
 
         z_sim = sim.get_val('phase0.timeseries.z')
@@ -542,7 +542,7 @@ def test_timeseries_expr_solve_ivp_polynomial_controls(self):
         phase0 = p.model._get_subsystem('phase0')
         sim = phase0.simulate()
 
-        z_computed = x * y
+        z_computed = x * y + x**2
         f_computed = 3 * g * np.cos(theta) ** 2
 
         z_sim = sim.get_val('phase0.timeseries.z')
diff --git a/dymos/transcriptions/pseudospectral/pseudospectral_base.py b/dymos/transcriptions/pseudospectral/pseudospectral_base.py
index d10884c01..744966b0d 100644
--- a/dymos/transcriptions/pseudospectral/pseudospectral_base.py
+++ b/dymos/transcriptions/pseudospectral/pseudospectral_base.py
@@ -483,6 +483,15 @@ def setup_timeseries_outputs(self, phase):
         gd = self.grid_data
 
         for name, options in phase._timeseries.items():
+            expr_ts = False
+            for _, output_options in options['outputs'].items():
+                if output_options['is_expr']:
+                    ts_exec_comp = om.ExecComp(has_diag_partials=True)
+                    phase.add_subsystem('timeseries_exec_comp', ts_exec_comp)
+                    expr_ts = True
+                if expr_ts:
+                    break
+
             if options['transcription'] is None:
                 ogd = None
             else:
diff --git a/dymos/utils/introspection.py b/dymos/utils/introspection.py
index f69380bbb..d6a8a5e97 100644
--- a/dymos/utils/introspection.py
+++ b/dymos/utils/introspection.py
@@ -2,6 +2,7 @@
 import fnmatch
 
 import openmdao.api as om
+import numpy as np
 from openmdao.utils.array_utils import shape_to_len
 from dymos.utils.misc import _unspecified
 from ..phase.options import StateOptionsDictionary, TimeseriesOutputOptionsDictionary
@@ -726,12 +727,16 @@ def configure_timeseries_output_introspection(phase):
         not_found = set()
 
         for output_name, output_options in ts_opts['outputs'].items():
-            try:
-                output_meta = transcription._get_timeseries_var_source(output_options['name'],
-                                                                       output_options['output_name'],
-                                                                       phase=phase)
-            except ValueError as e:
-                not_found.add(output_name)
+            if output_options['is_expr']:
+                output_meta = configure_timeseries_expr_introspection(phase, output_options)
+
+            else:
+                try:
+                    output_meta = transcription._get_timeseries_var_source(output_options['name'],
+                                                                           output_options['output_name'],
+                                                                           phase=phase)
+                except ValueError as e:
+                    not_found.add(output_name)
 
             output_options['src'] = output_meta['src']
             output_options['src_idxs'] = output_meta['src_idxs']
@@ -751,6 +756,86 @@ def configure_timeseries_output_introspection(phase):
             ts_opts['outputs'].pop(s)
 
 
+def configure_timeseries_expr_introspection(phase, expr_options):
+    """
+    Introspect the timeseries expressions, add expressions to the exec comp, and make connections to the input variables
+
+    Parameters
+    ----------
+    phase : Phase
+        The phase object whose timeseries outputs are to be introspected.
+    expr_options : dict
+        Options for the expression to be added to the timeseries exec comp
+    """
+    import re
+    timeseries_ec = phase._get_subsystem('timeseries_exec_comp')
+    transcription = phase.options['transcription']
+    gd = transcription.grid_data
+    num_output_nodes = gd.subset_num_nodes['all']
+
+    expr = expr_options['name']
+    expr_lhs = expr.split('=')[0]
+    if '.' in expr_lhs or ':' in expr_lhs:
+        raise ValueError(f'{expr_lhs} is an invalid name for the timeseries LHS. Must use a valid variable name')
+    expr_reduced = expr
+    units = expr_options['units'] if expr_options['units'] is not _unspecified else None
+    shape = expr_options['shape'] if expr_options['shape'] is not _unspecified else (1,)
+
+    abs_names = [x.strip() for x in re.findall(re.compile(r'([_a-zA-Z]\w*[ ]*\(?:?[.]?)'), expr)
+                 if not x.endswith('(') and not x.endswith(':')]
+    for name in abs_names:
+        if name.endswith('.'):
+            idx = abs_names.index(name)
+            abs_names[idx + 1] = name + abs_names[idx + 1]
+            abs_names.remove(name)
+    for name in abs_names:
+        var_name = name.split('.')[-1]
+        if var_name not in phase.timeseries_ec_vars.keys():
+            phase.timeseries_ec_vars[var_name] = {'added_to_ec': False, 'abs_name': name}
+    expr_vars = {}
+    expr_kwargs = {}
+    meta = {'units': units, 'shape': shape}
+    for var in phase.timeseries_ec_vars:
+        if phase.timeseries_ec_vars[var]['added_to_ec']:
+            continue
+        elif var == expr_lhs:
+            var_units = units
+            var_shape = shape
+            meta['src'] = f'timeseries_exec_comp.{var}'
+            phase.timeseries_ec_vars[var]['added_to_ec'] = True
+        else:
+            try:
+                expr_vars[var] = transcription._get_timeseries_var_source(var, output_name=var, phase=phase)
+            except ValueError as e:
+                expr_vars[var] = transcription._get_timeseries_var_source(phase.timeseries_ec_vars[var]['abs_name'],
+                                                                          output_name=var, phase=phase)
+
+            var_units = expr_vars[var]['units']
+            var_shape = expr_vars[var]['shape']
+
+        expr_kwargs[var] = {'units': var_units, 'shape': (num_output_nodes,) + var_shape}
+
+    for input_var in expr_vars:
+        abs_name = phase.timeseries_ec_vars[input_var]['abs_name']
+        if '.' in abs_name:
+            expr_reduced = expr_reduced.replace(abs_name, input_var)
+
+    timeseries_ec.add_expr(expr_reduced, **expr_kwargs)
+
+    for var, options in expr_vars.items():
+        if phase.timeseries_ec_vars[var]['added_to_ec']:
+            continue
+        else:
+            phase.connect(src_name=options['src'], tgt_name=f'timeseries_exec_comp.{var}',
+                          src_indices=options['src_idxs'])
+
+        phase.timeseries_ec_vars[var]['added_to_ec'] = True
+
+    meta['src_idxs'] = None
+
+    return meta
+
+
 def filter_outputs(patterns, sys):
     """
     Find all outputs of the given system that match one or more of the strings given in patterns.
@@ -1057,7 +1142,7 @@ def get_source_metadata(ode, src, user_units=_unspecified, user_shape=_unspecifi
     ode_outputs = ode if isinstance(ode, dict) else get_promoted_vars(ode, iotypes='output')
 
     if src not in ode_outputs:
-        raise ValueError(f'Unable to find the source {src} in the ODE.')
+        raise ValueError(f"Unable to find the source '{src}' in the ODE.")
 
     if user_units in {None, _unspecified}:
         meta['units'] = ode_outputs[src]['units']

From e21d70bacb67bd94110eb820af99109a8bd5e510 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Tue, 11 Oct 2022 16:08:16 -0500
Subject: [PATCH 05/34] modified configure method for gl transcription to check
 if a timeseries output is an expression. If it is, then the output name
 option is used instead of the name

---
 dymos/phase/phase.py                                 | 3 ++-
 dymos/transcriptions/pseudospectral/gauss_lobatto.py | 3 ++-
 dymos/utils/introspection.py                         | 8 +++++++-
 3 files changed, 11 insertions(+), 3 deletions(-)

diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index 527b67ad0..f7a387c65 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -1599,7 +1599,8 @@ def classify_var(self, var):
         return classify_var(var, state_options=self.state_options,
                             parameter_options=self.parameter_options,
                             control_options=self.control_options,
-                            polynomial_control_options=self.polynomial_control_options)
+                            polynomial_control_options=self.polynomial_control_options,
+                            timeseries_options=self._timeseries)
 
     def _check_ode(self):
         """
diff --git a/dymos/transcriptions/pseudospectral/gauss_lobatto.py b/dymos/transcriptions/pseudospectral/gauss_lobatto.py
index 1f7473e26..7695142ab 100644
--- a/dymos/transcriptions/pseudospectral/gauss_lobatto.py
+++ b/dymos/transcriptions/pseudospectral/gauss_lobatto.py
@@ -365,7 +365,8 @@ def configure_interleave_comp(self, phase):
 
         for timeseries_name, timeseries_options in phase._timeseries.items():
             for ts_output_name, ts_output in timeseries_options['outputs'].items():
-                var_type = phase.classify_var(ts_output['name'])
+                name = ts_output['name'] if not ts_output['is_expr'] else ts_output_name
+                var_type = phase.classify_var(name)
                 if var_type == 'ode':
                     units = ts_output['units']
                     shape = ts_output['shape']
diff --git a/dymos/utils/introspection.py b/dymos/utils/introspection.py
index d6a8a5e97..6e06f4fcc 100644
--- a/dymos/utils/introspection.py
+++ b/dymos/utils/introspection.py
@@ -9,7 +9,8 @@
 from .misc import get_rate_units
 
 
-def classify_var(var, state_options, parameter_options, control_options, polynomial_control_options):
+def classify_var(var, state_options, parameter_options, control_options,
+                 polynomial_control_options, timeseries_options=None):
     """
     Classifies a variable of the given name or path.
 
@@ -66,6 +67,11 @@ def classify_var(var, state_options, parameter_options, control_options, polynom
             return 'control_rate2'
         elif var[:-6] in polynomial_control_options:
             return 'polynomial_control_rate2'
+    elif timeseries_options is not None:
+        for timeseries in timeseries_options:
+            if var in timeseries_options[timeseries]['outputs']:
+                if timeseries_options[timeseries]['outputs'][var]['is_expr']:
+                    return 'timeseries_exec_comp_output'
 
     return 'ode'
 

From ec5e6bc7f6db32ea20efb3c15b99d4efbd5a8d69 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Wed, 12 Oct 2022 10:15:53 -0500
Subject: [PATCH 06/34] added timeseries exec comp for solve_ivp and explicit
 shooting

---
 .../explicit_shooting/explicit_shooting.py             | 10 +++++++++-
 dymos/transcriptions/solve_ivp/solve_ivp.py            |  2 +-
 2 files changed, 10 insertions(+), 2 deletions(-)

diff --git a/dymos/transcriptions/explicit_shooting/explicit_shooting.py b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
index 80cce487c..87dad464a 100644
--- a/dymos/transcriptions/explicit_shooting/explicit_shooting.py
+++ b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
@@ -361,7 +361,15 @@ def setup_timeseries_outputs(self, phase):
             The phase object to which this transcription instance applies.
         """
         gd = self.grid_data
-        for name in phase._timeseries:
+        for name, options in phase._timeseries.items():
+            expr_ts = False
+            for _, output_options in options['outputs'].items():
+                if output_options['is_expr']:
+                    ts_exec_comp = om.ExecComp(has_diag_partials=True)
+                    phase.add_subsystem('timeseries_exec_comp', ts_exec_comp)
+                    expr_ts = True
+                if expr_ts:
+                    break
             timeseries_comp = ExplicitTimeseriesComp(input_grid_data=gd,
                                                      output_subset='segment_ends')
             phase.add_subsystem(name, subsys=timeseries_comp)
diff --git a/dymos/transcriptions/solve_ivp/solve_ivp.py b/dymos/transcriptions/solve_ivp/solve_ivp.py
index 3759771ba..c985b20e2 100644
--- a/dymos/transcriptions/solve_ivp/solve_ivp.py
+++ b/dymos/transcriptions/solve_ivp/solve_ivp.py
@@ -512,7 +512,7 @@ def setup_timeseries_outputs(self, phase):
         phase.add_subsystem('timeseries', subsys=timeseries_comp)
 
         # Remove all subsequent timeseries
-        for ts_name in list(phase._timeseries.keys()):
+        for ts_name, options in phase._timeseries.items():
             if ts_name != 'timeseries':
                 phase._timeseries.pop(ts_name)
 

From ec6aba5941068ccf8629e7e637bac0378f873350 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Thu, 20 Oct 2022 10:36:30 -0500
Subject: [PATCH 07/34] changed introspection to determine shape from in-built
 methods. Also ensured no white spaces are left in variable names

---
 dymos/phase/phase.py                | 2 +-
 dymos/phase/test/test_timeseries.py | 2 +-
 dymos/utils/introspection.py        | 8 ++++----
 3 files changed, 6 insertions(+), 6 deletions(-)

diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index f7a387c65..f0a5c88fd 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -1343,7 +1343,7 @@ def _add_timeseries_output(self, name, output_name=None, units=_unspecified, sha
             raise ValueError(f'Timeseries {timeseries} does not exist in phase {self.pathname}')
 
         if expr:
-            output_name = name.split('=')[0]
+            output_name = name.split('=')[0].strip()
         elif '*' in name:
             output_name = name
         elif output_name is None:
diff --git a/dymos/phase/test/test_timeseries.py b/dymos/phase/test/test_timeseries.py
index b3d063db2..598417469 100644
--- a/dymos/phase/test/test_timeseries.py
+++ b/dymos/phase/test/test_timeseries.py
@@ -569,7 +569,7 @@ def test_timeseries_expr_solve_ivp(self):
         p.run_driver()
         x = p.get_val('phase0.timeseries.states:x')
         y = p.get_val('phase0.timeseries.states:y')
-        theta = p.get_val('phase0.timeseries.controls:theta')
+        theta = p.get_val('phase0.timeseries.polynomial_controls:theta')
         g = p.get_val('phase0.timeseries.parameters:g')
         t = p.get_val('phase0.timeseries.time')
 
diff --git a/dymos/utils/introspection.py b/dymos/utils/introspection.py
index 6e06f4fcc..3ce2d79aa 100644
--- a/dymos/utils/introspection.py
+++ b/dymos/utils/introspection.py
@@ -1,5 +1,6 @@
 from collections.abc import Iterable
 import fnmatch
+import re
 
 import openmdao.api as om
 import numpy as np
@@ -773,14 +774,13 @@ def configure_timeseries_expr_introspection(phase, expr_options):
     expr_options : dict
         Options for the expression to be added to the timeseries exec comp
     """
-    import re
+
     timeseries_ec = phase._get_subsystem('timeseries_exec_comp')
     transcription = phase.options['transcription']
-    gd = transcription.grid_data
-    num_output_nodes = gd.subset_num_nodes['all']
+    num_output_nodes = transcription._get_num_timeseries_nodes()
 
     expr = expr_options['name']
-    expr_lhs = expr.split('=')[0]
+    expr_lhs = expr.split('=')[0].strip()
     if '.' in expr_lhs or ':' in expr_lhs:
         raise ValueError(f'{expr_lhs} is an invalid name for the timeseries LHS. Must use a valid variable name')
     expr_reduced = expr

From 1650b1c56a44f8620c10754797dbabe4d188f1f9 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Thu, 20 Oct 2022 11:10:19 -0500
Subject: [PATCH 08/34] Added the polynomial controls and their rates as
 timeseries outputs in the explicit shooting transcription

---
 .../explicit_shooting/explicit_shooting.py          | 13 +++++++++++++
 .../test/test_explicit_shooting.py                  |  2 ++
 2 files changed, 15 insertions(+)

diff --git a/dymos/transcriptions/explicit_shooting/explicit_shooting.py b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
index ebc3f5a75..a87b43f9e 100644
--- a/dymos/transcriptions/explicit_shooting/explicit_shooting.py
+++ b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
@@ -255,6 +255,19 @@ def setup_polynomial_controls(self, phase):
             The phase object to which this transcription instance applies.
         """
         phase._check_polynomial_control_options()
+        if phase.control_options:
+            for name, options in phase.polynomial_control_options.items():
+                for ts_name, ts_options in phase._timeseries.items():
+                    if f'polynomial_controls:{name}' not in ts_options['outputs']:
+                        phase.add_timeseries_output(name, output_name=f'polynomial_controls:{name}',
+                                                    timeseries=ts_name)
+                    if f'polynomial_control_rates:{name}_rate' not in ts_options['outputs']:
+                        phase.add_timeseries_output(f'{name}_rate', output_name=f'polynomial_control_rates:{name}_rate',
+                                                    timeseries=ts_name)
+                    if f'polynomial_control_rates:{name}_rate2' not in ts_options['outputs']:
+                        phase.add_timeseries_output(f'{name}_rate2',
+                                                    output_name=f'polynomial_control_rates:{name}_rate2',
+                                                    timeseries=ts_name)
 
     def configure_polynomial_controls(self, phase):
         """
diff --git a/dymos/transcriptions/explicit_shooting/test/test_explicit_shooting.py b/dymos/transcriptions/explicit_shooting/test/test_explicit_shooting.py
index dd90f3f05..66d362bad 100644
--- a/dymos/transcriptions/explicit_shooting/test/test_explicit_shooting.py
+++ b/dymos/transcriptions/explicit_shooting/test/test_explicit_shooting.py
@@ -476,11 +476,13 @@ def test_brachistochrone_explicit_shooting_polynomial_control(self):
         y = prob.get_val('phase0.timeseries.states:y')
         t = prob.get_val('phase0.timeseries.time')
         tp = prob.get_val('phase0.timeseries.time_phase')
+        theta = prob.get_val('phase0.polynomial_controls:theta')
 
         assert_near_equal(x[-1, ...], 10.0, tolerance=1.0E-5)
         assert_near_equal(y[-1, ...], 5.0, tolerance=1.0E-5)
         assert_near_equal(t[-1, ...], 1.8016, tolerance=5.0E-3)
         assert_near_equal(tp[-1, ...], 1.8016, tolerance=5.0E-3)
+        assert_near_equal(theta[-1, ...], 90.2, tolerance=1e-1)
 
         with np.printoptions(linewidth=1024):
             cpd = prob.check_partials(compact_print=False, method='cs', out_stream=None)

From 44508f7097b1c807ec6451257c0e95086c9f8b79 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Thu, 20 Oct 2022 11:46:00 -0500
Subject: [PATCH 09/34] fixed a bug in introspection where units of rate2 was
 incorrect

---
 dymos/transcriptions/explicit_shooting/explicit_shooting.py   | 4 ++--
 .../explicit_shooting/test/test_explicit_shooting.py          | 2 +-
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/dymos/transcriptions/explicit_shooting/explicit_shooting.py b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
index a87b43f9e..abf842994 100644
--- a/dymos/transcriptions/explicit_shooting/explicit_shooting.py
+++ b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
@@ -255,7 +255,7 @@ def setup_polynomial_controls(self, phase):
             The phase object to which this transcription instance applies.
         """
         phase._check_polynomial_control_options()
-        if phase.control_options:
+        if phase.polynomial_control_options:
             for name, options in phase.polynomial_control_options.items():
                 for ts_name, ts_options in phase._timeseries.items():
                     if f'polynomial_controls:{name}' not in ts_options['outputs']:
@@ -667,7 +667,7 @@ def _get_timeseries_var_source(self, var, output_name, phase):
             control_name = var[:-6]
             path = f'integrator.polynomial_control_rates:{control_name}_rate2'
             control = phase.polynomial_control_options[control_name]
-            src_units = get_rate_units(control['units'], time_units, deriv=1)
+            src_units = get_rate_units(control['units'], time_units, deriv=2)
             src_shape = control['shape']
         elif var_type == 'parameter':
             path = f'parameter_vals:{var}'
diff --git a/dymos/transcriptions/explicit_shooting/test/test_explicit_shooting.py b/dymos/transcriptions/explicit_shooting/test/test_explicit_shooting.py
index 66d362bad..21aef7aa5 100644
--- a/dymos/transcriptions/explicit_shooting/test/test_explicit_shooting.py
+++ b/dymos/transcriptions/explicit_shooting/test/test_explicit_shooting.py
@@ -476,7 +476,7 @@ def test_brachistochrone_explicit_shooting_polynomial_control(self):
         y = prob.get_val('phase0.timeseries.states:y')
         t = prob.get_val('phase0.timeseries.time')
         tp = prob.get_val('phase0.timeseries.time_phase')
-        theta = prob.get_val('phase0.polynomial_controls:theta')
+        theta = prob.get_val('phase0.timeseries.polynomial_controls:theta')
 
         assert_near_equal(x[-1, ...], 10.0, tolerance=1.0E-5)
         assert_near_equal(y[-1, ...], 5.0, tolerance=1.0E-5)

From 44a8dffea3fb87eb1bb6ce997eb25881bf10a197 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Fri, 21 Oct 2022 13:24:40 -0500
Subject: [PATCH 10/34] changed order of instantiation for timeseries exec comp
 and timeseries_component

---
 dymos/phase/test/test_timeseries.py         | 29 ++++++++-------------
 dymos/transcriptions/solve_ivp/solve_ivp.py | 27 ++++++++++++++-----
 2 files changed, 32 insertions(+), 24 deletions(-)

diff --git a/dymos/phase/test/test_timeseries.py b/dymos/phase/test/test_timeseries.py
index 598417469..f57578db8 100644
--- a/dymos/phase/test/test_timeseries.py
+++ b/dymos/phase/test/test_timeseries.py
@@ -567,46 +567,39 @@ def test_timeseries_expr_solve_ivp(self):
 
         p = self.make_problem_brachistochrone(transcription=tx)
         p.run_driver()
-        x = p.get_val('phase0.timeseries.states:x')
-        y = p.get_val('phase0.timeseries.states:y')
-        theta = p.get_val('phase0.timeseries.polynomial_controls:theta')
-        g = p.get_val('phase0.timeseries.parameters:g')
+
         t = p.get_val('phase0.timeseries.time')
+        z_ts = p.get_val('phase0.timeseries.z')
+        f_ts = p.get_val('phase0.timeseries.f')
 
         phase0 = p.model._get_subsystem('phase0')
         sim = phase0.simulate()
 
-        z_computed = x * y + x**2
-        f_computed = 3 * g * np.cos(theta) ** 2
-
         z_sim = sim.get_val('phase0.timeseries.z')
         f_sim = sim.get_val('phase0.timeseries.f')
         t_sim = sim.get_val('phase0.timeseries.time')
-        assert_timeseries_near_equal(t, z_computed, t_sim, z_sim, tolerance=1e-3)
-        assert_timeseries_near_equal(t, f_computed, t_sim, f_sim, tolerance=1e-3)
+
+        assert_timeseries_near_equal(t, z_ts, t_sim, z_sim, abs_tolerance=1e-1)
+        assert_timeseries_near_equal(t, f_ts, t_sim, f_sim, abs_tolerance=1e-1)
 
     def test_timeseries_expr_solve_ivp_polynomial_controls(self):
         tx = dm.Radau(num_segments=5, order=3, compressed=True)
 
         p = self.make_problem_brachistochrone(transcription=tx, polynomial_control=True)
         p.run_driver()
-        x = p.get_val('phase0.timeseries.states:x')
-        y = p.get_val('phase0.timeseries.states:y')
-        theta = p.get_val('phase0.timeseries.polynomial_controls:theta')
-        g = p.get_val('phase0.timeseries.parameters:g')
+
         t = p.get_val('phase0.timeseries.time')
+        z_ts = p.get_val('phase0.timeseries.z')
+        f_ts = p.get_val('phase0.timeseries.f')
 
         phase0 = p.model._get_subsystem('phase0')
         sim = phase0.simulate()
 
-        z_computed = x * y + x**2
-        f_computed = 3 * g * np.cos(theta) ** 2
-
         z_sim = sim.get_val('phase0.timeseries.z')
         f_sim = sim.get_val('phase0.timeseries.f')
         t_sim = sim.get_val('phase0.timeseries.time')
-        assert_timeseries_near_equal(t, z_computed, t_sim, z_sim, tolerance=1e-3)
-        assert_timeseries_near_equal(t, f_computed, t_sim, f_sim, tolerance=1e-3)
+        assert_timeseries_near_equal(t, z_ts, t_sim, z_sim, abs_tolerance=1e-3)
+        assert_timeseries_near_equal(t, f_ts, t_sim, f_sim, abs_tolerance=1e-3)
 
 
 class TestTimeseriesMinTimeClimb(unittest.TestCase):
diff --git a/dymos/transcriptions/solve_ivp/solve_ivp.py b/dymos/transcriptions/solve_ivp/solve_ivp.py
index c985b20e2..7cca34916 100644
--- a/dymos/transcriptions/solve_ivp/solve_ivp.py
+++ b/dymos/transcriptions/solve_ivp/solve_ivp.py
@@ -503,6 +503,16 @@ def setup_timeseries_outputs(self, phase):
             The phase object to which this transcription instance applies.
         """
         gd = self.grid_data
+        # Remove all timeseries other than 'timeseries'
+        for ts_name, options in phase._timeseries.items():
+            if ts_name != 'timeseries':
+                phase._timeseries.pop(ts_name)
+                continue
+            ts_exec_comp = om.ExecComp(has_diag_partials=True)
+            phase.add_subsystem('timeseries_exec_comp', ts_exec_comp)
+            expr_ts = True
+            if expr_ts:
+                break
 
         timeseries_comp = \
             SolveIVPTimeseriesOutputComp(input_grid_data=gd,
@@ -511,11 +521,6 @@ def setup_timeseries_outputs(self, phase):
 
         phase.add_subsystem('timeseries', subsys=timeseries_comp)
 
-        # Remove all subsequent timeseries
-        for ts_name, options in phase._timeseries.items():
-            if ts_name != 'timeseries':
-                phase._timeseries.pop(ts_name)
-
     def get_parameter_connections(self, name, phase):
         """
         Returns info about a parameter's target connections in the phase.
@@ -588,7 +593,6 @@ def _get_timeseries_var_source(self, var, output_name, phase):
             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']
 
@@ -669,3 +673,14 @@ def _get_timeseries_var_source(self, var, output_name, phase):
         meta['shape'] = src_shape
 
         return meta
+
+    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_segments * self.options['output_nodes_per_seg']

From 81cbcb6e74167dd5a14e0b4617edc111862c25df Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Fri, 21 Oct 2022 14:12:28 -0500
Subject: [PATCH 11/34] Added capability for analytic phases

---
 dymos/phase/test/test_analytic_phase.py   | 29 +++++++++++++++++++++++
 dymos/transcriptions/analytic/analytic.py |  8 +++++++
 2 files changed, 37 insertions(+)

diff --git a/dymos/phase/test/test_analytic_phase.py b/dymos/phase/test/test_analytic_phase.py
index 18490d968..802406281 100644
--- a/dymos/phase/test/test_analytic_phase.py
+++ b/dymos/phase/test/test_analytic_phase.py
@@ -185,6 +185,35 @@ def test_set_polynomial_control_options(self):
 
         self.assertEqual('AnalyticPhase does not support polynomial controls.', str(e.exception))
 
+    def test_timeseries_expr(self):
+        p = om.Problem()
+        traj = p.model.add_subsystem('traj', dm.Trajectory())
+
+        phase = dm.AnalyticPhase(ode_class=SimpleIVPSolution, num_nodes=10)
+
+        traj.add_phase('phase', phase)
+
+        phase.set_time_options(units='s', targets=['t'], fix_initial=True, fix_duration=True)
+        phase.add_state('y')
+        phase.add_parameter('y0', opt=False)
+        phase.add_timeseries_output('z = y0 + y**2')
+
+        p.setup()
+
+        p.set_val('traj.phase.t_initial', 0.0, units='s')
+        p.set_val('traj.phase.t_duration', 2.0, units='s')
+        p.set_val('traj.phase.parameters:y0', 0.5, units='unitless')
+
+        p.run_model()
+
+        y = p.get_val('traj.phase.timeseries.states:y', units='unitless')
+        y0 = p.get_val('traj.phase.timeseries.parameters:y0', units='unitless')
+
+        expected_z = y0 + y**2
+        z = p.get_val('traj.phase.timeseries.z')
+
+        assert_near_equal(z, expected_z)
+
 
 class TestLinkedAnalyticPhases(unittest.TestCase):
 
diff --git a/dymos/transcriptions/analytic/analytic.py b/dymos/transcriptions/analytic/analytic.py
index a89e483fb..38b8a5da6 100644
--- a/dymos/transcriptions/analytic/analytic.py
+++ b/dymos/transcriptions/analytic/analytic.py
@@ -266,6 +266,14 @@ def setup_timeseries_outputs(self, phase):
         gd = self.grid_data
 
         for name, options in phase._timeseries.items():
+            expr_ts = False
+            for _, output_options in options['outputs'].items():
+                if output_options['is_expr']:
+                    ts_exec_comp = om.ExecComp(has_diag_partials=True)
+                    phase.add_subsystem('timeseries_exec_comp', ts_exec_comp)
+                    expr_ts = True
+                if expr_ts:
+                    break
             if options['transcription'] is None:
                 ogd = None
             else:

From 7d23cc038857771d5d20900d03d53a40471f0469 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Wed, 26 Oct 2022 10:50:23 -0500
Subject: [PATCH 12/34] deleted a test that checked a capability that had not
 been implemented

---
 dymos/phase/test/test_timeseries.py             | 13 -------------
 dymos/transcriptions/common/timeseries_group.py |  0
 2 files changed, 13 deletions(-)
 create mode 100644 dymos/transcriptions/common/timeseries_group.py

diff --git a/dymos/phase/test/test_timeseries.py b/dymos/phase/test/test_timeseries.py
index f57578db8..c6835d4a3 100644
--- a/dymos/phase/test/test_timeseries.py
+++ b/dymos/phase/test/test_timeseries.py
@@ -629,19 +629,6 @@ def test_timeseries_expr_gl(self):
         R_ts = p.get_val('traj.phase0.timeseries.R')
         assert_near_equal(R_computed, R_ts, tolerance=1e-12)
 
-    def test_timeseries_expr_explicit_shooting(self):
-        p = min_time_climb(transcription_class=dm.ExplicitShooting, num_seg=12,
-                           transcription_order=3, timeseries_expr=True)
-        p.run_model()
-
-        pres = p.get_val('traj.phase0.timeseries.pres')
-        temp = p.get_val('traj.phase0.timeseries.temp')
-        rho = p.get_val('traj.phase0.timeseries.rho')
-
-        R_computed = pres / (temp * rho)
-        R_ts = p.get_val('traj.phase0.timeseries.R')
-        assert_near_equal(R_computed, R_ts, tolerance=1e-12)
-
 
 if __name__ == '__main__':  # pragma: no cover
     unittest.main()
diff --git a/dymos/transcriptions/common/timeseries_group.py b/dymos/transcriptions/common/timeseries_group.py
new file mode 100644
index 000000000..e69de29bb

From 1a080cd80466c975f7486935a842eace8c7debdb Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Wed, 26 Oct 2022 14:03:34 -0500
Subject: [PATCH 13/34] Added a timeseries group and moved timeseries output
 comp and the exec comp inside it

---
 dymos/phase/phase.py                          |  4 +--
 dymos/phase/test/test_timeseries.py           | 28 +++++++++----------
 dymos/transcriptions/analytic/analytic.py     | 14 ++++------
 dymos/transcriptions/common/__init__.py       |  1 +
 .../transcriptions/common/timeseries_group.py | 22 +++++++++++++++
 .../explicit_shooting/explicit_shooting.py    | 11 ++++----
 .../pseudospectral/pseudospectral_base.py     | 14 ++++------
 dymos/transcriptions/solve_ivp/solve_ivp.py   | 25 +++++++++--------
 dymos/transcriptions/transcription_base.py    |  2 +-
 dymos/utils/introspection.py                  | 12 ++++----
 10 files changed, 77 insertions(+), 56 deletions(-)

diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index f0a5c88fd..8ef617e14 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -2071,13 +2071,13 @@ def initialize_values_from_phase(self, prob, from_phase, phase_path='', skip_par
             op_dict = MPI.COMM_WORLD.bcast(op_dict, root=0)
 
         # Set the integration times
-        op = op_dict['timeseries.time']
+        op = op_dict['timeseries.timeseries_comp.time']
         prob.set_val(f'{self_path}t_initial', op['val'][0, ...])
         prob.set_val(f'{self_path}t_duration', op['val'][-1, ...] - op['val'][0, ...])
 
         # Assign initial state values
         for name in phs.state_options:
-            op = op_dict[f'timeseries.states:{name}']
+            op = op_dict[f'timeseries.timeseries_comp.states:{name}']
             prob[f'{self_path}initial_states:{name}'][...] = op['val'][0, ...]
 
         # Assign control values
diff --git a/dymos/phase/test/test_timeseries.py b/dymos/phase/test/test_timeseries.py
index c6835d4a3..e11343fec 100644
--- a/dymos/phase/test/test_timeseries.py
+++ b/dymos/phase/test/test_timeseries.py
@@ -10,6 +10,7 @@
 
 import dymos as dm
 from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE
+from dymos.examples.min_time_climb.min_time_climb_ode import MinTimeClimbODE
 
 from dymos.models.atmosphere import USatm1976Comp
 from dymos.examples.min_time_climb.aero import AeroGroup
@@ -288,8 +289,8 @@ def setup(self):
                            subsys=FlightPathEOM2D(num_nodes=nn),
                            promotes_inputs=['m', 'v', 'gam', 'alpha'])
 
-        # foo = self.add_subsystem('foo', om.IndepVarComp())
-        # foo.add_output('rho', val=100 * np.ones(nn), units='g/cm**3')
+        foo = self.add_subsystem('foo', om.IndepVarComp())
+        foo.add_output('rho', val=100 * np.ones(nn), units='g/cm**3')
 
         self.connect('aero.f_drag', 'flight_dynamics.D')
         self.connect('aero.f_lift', 'flight_dynamics.L')
@@ -308,7 +309,9 @@ def min_time_climb(num_seg=3, transcription_class=dm.Radau, transcription_order=
 
     traj = dm.Trajectory()
 
-    phase = dm.Phase(ode_class=MinTimeClimbODEDuplicateOutput, transcription=tx)
+    ode_class = MinTimeClimbODEDuplicateOutput if not timeseries_expr else MinTimeClimbODE
+
+    phase = dm.Phase(ode_class=ode_class, transcription=tx)
     traj.add_phase('phase0', phase)
 
     p.model.add_subsystem('traj', traj)
@@ -432,7 +435,7 @@ def make_problem_brachistochrone(transcription, polynomial_control=False):
         phase.add_boundary_constraint('x', loc='final', equals=10.0)
         phase.add_boundary_constraint('y', loc='final', equals=5.0)
 
-        phase.add_parameter('g', opt=True, units='m/s**2', val=9.80665, include_timeseries=True)
+        phase.add_parameter('g', opt=False, units='m/s**2', val=9.80665, include_timeseries=True)
 
         phase.add_objective('time_phase', loc='final', scaler=10)
         phase.add_timeseries_output('z=x*y + x**2', units='m**2')
@@ -568,19 +571,16 @@ def test_timeseries_expr_solve_ivp(self):
         p = self.make_problem_brachistochrone(transcription=tx)
         p.run_driver()
 
-        t = p.get_val('phase0.timeseries.time')
         z_ts = p.get_val('phase0.timeseries.z')
         f_ts = p.get_val('phase0.timeseries.f')
 
         phase0 = p.model._get_subsystem('phase0')
-        sim = phase0.simulate()
+        sim = phase0.simulate(times_per_seg=30)
 
         z_sim = sim.get_val('phase0.timeseries.z')
         f_sim = sim.get_val('phase0.timeseries.f')
-        t_sim = sim.get_val('phase0.timeseries.time')
-
-        assert_timeseries_near_equal(t, z_ts, t_sim, z_sim, abs_tolerance=1e-1)
-        assert_timeseries_near_equal(t, f_ts, t_sim, f_sim, abs_tolerance=1e-1)
+        assert_near_equal(z_ts[-1], z_sim[-1], tolerance=1e-3)
+        assert_near_equal(f_ts[-1], f_sim[-1], tolerance=1e-3)
 
     def test_timeseries_expr_solve_ivp_polynomial_controls(self):
         tx = dm.Radau(num_segments=5, order=3, compressed=True)
@@ -588,18 +588,16 @@ def test_timeseries_expr_solve_ivp_polynomial_controls(self):
         p = self.make_problem_brachistochrone(transcription=tx, polynomial_control=True)
         p.run_driver()
 
-        t = p.get_val('phase0.timeseries.time')
         z_ts = p.get_val('phase0.timeseries.z')
         f_ts = p.get_val('phase0.timeseries.f')
 
         phase0 = p.model._get_subsystem('phase0')
-        sim = phase0.simulate()
+        sim = phase0.simulate(times_per_seg=30)
 
         z_sim = sim.get_val('phase0.timeseries.z')
         f_sim = sim.get_val('phase0.timeseries.f')
-        t_sim = sim.get_val('phase0.timeseries.time')
-        assert_timeseries_near_equal(t, z_ts, t_sim, z_sim, abs_tolerance=1e-3)
-        assert_timeseries_near_equal(t, f_ts, t_sim, f_sim, abs_tolerance=1e-3)
+        assert_near_equal(z_ts[-1], z_sim[-1], tolerance=1e-3)
+        assert_near_equal(f_ts[-1], f_sim[-1], tolerance=1e-3)
 
 
 class TestTimeseriesMinTimeClimb(unittest.TestCase):
diff --git a/dymos/transcriptions/analytic/analytic.py b/dymos/transcriptions/analytic/analytic.py
index 38b8a5da6..8d2d049b1 100644
--- a/dymos/transcriptions/analytic/analytic.py
+++ b/dymos/transcriptions/analytic/analytic.py
@@ -8,7 +8,7 @@
 from ...utils.indexing import get_src_indices_by_row
 from ..grid_data import GridData
 from .analytic_timeseries_output_comp import AnalyticTimeseriesOutputComp
-from ..common.time_comp import TimeComp
+from ..common import TimeComp, TimeseriesOutputGroup
 
 
 class Analytic(TranscriptionBase):
@@ -266,13 +266,10 @@ def setup_timeseries_outputs(self, phase):
         gd = self.grid_data
 
         for name, options in phase._timeseries.items():
-            expr_ts = False
+            has_expr = False
             for _, output_options in options['outputs'].items():
                 if output_options['is_expr']:
-                    ts_exec_comp = om.ExecComp(has_diag_partials=True)
-                    phase.add_subsystem('timeseries_exec_comp', ts_exec_comp)
-                    expr_ts = True
-                if expr_ts:
+                    has_expr = True
                     break
             if options['transcription'] is None:
                 ogd = None
@@ -284,9 +281,10 @@ def setup_timeseries_outputs(self, phase):
                                                            output_subset=options['subset'],
                                                            time_units=phase.time_options['units'])
 
-            phase.add_subsystem(name, subsys=timeseries_comp)
+            timeseries_group = TimeseriesOutputGroup(has_expr=has_expr, timeseries_output_comp=timeseries_comp)
+            phase.add_subsystem(name, subsys=timeseries_group)
 
-            phase.connect('dt_dstau', (f'{name}.dt_dstau'), flat_src_indices=True)
+            phase.connect('dt_dstau', f'{name}.dt_dstau', flat_src_indices=True)
 
     def configure_defects(self, phase):
         """
diff --git a/dymos/transcriptions/common/__init__.py b/dymos/transcriptions/common/__init__.py
index aac6319aa..edbf127fb 100644
--- a/dymos/transcriptions/common/__init__.py
+++ b/dymos/transcriptions/common/__init__.py
@@ -3,3 +3,4 @@
 from .parameter_comp import ParameterComp
 from .polynomial_control_group import PolynomialControlGroup
 from .time_comp import TimeComp
+from .timeseries_group import TimeseriesOutputGroup
diff --git a/dymos/transcriptions/common/timeseries_group.py b/dymos/transcriptions/common/timeseries_group.py
index e69de29bb..daf148aa4 100644
--- a/dymos/transcriptions/common/timeseries_group.py
+++ b/dymos/transcriptions/common/timeseries_group.py
@@ -0,0 +1,22 @@
+import openmdao.api as om
+import numpy as np
+
+from .timeseries_output_comp import TimeseriesOutputCompBase
+
+
+class TimeseriesOutputGroup(om.Group):
+    def initialize(self):
+        self.options.declare('timeseries_output_comp', types=TimeseriesOutputCompBase,
+                             desc='Timeseries component specific to the transcription of the optimal control problem.')
+
+        self.options.declare('has_expr', types=bool,
+                             desc='If true, timeseries group has an expression to be computed')
+
+    def setup(self):
+        timeseries_output_comp = self.options['timeseries_output_comp']
+        has_expr = self.options['has_expr']
+        if has_expr:
+            self.add_subsystem('timeseries_exec_comp', om.ExecComp(has_diag_partials=True))
+
+        self.add_subsystem('timeseries_comp', timeseries_output_comp,
+                           promotes_inputs=['*'], promotes_outputs=['*'])
diff --git a/dymos/transcriptions/explicit_shooting/explicit_shooting.py b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
index 854844296..ac5e00097 100644
--- a/dymos/transcriptions/explicit_shooting/explicit_shooting.py
+++ b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
@@ -10,6 +10,7 @@
 from ...utils.misc import get_rate_units, CoerceDesvar
 from ...utils.introspection import get_promoted_vars, get_source_metadata
 from ...utils.constants import INF_BOUND
+from ..common import TimeseriesOutputGroup
 
 
 class ExplicitShooting(TranscriptionBase):
@@ -396,17 +397,15 @@ def setup_timeseries_outputs(self, phase):
         """
         gd = self.grid_data
         for name, options in phase._timeseries.items():
-            expr_ts = False
+            has_expr = False
             for _, output_options in options['outputs'].items():
                 if output_options['is_expr']:
-                    ts_exec_comp = om.ExecComp(has_diag_partials=True)
-                    phase.add_subsystem('timeseries_exec_comp', ts_exec_comp)
-                    expr_ts = True
-                if expr_ts:
+                    has_expr = True
                     break
             timeseries_comp = ExplicitTimeseriesComp(input_grid_data=gd,
                                                      output_subset='segment_ends')
-            phase.add_subsystem(name, subsys=timeseries_comp)
+            timeseries_group = TimeseriesOutputGroup(has_expr=has_expr, timeseries_output_comp=timeseries_comp)
+            phase.add_subsystem(name, subsys=timeseries_group)
 
     def configure_timeseries_outputs(self, phase):
         """
diff --git a/dymos/transcriptions/pseudospectral/pseudospectral_base.py b/dymos/transcriptions/pseudospectral/pseudospectral_base.py
index 744966b0d..c2ea41c09 100644
--- a/dymos/transcriptions/pseudospectral/pseudospectral_base.py
+++ b/dymos/transcriptions/pseudospectral/pseudospectral_base.py
@@ -2,7 +2,7 @@
 
 import openmdao.api as om
 from ..transcription_base import TranscriptionBase
-from ..common import TimeComp
+from ..common import TimeComp, TimeseriesOutputGroup
 from .components import StateIndependentsComp, StateInterpComp, CollocationComp, \
     PseudospectralTimeseriesOutputComp
 from ...utils.misc import CoerceDesvar, get_rate_units, reshape_val
@@ -483,13 +483,10 @@ def setup_timeseries_outputs(self, phase):
         gd = self.grid_data
 
         for name, options in phase._timeseries.items():
-            expr_ts = False
+            has_expr = False
             for _, output_options in options['outputs'].items():
                 if output_options['is_expr']:
-                    ts_exec_comp = om.ExecComp(has_diag_partials=True)
-                    phase.add_subsystem('timeseries_exec_comp', ts_exec_comp)
-                    expr_ts = True
-                if expr_ts:
+                    has_expr = True
                     break
 
             if options['transcription'] is None:
@@ -501,9 +498,10 @@ def setup_timeseries_outputs(self, phase):
                                                                  output_grid_data=ogd,
                                                                  output_subset=options['subset'],
                                                                  time_units=phase.time_options['units'])
-            phase.add_subsystem(name, subsys=timeseries_comp)
+            timeseries_group = TimeseriesOutputGroup(has_expr=has_expr, timeseries_output_comp=timeseries_comp)
+            phase.add_subsystem(name, subsys=timeseries_group)
 
-            phase.connect('dt_dstau', (f'{name}.dt_dstau'), flat_src_indices=True)
+            phase.connect('dt_dstau', f'{name}.dt_dstau', flat_src_indices=True)
 
     def _get_objective_src(self, var, loc, phase, ode_outputs=None):
         """
diff --git a/dymos/transcriptions/solve_ivp/solve_ivp.py b/dymos/transcriptions/solve_ivp/solve_ivp.py
index 7cca34916..95aacced4 100644
--- a/dymos/transcriptions/solve_ivp/solve_ivp.py
+++ b/dymos/transcriptions/solve_ivp/solve_ivp.py
@@ -8,7 +8,7 @@
 from ..transcription_base import TranscriptionBase
 from .components import SegmentSimulationComp, SegmentStateMuxComp, \
     SolveIVPControlGroup, SolveIVPPolynomialControlGroup, SolveIVPTimeseriesOutputComp
-from ..common import TimeComp
+from ..common import TimeComp, TimeseriesOutputGroup
 from ...utils.misc import get_rate_units
 from ...utils.introspection import get_promoted_vars, get_targets, get_source_metadata, get_target_metadata
 from ...utils.indexing import get_src_indices_by_row
@@ -503,23 +503,26 @@ def setup_timeseries_outputs(self, phase):
             The phase object to which this transcription instance applies.
         """
         gd = self.grid_data
-        # Remove all timeseries other than 'timeseries'
-        for ts_name, options in phase._timeseries.items():
-            if ts_name != 'timeseries':
-                phase._timeseries.pop(ts_name)
-                continue
-            ts_exec_comp = om.ExecComp(has_diag_partials=True)
-            phase.add_subsystem('timeseries_exec_comp', ts_exec_comp)
-            expr_ts = True
-            if expr_ts:
+        # Check if timeseries contains an expression that needs to be evaluated
+        for _, output_options in phase._timeseries['timeseries']['outputs'].items():
+            if output_options['is_expr']:
+                has_expr = True
                 break
+            else:
+                has_expr = False
 
         timeseries_comp = \
             SolveIVPTimeseriesOutputComp(input_grid_data=gd,
                                          output_nodes_per_seg=self.options['output_nodes_per_seg'],
                                          time_units=phase.time_options['units'])
 
-        phase.add_subsystem('timeseries', subsys=timeseries_comp)
+        timeseries_group = TimeseriesOutputGroup(has_expr=has_expr, timeseries_output_comp=timeseries_comp)
+        phase.add_subsystem('timeseries', subsys=timeseries_group)
+
+        # Remove all subsequent timeseries
+        for ts_name in list(phase._timeseries.keys()):
+            if ts_name != 'timeseries':
+                phase._timeseries.pop(ts_name)
 
     def get_parameter_connections(self, name, phase):
         """
diff --git a/dymos/transcriptions/transcription_base.py b/dymos/transcriptions/transcription_base.py
index 4f1be9ed5..e994afb22 100644
--- a/dymos/transcriptions/transcription_base.py
+++ b/dymos/transcriptions/transcription_base.py
@@ -381,7 +381,7 @@ def configure_timeseries_outputs(self, 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)
+            timeseries_comp = phase._get_subsystem(f'{timeseries_name}.timeseries_comp')
 
             for ts_output_name, ts_output in timeseries_options['outputs'].items():
                 name = ts_output['output_name'] if ts_output['output_name'] is not None else ts_output['name']
diff --git a/dymos/utils/introspection.py b/dymos/utils/introspection.py
index 3ce2d79aa..4a69c85a3 100644
--- a/dymos/utils/introspection.py
+++ b/dymos/utils/introspection.py
@@ -735,7 +735,7 @@ def configure_timeseries_output_introspection(phase):
 
         for output_name, output_options in ts_opts['outputs'].items():
             if output_options['is_expr']:
-                output_meta = configure_timeseries_expr_introspection(phase, output_options)
+                output_meta = configure_timeseries_expr_introspection(phase, ts_name, output_options)
 
             else:
                 try:
@@ -763,7 +763,7 @@ def configure_timeseries_output_introspection(phase):
             ts_opts['outputs'].pop(s)
 
 
-def configure_timeseries_expr_introspection(phase, expr_options):
+def configure_timeseries_expr_introspection(phase, timeseries_name, expr_options):
     """
     Introspect the timeseries expressions, add expressions to the exec comp, and make connections to the input variables
 
@@ -771,11 +771,13 @@ def configure_timeseries_expr_introspection(phase, expr_options):
     ----------
     phase : Phase
         The phase object whose timeseries outputs are to be introspected.
+    timeseries_name : str
+        The name of the timeseries that the expr is to be added to
     expr_options : dict
         Options for the expression to be added to the timeseries exec comp
     """
 
-    timeseries_ec = phase._get_subsystem('timeseries_exec_comp')
+    timeseries_ec = phase._get_subsystem(f'{timeseries_name}.timeseries_exec_comp')
     transcription = phase.options['transcription']
     num_output_nodes = transcription._get_num_timeseries_nodes()
 
@@ -807,7 +809,7 @@ def configure_timeseries_expr_introspection(phase, expr_options):
         elif var == expr_lhs:
             var_units = units
             var_shape = shape
-            meta['src'] = f'timeseries_exec_comp.{var}'
+            meta['src'] = f'{timeseries_name}.timeseries_exec_comp.{var}'
             phase.timeseries_ec_vars[var]['added_to_ec'] = True
         else:
             try:
@@ -832,7 +834,7 @@ def configure_timeseries_expr_introspection(phase, expr_options):
         if phase.timeseries_ec_vars[var]['added_to_ec']:
             continue
         else:
-            phase.connect(src_name=options['src'], tgt_name=f'timeseries_exec_comp.{var}',
+            phase.connect(src_name=options['src'], tgt_name=f'{timeseries_name}.timeseries_exec_comp.{var}',
                           src_indices=options['src_idxs'])
 
         phase.timeseries_ec_vars[var]['added_to_ec'] = True

From 9b0c1b0e30b6a73446bc039d0f88f0a4cecf0f90 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Wed, 26 Oct 2022 14:20:39 -0500
Subject: [PATCH 14/34] Changed timeseries plotting to go one level deeper to
 extract timeseries metadata

---
 dymos/visualization/timeseries_plots.py | 18 ++++++++++--------
 1 file changed, 10 insertions(+), 8 deletions(-)

diff --git a/dymos/visualization/timeseries_plots.py b/dymos/visualization/timeseries_plots.py
index 19ed7bcf3..1869d9d4b 100644
--- a/dymos/visualization/timeseries_plots.py
+++ b/dymos/visualization/timeseries_plots.py
@@ -396,17 +396,19 @@ def timeseries_plots(solution_recorder_filename, simulation_record_file=None, pl
                 break
         # get the time units first so they can be associated with all the variables
         for timeseries_node_child in timeseries_node['children']:
-            if timeseries_node_child['name'] == 'time':
-                units_for_time = timeseries_node_child['units']
-                break
+            for timeseries_node_g_child in timeseries_node_child['children']:
+                if timeseries_node_g_child['name'] == 'time':
+                    units_for_time = timeseries_node_g_child['units']
+                    break
         for timeseries_node_child in timeseries_node['children']:
             # plot everything in the timeseries except input_values. Also not time since that
             #   is the independent variable in these plot
-            if not timeseries_node_child['name'].startswith("input_values:")\
-                    and not timeseries_node_child['name'] == 'time':
-                varname = timeseries_node_child['name']
-                var_units[varname] = timeseries_node_child['units']
-                time_units[varname] = units_for_time
+            for timeseries_node_g_child in timeseries_node_child['children']:
+                if not timeseries_node_g_child['name'].startswith("input_values:")\
+                        and not timeseries_node_g_child['name'] == 'time':
+                    varname = timeseries_node_g_child['name']
+                    var_units[varname] = timeseries_node_g_child['units']
+                    time_units[varname] = units_for_time
 
     # Check to see if there is anything to plot
     if len(var_units) == 0:

From acea66f193bd98fd994f5991a7c2fd63ccb71dd0 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Wed, 26 Oct 2022 16:25:46 -0500
Subject: [PATCH 15/34] Updated name of the timesereies comp subsystem to match
 new name

---
 ...test_brachistochrone_control_rate_targets.py |  1 +
 dymos/phase/test/test_timeseries.py             |  2 +-
 dymos/transcriptions/common/timeseries_group.py | 17 +++++++++++++++++
 .../explicit_shooting/explicit_shooting.py      |  2 +-
 dymos/utils/introspection.py                    | 13 ++++++++++---
 5 files changed, 30 insertions(+), 5 deletions(-)

diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py b/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py
index 2edd1e06e..2265dbab8 100644
--- a/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py
+++ b/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py
@@ -193,6 +193,7 @@ def test_brachistochrone_control_rate_targets_gauss_lobatto(self):
         theta_exp = exp_out.get_val('phase0.timeseries.control_rates:theta_rate')
 
         io_meta = p.model.phase0.timeseries.get_io_metadata(iotypes=('input', 'output'), get_remote=True)
+        print(io_meta)
         self.assertEqual(io_meta['controls:theta']['units'], 'rad*s')
 
         ax.plot(t_imp, theta_imp, 'ro', label='solution')
diff --git a/dymos/phase/test/test_timeseries.py b/dymos/phase/test/test_timeseries.py
index e11343fec..91eb7da20 100644
--- a/dymos/phase/test/test_timeseries.py
+++ b/dymos/phase/test/test_timeseries.py
@@ -16,7 +16,7 @@
 from dymos.examples.min_time_climb.aero import AeroGroup
 from dymos.examples.min_time_climb.prop import PropGroup
 from dymos.models.eom import FlightPathEOM2D
-from dymos.utils.testing_utils import assert_timeseries_near_equal
+
 
 @use_tempdirs
 class TestTimeseriesOutput(unittest.TestCase):
diff --git a/dymos/transcriptions/common/timeseries_group.py b/dymos/transcriptions/common/timeseries_group.py
index daf148aa4..b5374cb0c 100644
--- a/dymos/transcriptions/common/timeseries_group.py
+++ b/dymos/transcriptions/common/timeseries_group.py
@@ -5,7 +5,21 @@
 
 
 class TimeseriesOutputGroup(om.Group):
+    """
+    Class definition of the TimeComp.
+
+    Parameters
+    ----------
+    **kwargs : dict
+        Dictionary of optional arguments.
+    """
+    def __init__(self, **kwargs):
+        super().__init__(**kwargs)
+
     def initialize(self):
+        """
+        Declare component options.
+        """
         self.options.declare('timeseries_output_comp', types=TimeseriesOutputCompBase,
                              desc='Timeseries component specific to the transcription of the optimal control problem.')
 
@@ -13,6 +27,9 @@ def initialize(self):
                              desc='If true, timeseries group has an expression to be computed')
 
     def setup(self):
+        """
+        Define the structure of the timeseries group.
+        """
         timeseries_output_comp = self.options['timeseries_output_comp']
         has_expr = self.options['has_expr']
         if has_expr:
diff --git a/dymos/transcriptions/explicit_shooting/explicit_shooting.py b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
index ac5e00097..9e99ca3a6 100644
--- a/dymos/transcriptions/explicit_shooting/explicit_shooting.py
+++ b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
@@ -419,7 +419,7 @@ def configure_timeseries_outputs(self, phase):
         integrator_comp = phase._get_subsystem('integrator')
         integrator_comp._configure_timeseries_outputs()
         for timeseries_name, timeseries_options in phase._timeseries.items():
-            timeseries_comp = phase._get_subsystem(timeseries_name)
+            timeseries_comp = phase._get_subsystem(f'{timeseries_name}.timeseries_comp')
 
             for output_name, options in integrator_comp._filtered_timeseries_outputs.items():
                 added_src = timeseries_comp._add_output_configure(output_name,
diff --git a/dymos/utils/introspection.py b/dymos/utils/introspection.py
index 4a69c85a3..f410acbfc 100644
--- a/dymos/utils/introspection.py
+++ b/dymos/utils/introspection.py
@@ -32,6 +32,8 @@ def classify_var(var, state_options, parameter_options, control_options,
         For each control variable, a dictionary of its options, keyed by name.
     polynomial_control_options : dict of {str: OptionsDictionary}
         For each polynomial variable, a dictionary of its options, keyed by name.
+    timeseries_options : {str: OptionsDictionary}
+        For each timeseries, a dictionary of its options, keyed by name.
 
     Returns
     -------
@@ -765,16 +767,21 @@ def configure_timeseries_output_introspection(phase):
 
 def configure_timeseries_expr_introspection(phase, timeseries_name, expr_options):
     """
-    Introspect the timeseries expressions, add expressions to the exec comp, and make connections to the input variables
+    Introspect the timeseries expressions, add expressions to the exec comp, and make connections to the variables.
 
     Parameters
     ----------
     phase : Phase
         The phase object whose timeseries outputs are to be introspected.
     timeseries_name : str
-        The name of the timeseries that the expr is to be added to
+        The name of the timeseries that the expr is to be added to.
     expr_options : dict
-        Options for the expression to be added to the timeseries exec comp
+        Options for the expression to be added to the timeseries exec comp.
+
+    Returns
+    -------
+    dict
+        A dictionary containing the metadata for the source. This consists of shape, units, and tags.
     """
 
     timeseries_ec = phase._get_subsystem(f'{timeseries_name}.timeseries_exec_comp')

From c78f1366e3925a4c3e08517b8d460b9fdd711717 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Thu, 27 Oct 2022 09:40:52 -0500
Subject: [PATCH 16/34] Added recordable=False to the option to supress the
 warnings

---
 dymos/transcriptions/common/timeseries_group.py | 3 +--
 1 file changed, 1 insertion(+), 2 deletions(-)

diff --git a/dymos/transcriptions/common/timeseries_group.py b/dymos/transcriptions/common/timeseries_group.py
index b5374cb0c..fcaa06004 100644
--- a/dymos/transcriptions/common/timeseries_group.py
+++ b/dymos/transcriptions/common/timeseries_group.py
@@ -1,5 +1,4 @@
 import openmdao.api as om
-import numpy as np
 
 from .timeseries_output_comp import TimeseriesOutputCompBase
 
@@ -20,7 +19,7 @@ def initialize(self):
         """
         Declare component options.
         """
-        self.options.declare('timeseries_output_comp', types=TimeseriesOutputCompBase,
+        self.options.declare('timeseries_output_comp', types=TimeseriesOutputCompBase, recordable=False,
                              desc='Timeseries component specific to the transcription of the optimal control problem.')
 
         self.options.declare('has_expr', types=bool,

From 20471785888c485f97a507e5db72eb53fa429d8a Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Thu, 27 Oct 2022 13:25:16 -0500
Subject: [PATCH 17/34] Modified test problems to match new API

---
 ...st_brachistochrone_control_rate_targets.py |  4 +--
 .../test_two_burn_orbit_raise_linkages.py     | 13 ++++----
 .../test/test_ex_min_time_climb.py            | 31 ++++++++++---------
 dymos/test/test_load_case.py                  | 13 ++++----
 .../pseudospectral/test/test_ps_base.py       |  3 +-
 dymos/utils/introspection.py                  |  3 +-
 6 files changed, 35 insertions(+), 32 deletions(-)

diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py b/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py
index 2265dbab8..86d294bc1 100644
--- a/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py
+++ b/dymos/examples/brachistochrone/test/test_brachistochrone_control_rate_targets.py
@@ -192,8 +192,8 @@ def test_brachistochrone_control_rate_targets_gauss_lobatto(self):
         t_exp = exp_out.get_val('phase0.timeseries.time')
         theta_exp = exp_out.get_val('phase0.timeseries.control_rates:theta_rate')
 
-        io_meta = p.model.phase0.timeseries.get_io_metadata(iotypes=('input', 'output'), get_remote=True)
-        print(io_meta)
+        io_meta = p.model.phase0.timeseries.timeseries_comp.get_io_metadata(
+            iotypes=('input', 'output'), get_remote=True)
         self.assertEqual(io_meta['controls:theta']['units'], 'rad*s')
 
         ax.plot(t_imp, theta_imp, 'ro', label='solution')
diff --git a/dymos/examples/finite_burn_orbit_raise/test/test_two_burn_orbit_raise_linkages.py b/dymos/examples/finite_burn_orbit_raise/test/test_two_burn_orbit_raise_linkages.py
index c0b6a3948..1a51ea977 100644
--- a/dymos/examples/finite_burn_orbit_raise/test/test_two_burn_orbit_raise_linkages.py
+++ b/dymos/examples/finite_burn_orbit_raise/test/test_two_burn_orbit_raise_linkages.py
@@ -2,6 +2,7 @@
 import warnings
 
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
+from dymos.utils.introspection import get_promoted_vars
 
 import matplotlib.pyplot as plt
 plt.switch_backend('Agg')
@@ -676,18 +677,16 @@ def test_two_burn_orbit_raise_gl_list_add_timeseries_output(self):
                           tolerance=2.0E-3)
 
         # check units
-        un = dict(p.model.list_outputs(units=True))
-        assert(un['traj.phases.burn1.timeseries.pos_x']['units'] == 'm')
-        assert(un['traj.phases.burn1.timeseries.pos_y']['units'] == 'DU')
-        assert(un['traj.phases.burn2.timeseries.pos_x']['units'] == 'm')
-        assert(un['traj.phases.burn2.timeseries.pos_y']['units'] == 'm')
+        un = get_promoted_vars(p.model, iotypes=('input', 'output'))
+        assert(un['traj.burn1.timeseries.pos_x']['units'] == 'm')
+        assert(un['traj.burn1.timeseries.pos_y']['units'] == 'DU')
+        assert(un['traj.burn2.timeseries.pos_x']['units'] == 'm')
+        assert(un['traj.burn2.timeseries.pos_y']['units'] == 'm')
 
     @use_tempdirs
     def test_two_burn_orbit_raise_gl_wildcard_add_timeseries_output(self):
         import numpy as np
 
-        import matplotlib.pyplot as plt
-
         import openmdao.api as om
         from openmdao.utils.assert_utils import assert_near_equal
 
diff --git a/dymos/examples/min_time_climb/test/test_ex_min_time_climb.py b/dymos/examples/min_time_climb/test/test_ex_min_time_climb.py
index 63f223be3..ea6030670 100644
--- a/dymos/examples/min_time_climb/test/test_ex_min_time_climb.py
+++ b/dymos/examples/min_time_climb/test/test_ex_min_time_climb.py
@@ -7,6 +7,8 @@
 import openmdao.api as om
 from openmdao.utils.assert_utils import assert_near_equal
 from dymos.utils.testing_utils import assert_timeseries_near_equal
+from dymos.utils.introspection import get_promoted_vars
+
 import dymos as dm
 from dymos.examples.min_time_climb.min_time_climb_ode import MinTimeClimbODE
 
@@ -133,37 +135,37 @@ def _test_results(self, p):
 
     def _test_wilcard_outputs(self, p):
         """ Test that all wilcard outputs are provided. """
-        output_dict = dict(p.model.list_outputs(units=True, out_stream=None))
+        output_dict = get_promoted_vars(p.model, iotypes=('output',))
         ts = {k: v for k, v in output_dict.items() if 'timeseries.' in k}
         for c in ['mach', 'CD0', 'kappa', 'CLa', 'CL', 'CD', 'q', 'f_lift', 'f_drag', 'thrust', 'm_dot']:
             assert(any([True for t in ts if 'timeseries.' + c in t]))
 
     def _test_timeseries_units(self, p):
         """ Test that the units from the timeseries are correct. """
-        output_dict = dict(p.model.list_outputs(units=True, out_stream=None))
-        assert(output_dict['traj.phases.phase0.timeseries.thrust']['units'] == 'lbf')  # no wildcard, from units dict
-        assert(output_dict['traj.phases.phase0.timeseries.m_dot']['units'] == 'kg/s')  # no wildcard, from ODE
-        assert(output_dict['traj.phases.phase0.timeseries.f_drag']['units'] == 'N')    # wildcard, from ODE
-        assert(output_dict['traj.phases.phase0.timeseries.f_lift']['units'] == 'lbf')  # wildcard, from units dict
+        output_dict = get_promoted_vars(p.model, iotypes=('output',))
+        assert(output_dict['traj.phase0.timeseries.thrust']['units'] == 'lbf')  # no wildcard, from units dict
+        assert(output_dict['traj.phase0.timeseries.m_dot']['units'] == 'kg/s')  # no wildcard, from ODE
+        assert(output_dict['traj.phase0.timeseries.f_drag']['units'] == 'N')    # wildcard, from ODE
+        assert(output_dict['traj.phase0.timeseries.f_lift']['units'] == 'lbf')  # wildcard, from units dict
 
     def _test_mach_rate(self, p, plot=False):
         """ Test that the mach rate is provided by the timeseries and is accurate. """
         # Verify correct timeseries output of mach_rate
-        output_dict = dict(p.model.list_outputs(units=True, out_stream=None))
+        output_dict = get_promoted_vars(p.model, iotypes=('output',))
         ts = {k: v for k, v in output_dict.items() if 'timeseries.' in k}
-        self.assertTrue('traj.phases.phase0.timeseries.mach_rate' in ts)
+        self.assertTrue('traj.phase0.timeseries.mach_rate' in ts)
 
         case = om.CaseReader('dymos_solution.db').get_case('final')
 
-        time = case['traj.phases.phase0.timeseries.time'][:, 0]
-        mach = case['traj.phases.phase0.timeseries.mach'][:, 0]
-        mach_rate = case['traj.phases.phase0.timeseries.mach_rate'][:, 0]
+        time = case['traj.phase0.timeseries.time'][:, 0]
+        mach = case['traj.phase0.timeseries.mach'][:, 0]
+        mach_rate = case['traj.phase0.timeseries.mach_rate'][:, 0]
 
         sim_case = om.CaseReader('dymos_simulation.db').get_case('final')
 
-        sim_time = sim_case['traj.phases.phase0.timeseries.time'][:, 0]
-        sim_mach = sim_case['traj.phases.phase0.timeseries.mach'][:, 0]
-        sim_mach_rate = sim_case['traj.phases.phase0.timeseries.mach_rate'][:, 0]
+        sim_time = sim_case['traj.phase0.timeseries.time'][:, 0]
+        sim_mach = sim_case['traj.phase0.timeseries.mach'][:, 0]
+        sim_mach_rate = sim_case['traj.phase0.timeseries.mach_rate'][:, 0]
 
         # Fit a numpy polynomial segment by segment to mach vs time, and compare the derivatives to mach_rate
         gd = p.model.traj.phases.phase0.options['transcription'].grid_data
@@ -249,5 +251,6 @@ def test_results_radau(self):
 
         self._test_mach_rate(p, plot=False)
 
+
 if __name__ == '__main__':  # pragma: no cover
     unittest.main()
diff --git a/dymos/test/test_load_case.py b/dymos/test/test_load_case.py
index c86dd3b45..0c3725b96 100644
--- a/dymos/test/test_load_case.py
+++ b/dymos/test/test_load_case.py
@@ -149,10 +149,11 @@ def test_load_case_lgl_to_radau(self):
         outputs = dict([(o[0], o[1]) for o in case.list_outputs(units=True, shape=True,
                                                                 out_stream=None)])
 
-        time_val = outputs['phase0.timeseries.time']['val']
-        theta_val = outputs['phase0.timeseries.controls:theta']['val']
+        print(outputs)
+        time_val = outputs['phase0.timeseries.timeseries_comp.time']['val']
+        theta_val = outputs['phase0.timeseries.timeseries_comp.controls:theta']['val']
 
-        assert_near_equal(q['phase0.timeseries.controls:theta'],
+        assert_near_equal(q['phase0.timeseries.timeseries_comp.controls:theta'],
                           q.model.phase0.interp(xs=time_val, ys=theta_val, nodes='all'),
                           tolerance=1.0E-3)
 
@@ -181,10 +182,10 @@ def test_load_case_radau_to_lgl(self):
         outputs = dict([(o[0], o[1]) for o in case.list_outputs(units=True, shape=True,
                                                                 out_stream=None)])
 
-        time_val = outputs['phase0.timeseries.time']['val']
-        theta_val = outputs['phase0.timeseries.controls:theta']['val']
+        time_val = outputs['phase0.timeseries.timeseries_comp.time']['val']
+        theta_val = outputs['phase0.timeseries.timeseries_comp.controls:theta']['val']
 
-        assert_near_equal(q['phase0.timeseries.controls:theta'],
+        assert_near_equal(q['phase0.timeseries.timeseries_comp.controls:theta'],
                           q.model.phase0.interp(xs=time_val, ys=theta_val, nodes='all'),
                           tolerance=1.0E-2)
 
diff --git a/dymos/transcriptions/pseudospectral/test/test_ps_base.py b/dymos/transcriptions/pseudospectral/test/test_ps_base.py
index 50e8d6cc0..94e8c012e 100644
--- a/dymos/transcriptions/pseudospectral/test/test_ps_base.py
+++ b/dymos/transcriptions/pseudospectral/test/test_ps_base.py
@@ -88,7 +88,8 @@ def test_add_state_units_none(self):
         p = make_problem(transcription=GaussLobatto, num_segments=10, order=3, compressed=True)
         p.run_model()
 
-        io_meta = p.model.traj.phases.phase.timeseries.get_io_metadata(iotypes=('output'), get_remote=True)
+        io_meta = p.model.traj.phases.phase.timeseries.timeseries_comp.get_io_metadata(
+            iotypes=('output',), get_remote=True)
         self.assertEqual(io_meta['states:x']['units'], None)
 
 
diff --git a/dymos/utils/introspection.py b/dymos/utils/introspection.py
index f410acbfc..dcfb07b94 100644
--- a/dymos/utils/introspection.py
+++ b/dymos/utils/introspection.py
@@ -144,8 +144,7 @@ def get_targets(ode, name, user_targets, control_rates=False):
     if isinstance(ode, dict):
         ode_inputs = ode
     else:
-        ode_inputs = {opts['prom_name']: opts for opts in
-                      ode.get_io_metadata(iotypes=('input',), get_remote=True).values()}
+        ode_inputs = get_promoted_vars(ode, iotypes=('input',))
     if user_targets is _unspecified:
         if name in ode_inputs and control_rates not in {1, 2}:
             return [name]

From 3e3baad2167dc0aef980dede9d6784a204c2ef55 Mon Sep 17 00:00:00 2001
From: Rob Falck <rfalck@nasa.gov>
Date: Fri, 28 Oct 2022 10:57:20 -0400
Subject: [PATCH 18/34] Added more tests for timeseries expressions and added
 the ability to pass keyword arguments to the enclosed ExecComp.

---
 dymos/phase/options.py              |   3 +
 dymos/phase/phase.py                |  20 ++--
 dymos/phase/test/test_timeseries.py | 154 ++++++++++++++++++++++++++++
 dymos/utils/introspection.py        |  33 +++++-
 4 files changed, 200 insertions(+), 10 deletions(-)

diff --git a/dymos/phase/options.py b/dymos/phase/options.py
index ee2220ed6..c6b021f5d 100644
--- a/dymos/phase/options.py
+++ b/dymos/phase/options.py
@@ -711,3 +711,6 @@ def __init__(self, read_only=False):
 
         self.declare(name='is_expr', default=False, allow_none=False,
                      desc='If true the requested timeseries is a mathematical expression')
+
+        self.declare(name='expr_kwargs', default={}, allow_none=False,
+                     desc='Options to be passed to the timeseries expression comp when adding the expression.')
diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index 8ef617e14..c5ab48938 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -1191,16 +1191,19 @@ def add_path_constraint(self, name, constraint_name=None, units=None, shape=None
                 self.add_timeseries_output(name, output_name=constraint_name, units=units, shape=shape)
 
     def add_timeseries_output(self, name, output_name=None, units=_unspecified, shape=_unspecified,
-                              timeseries='timeseries'):
+                              timeseries='timeseries', **kwargs):
         r"""
         Add a variable to the timeseries outputs of the phase.
 
+        If name is given as an expression, this expression will be passed to an OpenMDAO ExecComp and the result
+        computed and stored in the timeseries output under the variable name to the left of the equal sign.
+
         Parameters
         ----------
         name : str, or list of str
-            The name(s) of the variable to be used as a timeseries output.  Must be one of
-            'time', 'time_phase', one of the states, controls, control rates, or parameters,
-            in the phase, the path to an output variable in the ODE, or a glob pattern
+            The name(s) of the variable to be used as a timeseries output, or a mathematical expression to be used
+            as a timeseries output. If a name, it be one of 'time', 'time_phase', one of the states, controls,
+            control rates, or parameters, in the phase, the path to an output variable in the ODE, or a glob pattern
             matching some outputs in the ODE.
         output_name : str or None or list or dict
             The name of the variable as listed in the phase timeseries outputs.  By
@@ -1216,6 +1219,8 @@ def add_timeseries_output(self, name, output_name=None, units=_unspecified, shap
             since Dymos doesn't necessarily know the shape of ODE outputs until setup time.
         timeseries : str or None
             The name of the timeseries to which the output is being added.
+        **kwargs
+            Additional arguments passed to the exec
         """
         if type(name) is list:
             for i, name_i in enumerate(name):
@@ -1245,7 +1250,8 @@ def add_timeseries_output(self, name, output_name=None, units=_unspecified, shap
                                         shape=shape,
                                         timeseries=timeseries,
                                         rate=False,
-                                        expr=expr)
+                                        expr=expr,
+                                        expr_kwargs=kwargs)
 
     def add_timeseries_rate_output(self, name, output_name=None, units=_unspecified, shape=_unspecified,
                                    timeseries='timeseries'):
@@ -1303,7 +1309,7 @@ def add_timeseries_rate_output(self, name, output_name=None, units=_unspecified,
                                         rate=True)
 
     def _add_timeseries_output(self, name, output_name=None, units=_unspecified, shape=_unspecified,
-                               timeseries='timeseries', rate=False, expr=False):
+                               timeseries='timeseries', rate=False, expr=False, expr_kwargs=None):
         r"""
         Add a single variable or rate to the timeseries outputs of the phase.
 
@@ -1333,6 +1339,7 @@ def _add_timeseries_output(self, name, output_name=None, units=_unspecified, sha
         rate : bool
             If True, add the rate of change of the named variable to the timeseries outputs of the
             phase.  The rate variable will be named f'{name}_rate'.  Defaults to False.
+        expr :
 
         Returns
         -------
@@ -1364,6 +1371,7 @@ def _add_timeseries_output(self, name, output_name=None, units=_unspecified, sha
             ts_output['shape'] = shape
             ts_output['is_rate'] = rate
             ts_output['is_expr'] = expr
+            ts_output['expr_kwargs'] = expr_kwargs
 
             self._timeseries[timeseries]['outputs'][output_name] = ts_output
 
diff --git a/dymos/phase/test/test_timeseries.py b/dymos/phase/test/test_timeseries.py
index 91eb7da20..d49ae77ac 100644
--- a/dymos/phase/test/test_timeseries.py
+++ b/dymos/phase/test/test_timeseries.py
@@ -1,4 +1,5 @@
 import unittest
+import warnings
 
 import numpy as np
 
@@ -455,6 +456,159 @@ def make_problem_brachistochrone(transcription, polynomial_control=False):
         p['phase0.parameters:g'] = 9.80665
         return p
 
+    def test_invalid_expr_var(self):
+        p = om.Problem(model=om.Group())
+
+        tx = dm.Radau(num_segments=5, order=3, compressed=True)
+
+        phase = dm.Phase(ode_class=BrachistochroneODE, transcription=tx)
+
+        p.model.add_subsystem('phase0', phase)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+
+        phase.set_state_options('x', fix_initial=True)
+        phase.set_state_options('y', fix_initial=True)
+        phase.set_state_options('v', fix_initial=True)
+
+        phase.add_control('theta', continuity=True, rate_continuity=True, opt=True,
+                          units='deg', lower=0.01, upper=179.9, ref=1, ref0=0)
+
+        phase.add_parameter('g', opt=False, units='m/s**2', val=9.80665, include_timeseries=True)
+
+        phase.add_timeseries_output('k=units', units='m**2')
+
+        with self.assertRaises(RuntimeError) as e:
+            p.setup(check=True)
+
+        expected = "phase0: Timeseries expression `k=units` uses the following disallowed variable " \
+                   "name in its definition: 'units'."
+        self.assertEqual(str(e.exception.__cause__), expected)
+
+    def test_input_units(self):
+        p = om.Problem(model=om.Group())
+
+        tx = dm.Radau(num_segments=5, order=3, compressed=True)
+
+        phase = dm.Phase(ode_class=BrachistochroneODE, transcription=tx)
+
+        p.model.add_subsystem('phase0', phase)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+
+        phase.set_state_options('x', fix_initial=True)
+        phase.set_state_options('y', fix_initial=True)
+        phase.set_state_options('v', fix_initial=True)
+
+        phase.add_control('theta', continuity=True, rate_continuity=True, opt=True,
+                          units='deg', lower=0.01, upper=179.9, ref=1, ref0=0)
+
+        phase.add_parameter('g', opt=False, units='m/s**2', val=9.80665, include_timeseries=True)
+
+        phase.add_timeseries_output('sin_theta=sin(theta)', units='unitless', theta={'units': 'rad'})
+
+        p.setup(check=True)
+
+        p['phase0.t_initial'] = 0.0
+        p['phase0.t_duration'] = 2.0
+
+        p['phase0.states:x'] = phase.interp('x', [0, 10])
+        p['phase0.states:y'] = phase.interp('y', [10, 5])
+        p['phase0.states:v'] = phase.interp('v', [0, 9.9])
+        p[f'phase0.controls:theta'] = phase.interp('theta', [5, 100])
+        p['phase0.parameters:g'] = 9.80665
+
+        p.run_model()
+
+        sin_theta = p.get_val('phase0.timeseries.sin_theta')
+        theta = p.get_val('phase0.timeseries.controls:theta', units='rad')
+
+        assert_near_equal(np.sin(theta), sin_theta)
+
+    def test_output_units(self):
+        p = om.Problem(model=om.Group())
+
+        tx = dm.Radau(num_segments=5, order=3, compressed=True)
+
+        phase = dm.Phase(ode_class=BrachistochroneODE, transcription=tx)
+
+        p.model.add_subsystem('phase0', phase)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+
+        phase.set_state_options('x', fix_initial=True)
+        phase.set_state_options('y', fix_initial=True)
+        phase.set_state_options('v', fix_initial=True)
+
+        phase.add_control('theta', continuity=True, rate_continuity=True, opt=True,
+                          units='deg', lower=0.01, upper=179.9, ref=1, ref0=0)
+
+        phase.add_parameter('g', opt=False, units='m/s**2', val=9.80665, include_timeseries=True)
+
+        phase.add_timeseries_output('sin_theta=sin(theta)', units='unitless', theta={'units': 'rad'})
+        phase.add_timeseries_output('sin_theta2=sin(theta)', sin_theta2={'units': 'unitless'})
+        phase.add_timeseries_output('sin_theta3=sin(theta)', shape=(1,), sin_theta3={'shape': (1,)})
+
+        expected = 'phase0: User-provided shape for timeseries expression output `sin_theta3`using ' \
+                   'both the\n`shape` keyword argument and the `sin_theta3` keyword argument dictionary.\n' \
+                   'The latter will take precedence.'
+
+        with warnings.catch_warnings(record=True) as ctx:
+            warnings.simplefilter('always')
+            p.setup(check=True)
+
+        self.assertIn(expected, [str(warning.message) for warning in ctx])
+
+        p['phase0.t_initial'] = 0.0
+        p['phase0.t_duration'] = 2.0
+
+        p['phase0.states:x'] = phase.interp('x', [0, 10])
+        p['phase0.states:y'] = phase.interp('y', [10, 5])
+        p['phase0.states:v'] = phase.interp('v', [0, 9.9])
+        p[f'phase0.controls:theta'] = phase.interp('theta', [5, 100])
+        p['phase0.parameters:g'] = 9.80665
+
+        p.run_model()
+
+        from dymos.utils.introspection import get_promoted_vars
+
+        sin_theta_meta = get_promoted_vars(phase, iotypes='output', metadata_keys=('units',))['timeseries.sin_theta']
+        sin_theta2_meta = get_promoted_vars(phase, iotypes='output', metadata_keys=('units',))['timeseries.sin_theta2']
+
+        self.assertEqual('unitless', sin_theta_meta['units'])
+        self.assertEqual('unitless', sin_theta2_meta['units'])
+
+    def test_ignore_shape(self):
+        p = om.Problem(model=om.Group())
+
+        tx = dm.Radau(num_segments=5, order=3, compressed=True)
+
+        phase = dm.Phase(ode_class=BrachistochroneODE, transcription=tx)
+
+        p.model.add_subsystem('phase0', phase)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+
+        phase.set_state_options('x', fix_initial=True)
+        phase.set_state_options('y', fix_initial=True)
+        phase.set_state_options('v', fix_initial=True)
+
+        phase.add_control('theta', continuity=True, rate_continuity=True, opt=True,
+                          units='deg', lower=0.01, upper=179.9, ref=1, ref0=0)
+
+        phase.add_parameter('g', opt=False, units='m/s**2', val=9.80665, include_timeseries=True)
+
+        phase.add_timeseries_output('sin_theta=sin(theta)', units='unitless', theta={'units': 'rad', 'shape': (1,)})
+
+        expected = 'User-provided shape for timeseries expression input `theta` ignored.\n' \
+                   'Using automatically determined shape (20, 1).'
+
+        with warnings.catch_warnings(record=True) as ctx:
+            warnings.simplefilter('always')
+            p.setup(check=True)
+
+        self.assertIn(expected, [str(warning.message) for warning in ctx])
+
     def test_timeseries_expr_radau(self):
         tx = dm.Radau(num_segments=5, order=3, compressed=True)
         p = self.make_problem_brachistochrone(transcription=tx)
diff --git a/dymos/utils/introspection.py b/dymos/utils/introspection.py
index dcfb07b94..961e96aba 100644
--- a/dymos/utils/introspection.py
+++ b/dymos/utils/introspection.py
@@ -737,7 +737,6 @@ def configure_timeseries_output_introspection(phase):
         for output_name, output_options in ts_opts['outputs'].items():
             if output_options['is_expr']:
                 output_meta = configure_timeseries_expr_introspection(phase, ts_name, output_options)
-
             else:
                 try:
                     output_meta = transcription._get_timeseries_var_source(output_options['name'],
@@ -782,7 +781,6 @@ def configure_timeseries_expr_introspection(phase, timeseries_name, expr_options
     dict
         A dictionary containing the metadata for the source. This consists of shape, units, and tags.
     """
-
     timeseries_ec = phase._get_subsystem(f'{timeseries_name}.timeseries_exec_comp')
     transcription = phase.options['transcription']
     num_output_nodes = transcription._get_num_timeseries_nodes()
@@ -792,6 +790,7 @@ def configure_timeseries_expr_introspection(phase, timeseries_name, expr_options
     if '.' in expr_lhs or ':' in expr_lhs:
         raise ValueError(f'{expr_lhs} is an invalid name for the timeseries LHS. Must use a valid variable name')
     expr_reduced = expr
+
     units = expr_options['units'] if expr_options['units'] is not _unspecified else None
     shape = expr_options['shape'] if expr_options['shape'] is not _unspecified else (1,)
 
@@ -807,9 +806,14 @@ def configure_timeseries_expr_introspection(phase, timeseries_name, expr_options
         if var_name not in phase.timeseries_ec_vars.keys():
             phase.timeseries_ec_vars[var_name] = {'added_to_ec': False, 'abs_name': name}
     expr_vars = {}
-    expr_kwargs = {}
+    expr_kwargs = expr_options['expr_kwargs']
     meta = {'units': units, 'shape': shape}
+
     for var in phase.timeseries_ec_vars:
+        if var in ('output_name', 'units', 'shape', 'timeseries'):
+            raise RuntimeError(f'{phase.pathname}: Timeseries expression `{expr}` uses the following disallowed '
+                               f'variable name in its definition: \'{var}\'.')
+
         if phase.timeseries_ec_vars[var]['added_to_ec']:
             continue
         elif var == expr_lhs:
@@ -827,7 +831,28 @@ def configure_timeseries_expr_introspection(phase, timeseries_name, expr_options
             var_units = expr_vars[var]['units']
             var_shape = expr_vars[var]['shape']
 
-        expr_kwargs[var] = {'units': var_units, 'shape': (num_output_nodes,) + var_shape}
+        if var not in expr_kwargs:
+            expr_kwargs[var] = {}
+
+        if 'units' not in expr_kwargs[var]:
+            expr_kwargs[var]['units'] = var_units
+        elif var == expr_lhs:
+            meta['units'] = expr_kwargs[var]['units']
+
+        if 'shape' in expr_kwargs[var] and var != expr_lhs:
+            om.issue_warning(f'{phase.pathname}: User-provided shape for timeseries expression input `{var}` ignored.\n'
+                             f'Using automatically determined shape {(num_output_nodes,) + var_shape}.',
+                             category=om.UnusedOptionWarning)
+        elif 'shape' in expr_kwargs[var] and var == expr_lhs:
+            if expr_options['shape'] is not _unspecified:
+                om.issue_warning(f'{phase.pathname}: User-provided shape for timeseries expression output `{var}`'
+                                 f'using both the\n`shape` keyword argument and the `{var}` keyword argument '
+                                 f'dictionary.\nThe latter will take precedence.',
+                                 category=om.SetupWarning)
+            var_shape = expr_kwargs[var]['shape']
+            meta['shape'] = var_shape
+
+        expr_kwargs[var]['shape'] = (num_output_nodes,) + var_shape
 
     for input_var in expr_vars:
         abs_name = phase.timeseries_ec_vars[input_var]['abs_name']

From d0f3af2c2710ca33f083bddaff2d8c19c6a54855 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Thu, 3 Nov 2022 08:31:45 -0500
Subject: [PATCH 19/34] added a kwargs dict thaat is redeclared each time in
 introspection to avoid mangling absolute shape and relative shape

---
 dymos/phase/phase.py                       | 2 +-
 dymos/phase/test/test_timeseries.py        | 2 +-
 dymos/transcriptions/transcription_base.py | 2 ++
 dymos/utils/introspection.py               | 7 +++++--
 4 files changed, 9 insertions(+), 4 deletions(-)

diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index c5ab48938..023549b3e 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -1220,7 +1220,7 @@ def add_timeseries_output(self, name, output_name=None, units=_unspecified, shap
         timeseries : str or None
             The name of the timeseries to which the output is being added.
         **kwargs
-            Additional arguments passed to the exec
+            Additional arguments passed to the exec comp.
         """
         if type(name) is list:
             for i, name_i in enumerate(name):
diff --git a/dymos/phase/test/test_timeseries.py b/dymos/phase/test/test_timeseries.py
index d49ae77ac..f317c9240 100644
--- a/dymos/phase/test/test_timeseries.py
+++ b/dymos/phase/test/test_timeseries.py
@@ -600,7 +600,7 @@ def test_ignore_shape(self):
 
         phase.add_timeseries_output('sin_theta=sin(theta)', units='unitless', theta={'units': 'rad', 'shape': (1,)})
 
-        expected = 'User-provided shape for timeseries expression input `theta` ignored.\n' \
+        expected = 'phase0: User-provided shape for timeseries expression input `theta` ignored.\n' \
                    'Using automatically determined shape (20, 1).'
 
         with warnings.catch_warnings(record=True) as ctx:
diff --git a/dymos/transcriptions/transcription_base.py b/dymos/transcriptions/transcription_base.py
index e994afb22..4d4ac8709 100644
--- a/dymos/transcriptions/transcription_base.py
+++ b/dymos/transcriptions/transcription_base.py
@@ -389,6 +389,8 @@ def configure_timeseries_outputs(self, phase):
                 shape = ts_output['shape']
                 src = ts_output['src']
                 is_rate = ts_output['is_rate']
+                print(name)
+                print(shape)
 
                 added_src = timeseries_comp._add_output_configure(name,
                                                                   shape=shape,
diff --git a/dymos/utils/introspection.py b/dymos/utils/introspection.py
index 961e96aba..0a35e89c6 100644
--- a/dymos/utils/introspection.py
+++ b/dymos/utils/introspection.py
@@ -807,6 +807,7 @@ def configure_timeseries_expr_introspection(phase, timeseries_name, expr_options
             phase.timeseries_ec_vars[var_name] = {'added_to_ec': False, 'abs_name': name}
     expr_vars = {}
     expr_kwargs = expr_options['expr_kwargs']
+    kwargs_rel_shape = {}
     meta = {'units': units, 'shape': shape}
 
     for var in phase.timeseries_ec_vars:
@@ -852,14 +853,16 @@ def configure_timeseries_expr_introspection(phase, timeseries_name, expr_options
             var_shape = expr_kwargs[var]['shape']
             meta['shape'] = var_shape
 
-        expr_kwargs[var]['shape'] = (num_output_nodes,) + var_shape
+        expr_kwargs[var]['shape'] = var_shape
+        kwargs_rel_shape[var] = {'units': expr_kwargs[var]['units'],
+                                 'shape': (num_output_nodes,) + expr_kwargs[var]['shape']}
 
     for input_var in expr_vars:
         abs_name = phase.timeseries_ec_vars[input_var]['abs_name']
         if '.' in abs_name:
             expr_reduced = expr_reduced.replace(abs_name, input_var)
 
-    timeseries_ec.add_expr(expr_reduced, **expr_kwargs)
+    timeseries_ec.add_expr(expr_reduced, **kwargs_rel_shape)
 
     for var, options in expr_vars.items():
         if phase.timeseries_ec_vars[var]['added_to_ec']:

From a9920379f587425c068a909f0f5ca18d34c753ae Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Thu, 3 Nov 2022 09:54:08 -0500
Subject: [PATCH 20/34] Added a cannonball problem test case

---
 .../test/test_cannonball_expr_constraint.py   | 132 ++++++++++++++++++
 1 file changed, 132 insertions(+)
 create mode 100644 dymos/examples/cannonball/test/test_cannonball_expr_constraint.py

diff --git a/dymos/examples/cannonball/test/test_cannonball_expr_constraint.py b/dymos/examples/cannonball/test/test_cannonball_expr_constraint.py
new file mode 100644
index 000000000..dc4ab58c8
--- /dev/null
+++ b/dymos/examples/cannonball/test/test_cannonball_expr_constraint.py
@@ -0,0 +1,132 @@
+import unittest
+import matplotlib.pyplot as plt
+
+from openmdao.utils.testing_utils import use_tempdirs
+import openmdao.api as om
+from openmdao.utils.assert_utils import assert_near_equal
+
+import dymos as dm
+from dymos.examples.cannonball.size_comp import CannonballSizeComp
+from dymos.examples.cannonball.cannonball_ode import CannonballODE
+
+
+plt.switch_backend('Agg')
+
+
+@use_tempdirs
+class TestCannonballBoundaryConstraint(unittest.TestCase):
+
+    def test_cannonball_design_boundary_constraint_expression(self):
+
+        p = om.Problem(model=om.Group())
+
+        p.driver = om.pyOptSparseDriver()
+        p.driver.options['optimizer'] = 'SLSQP'
+        p.driver.declare_coloring()
+
+        p.model.add_subsystem('size_comp', CannonballSizeComp(),
+                              promotes_inputs=['radius', 'dens'])
+        p.model.set_input_defaults('dens', val=7.87, units='g/cm**3')
+        p.model.add_design_var('radius', lower=0.01, upper=0.10,
+                               ref0=0.01, ref=0.10, units='m')
+
+        traj = p.model.add_subsystem('traj', dm.Trajectory())
+
+        transcription = dm.Radau(num_segments=5, order=3, compressed=True)
+        ascent = dm.Phase(ode_class=CannonballODE, transcription=transcription)
+
+        ascent = traj.add_phase('ascent', ascent)
+
+        ascent.set_time_options(fix_initial=True, duration_bounds=(1, 100),
+                                duration_ref=100, units='s')
+
+        ascent.set_time_options(fix_initial=True, duration_bounds=(1, 100), duration_ref=100, units='s')
+        ascent.add_state('r', fix_initial=True, fix_final=False, rate_source='r_dot', units='m')
+        ascent.add_state('h', fix_initial=True, fix_final=False, units='m', rate_source='h_dot')
+        ascent.add_state('gam', fix_initial=False, fix_final=True, units='rad', rate_source='gam_dot')
+        ascent.add_state('v', fix_initial=False, fix_final=False, units='m/s', rate_source='v_dot')
+
+        ascent.add_parameter('S', units='m**2', dynamic=False)
+        ascent.add_parameter('m', units='kg', dynamic=False)
+
+        # Limit the muzzle energy
+        ascent.add_boundary_constraint('ke = 0.5 * m * v**2', loc='initial',
+                                       upper=400000, lower=0, ref=100000)
+
+        # Second Phase (descent)
+        transcription = dm.GaussLobatto(num_segments=5, order=3, compressed=True)
+        descent = dm.Phase(ode_class=CannonballODE, transcription=transcription)
+
+        traj.add_phase('descent', descent)
+
+        descent.set_time_options(initial_bounds=(.5, 100), duration_bounds=(.5, 100),
+                                 duration_ref=100, units='s')
+        descent.add_state('r', units='m', rate_source='r_dot')
+        descent.add_state('h', units='m', rate_source='h_dot', fix_initial=False, fix_final=True)
+        descent.add_state('gam', units='rad', rate_source='gam_dot', fix_initial=False, fix_final=False)
+        descent.add_state('v', units='m/s', rate_source='v_dot', fix_initial=False, fix_final=False)
+
+        descent.add_parameter('S', units='m**2', dynamic=False)
+        descent.add_parameter('m', units='kg', dynamic=False)
+
+        descent.add_objective('r', loc='final', scaler=-1.0)
+
+        traj.add_parameter('CD',
+                           targets={'ascent': ['CD'], 'descent': ['CD']},
+                           val=0.5, units=None, opt=False, dynamic=False)
+
+        traj.add_parameter('m', units='kg', val=1.0,
+                           targets={'ascent': 'mass', 'descent': 'mass'}, dynamic=False)
+
+        traj.add_parameter('S', units='m**2', val=0.005, dynamic=False)
+
+        # Link Phases (link time and all state variables)
+        traj.link_phases(phases=['ascent', 'descent'], vars=['*'])
+
+        p.model.connect('size_comp.mass', 'traj.parameters:m')
+        p.model.connect('size_comp.S', 'traj.parameters:S')
+
+        p.model.linear_solver = om.DirectSolver()
+
+        # Finish Problem Setup
+        p.setup()
+
+        # Set constants and initial guesses
+        p.set_val('radius', 0.05, units='m')
+        p.set_val('dens', 7.87, units='g/cm**3')
+
+        p.set_val('traj.parameters:CD', 0.5)
+
+        p.set_val('traj.ascent.t_initial', 0.0)
+        p.set_val('traj.ascent.t_duration', 10.0)
+
+        p.set_val('traj.ascent.states:r', ascent.interpolate(ys=[0, 100],
+                                                             nodes='state_input'))
+        p.set_val('traj.ascent.states:h', ascent.interpolate(ys=[0, 100],
+                                                             nodes='state_input'))
+        p.set_val('traj.ascent.states:v', ascent.interpolate(ys=[200, 150],
+                                                             nodes='state_input'))
+        p.set_val('traj.ascent.states:gam', ascent.interpolate(ys=[25, 0],
+                                                               nodes='state_input'), units='deg')
+
+        p.set_val('traj.descent.t_initial', 10.0)
+        p.set_val('traj.descent.t_duration', 10.0)
+
+        p.set_val('traj.descent.states:r', descent.interpolate(ys=[100, 200],
+                                                               nodes='state_input'))
+        p.set_val('traj.descent.states:h', descent.interpolate(ys=[100, 0],
+                                                               nodes='state_input'))
+        p.set_val('traj.descent.states:v', descent.interpolate(ys=[150, 200],
+                                                               nodes='state_input'))
+        p.set_val('traj.descent.states:gam', descent.interpolate(ys=[0, -45],
+                                                                 nodes='state_input'), units='deg')
+
+        # Run the optimization and final explicit simulation
+        dm.run_problem(p)
+
+        assert_near_equal(p.get_val('traj.descent.states:r')[-1],
+                          3183.25, tolerance=1.0E-2)
+
+
+if __name__ == '__main__':
+    unittest.main()
\ No newline at end of file

From 5517178839ac92e94c5af8928d932e0fc0796ea6 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Thu, 3 Nov 2022 10:25:35 -0500
Subject: [PATCH 21/34] Added checks to see if a constraint is an expression

---
 dymos/phase/options.py |  4 ++++
 dymos/phase/phase.py   | 24 +++++++++++++++++++++---
 2 files changed, 25 insertions(+), 3 deletions(-)

diff --git a/dymos/phase/options.py b/dymos/phase/options.py
index c6b021f5d..66ca173cb 100644
--- a/dymos/phase/options.py
+++ b/dymos/phase/options.py
@@ -666,6 +666,10 @@ def __init__(self, read_only=False):
                      desc='If True, the given indices will be treated as indices into a C-order flattened array based '
                           'on the shaped of the constrained variable at a point in time.')
 
+        self.declare(name='is_expr', types=bool, default=False,
+                     desc='If True, the given constraint is an expression that must be evaluated rather than a'
+                          ' single variable.')
+
 
 class TimeseriesOutputOptionsDictionary(om.OptionsDictionary):
     """
diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index 023549b3e..20fac6af2 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -1028,9 +1028,9 @@ def add_boundary_constraint(self, name, loc, constraint_name=None, units=None,
         loc : str
             The location of the boundary constraint ('initial' or 'final').
         constraint_name : str or None
-            The name of the variable as provided to the boundary constraint comp.  By
-            default this is the last element in `name` when split by dots.  The user may
-            override the constraint name if splitting the path causes name collisions.
+            The name of the boundary constraint. By default, this is the left-hand side of
+            the given expression or "var_constraint" if var is a single variable. The user may
+            override the constraint name if desired.
         units : str or None
             The units in which the boundary constraint is to be applied.  If None, use the
             units associated with the constrained output.  If provided, must be compatible with
@@ -1071,6 +1071,16 @@ def add_boundary_constraint(self, name, loc, constraint_name=None, units=None,
             raise ValueError(f'Invalid boundary constraint location "{loc}". Must be '
                              '"initial" or "final".')
 
+        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
+
         if constraint_name is None:
             constraint_name = name.rpartition('.')[-1]
 
@@ -1082,6 +1092,13 @@ def add_boundary_constraint(self, name, loc, constraint_name=None, units=None,
             raise ValueError(f'Cannot add new {loc} boundary constraint for variable `{name}` and indices {indices}. '
                              f'One already exists.')
 
+        existing_bc = [bc for bc in bc_list if bc['name'] == constraint_name and
+                       bc['indices'] is None and indices is None]
+
+        if existing_bc:
+            raise ValueError(f'Cannot add new {loc} boundary constraint named `{constraint_name}` and indices{indices}.'
+                             f' `{constraint_name}` is already in use as a {loc} boundary constraint')
+
         bc = ConstraintOptionsDictionary()
         bc_list.append(bc)
 
@@ -1099,6 +1116,7 @@ def add_boundary_constraint(self, name, loc, constraint_name=None, units=None,
         bc['linear'] = linear
         bc['units'] = units
         bc['flat_indices'] = flat_indices
+        bc['is_expr'] = is_expr
 
         # Automatically add the requested variable to the timeseries outputs if it's an ODE output.
         var_type = self.classify_var(name)

From 396f2d37c59674dd1eca4dcb747e06b30c73d442 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Thu, 3 Nov 2022 10:26:27 -0500
Subject: [PATCH 22/34] Removed print statements from debugging

---
 dymos/transcriptions/transcription_base.py | 2 --
 1 file changed, 2 deletions(-)

diff --git a/dymos/transcriptions/transcription_base.py b/dymos/transcriptions/transcription_base.py
index 4d4ac8709..e994afb22 100644
--- a/dymos/transcriptions/transcription_base.py
+++ b/dymos/transcriptions/transcription_base.py
@@ -389,8 +389,6 @@ def configure_timeseries_outputs(self, phase):
                 shape = ts_output['shape']
                 src = ts_output['src']
                 is_rate = ts_output['is_rate']
-                print(name)
-                print(shape)
 
                 added_src = timeseries_comp._add_output_configure(name,
                                                                   shape=shape,

From b374b6cf8f7607b19910c0ce315f6f475505c6be Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Tue, 29 Nov 2022 13:04:02 -0600
Subject: [PATCH 23/34] Moved introspection of expressions outside of
 timeseries introspection. Added capability for boundary constraints to find
 timeseries expr variables

---
 .../test/test_cannonball_expr_constraint.py   |  28 +--
 dymos/phase/phase.py                          |  10 +-
 dymos/phase/test/test_timeseries.py           |   4 +-
 dymos/trajectory/trajectory.py                |   2 +-
 .../pseudospectral/pseudospectral_base.py     |   5 +
 dymos/transcriptions/solve_ivp/solve_ivp.py   |   2 +-
 dymos/transcriptions/transcription_base.py    |   5 +-
 dymos/utils/introspection.py                  | 208 +++++++++---------
 8 files changed, 135 insertions(+), 129 deletions(-)

diff --git a/dymos/examples/cannonball/test/test_cannonball_expr_constraint.py b/dymos/examples/cannonball/test/test_cannonball_expr_constraint.py
index dc4ab58c8..cdf0015f4 100644
--- a/dymos/examples/cannonball/test/test_cannonball_expr_constraint.py
+++ b/dymos/examples/cannonball/test/test_cannonball_expr_constraint.py
@@ -76,7 +76,7 @@ def test_cannonball_design_boundary_constraint_expression(self):
                            val=0.5, units=None, opt=False, dynamic=False)
 
         traj.add_parameter('m', units='kg', val=1.0,
-                           targets={'ascent': 'mass', 'descent': 'mass'}, dynamic=False)
+                           targets={'ascent': 'm', 'descent': 'm'}, dynamic=False)
 
         traj.add_parameter('S', units='m**2', val=0.005, dynamic=False)
 
@@ -100,26 +100,18 @@ def test_cannonball_design_boundary_constraint_expression(self):
         p.set_val('traj.ascent.t_initial', 0.0)
         p.set_val('traj.ascent.t_duration', 10.0)
 
-        p.set_val('traj.ascent.states:r', ascent.interpolate(ys=[0, 100],
-                                                             nodes='state_input'))
-        p.set_val('traj.ascent.states:h', ascent.interpolate(ys=[0, 100],
-                                                             nodes='state_input'))
-        p.set_val('traj.ascent.states:v', ascent.interpolate(ys=[200, 150],
-                                                             nodes='state_input'))
-        p.set_val('traj.ascent.states:gam', ascent.interpolate(ys=[25, 0],
-                                                               nodes='state_input'), units='deg')
+        p.set_val('traj.ascent.states:r', ascent.interp(ys=[0, 100], nodes='state_input'))
+        p.set_val('traj.ascent.states:h', ascent.interp(ys=[0, 100], nodes='state_input'))
+        p.set_val('traj.ascent.states:v', ascent.interp(ys=[200, 150], nodes='state_input'))
+        p.set_val('traj.ascent.states:gam', ascent.interp(ys=[25, 0], nodes='state_input'), units='deg')
 
         p.set_val('traj.descent.t_initial', 10.0)
         p.set_val('traj.descent.t_duration', 10.0)
 
-        p.set_val('traj.descent.states:r', descent.interpolate(ys=[100, 200],
-                                                               nodes='state_input'))
-        p.set_val('traj.descent.states:h', descent.interpolate(ys=[100, 0],
-                                                               nodes='state_input'))
-        p.set_val('traj.descent.states:v', descent.interpolate(ys=[150, 200],
-                                                               nodes='state_input'))
-        p.set_val('traj.descent.states:gam', descent.interpolate(ys=[0, -45],
-                                                                 nodes='state_input'), units='deg')
+        p.set_val('traj.descent.states:r', descent.interp(ys=[100, 200], nodes='state_input'))
+        p.set_val('traj.descent.states:h', descent.interp(ys=[100, 0], nodes='state_input'))
+        p.set_val('traj.descent.states:v', descent.interp(ys=[150, 200], nodes='state_input'))
+        p.set_val('traj.descent.states:gam', descent.interp(ys=[0, -45], nodes='state_input'), units='deg')
 
         # Run the optimization and final explicit simulation
         dm.run_problem(p)
@@ -129,4 +121,4 @@ def test_cannonball_design_boundary_constraint_expression(self):
 
 
 if __name__ == '__main__':
-    unittest.main()
\ No newline at end of file
+    unittest.main()
diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index 20fac6af2..bf62c58e7 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -22,7 +22,7 @@
 from ..utils.indexing import get_constraint_flat_idxs
 from ..utils.introspection import configure_time_introspection, _configure_constraint_introspection, \
     configure_controls_introspection, configure_parameters_introspection, \
-    configure_timeseries_output_introspection, classify_var, get_promoted_vars
+    configure_timeseries_output_introspection, classify_var, configure_timeseries_expr_introspection
 from ..utils.misc import _unspecified
 from ..utils.lgl import lgl
 
@@ -1081,7 +1081,9 @@ def add_boundary_constraint(self, name, loc, constraint_name=None, units=None,
         else:
             is_expr = False
 
-        if constraint_name is None:
+        if is_expr:
+            constraint_name = name.split('=')[0].strip()
+        elif constraint_name is None:
             constraint_name = name.rpartition('.')[-1]
 
         bc_list = self._initial_boundary_constraints if loc == 'initial' else self._final_boundary_constraints
@@ -1721,7 +1723,9 @@ def configure(self):
 
         _configure_constraint_introspection(self)
 
-        transcription._configure_boundary_constraints(self)
+        configure_timeseries_expr_introspection(self)
+
+        transcription.configure_boundary_constraints(self)
 
         transcription.configure_path_constraints(self)
 
diff --git a/dymos/phase/test/test_timeseries.py b/dymos/phase/test/test_timeseries.py
index f317c9240..db1e03ed8 100644
--- a/dymos/phase/test/test_timeseries.py
+++ b/dymos/phase/test/test_timeseries.py
@@ -483,7 +483,7 @@ def test_invalid_expr_var(self):
 
         expected = "phase0: Timeseries expression `k=units` uses the following disallowed variable " \
                    "name in its definition: 'units'."
-        self.assertEqual(str(e.exception.__cause__), expected)
+        self.assertEqual(str(e.exception), expected)
 
     def test_input_units(self):
         p = om.Problem(model=om.Group())
@@ -549,7 +549,7 @@ def test_output_units(self):
         phase.add_timeseries_output('sin_theta2=sin(theta)', sin_theta2={'units': 'unitless'})
         phase.add_timeseries_output('sin_theta3=sin(theta)', shape=(1,), sin_theta3={'shape': (1,)})
 
-        expected = 'phase0: User-provided shape for timeseries expression output `sin_theta3`using ' \
+        expected = 'phase0: User-provided shape for timeseries expression output `sin_theta3` using ' \
                    'both the\n`shape` keyword argument and the `sin_theta3` keyword argument dictionary.\n' \
                    'The latter will take precedence.'
 
diff --git a/dymos/trajectory/trajectory.py b/dymos/trajectory/trajectory.py
index 73ce8335e..d7f4f05d5 100644
--- a/dymos/trajectory/trajectory.py
+++ b/dymos/trajectory/trajectory.py
@@ -1082,7 +1082,7 @@ def _print_constraints(phase, outstream):
             for loc, d in ds.items():
                 str_loc = f'[{loc}]'
                 for options in d:
-                    expr = options['name']
+                    expr = options['constraint_name']
                     _, shape, units, linear = tx._get_objective_src(expr, loc, phase, ode_outputs=ode_outputs)
 
                     equals = options['equals']
diff --git a/dymos/transcriptions/pseudospectral/pseudospectral_base.py b/dymos/transcriptions/pseudospectral/pseudospectral_base.py
index c2ea41c09..b52090813 100644
--- a/dymos/transcriptions/pseudospectral/pseudospectral_base.py
+++ b/dymos/transcriptions/pseudospectral/pseudospectral_base.py
@@ -602,6 +602,11 @@ def _get_objective_src(self, var, loc, phase, ode_outputs=None):
             units = control_rate_units
             linear = False
             constraint_path = f'polynomial_control_rates:{var}'
+        elif var_type == 'timeseries_exec_comp_output':
+            shape = (1,)
+            units = None
+            constraint_path = f'timeseries.timeseries_exec_comp.{var}'
+            linear = False
         else:
             # Failed to find variable, assume it is in the ODE. This requires introspection.
             constraint_path = f'{self._rhs_source}.{var}'
diff --git a/dymos/transcriptions/solve_ivp/solve_ivp.py b/dymos/transcriptions/solve_ivp/solve_ivp.py
index 95aacced4..bfab091fa 100644
--- a/dymos/transcriptions/solve_ivp/solve_ivp.py
+++ b/dymos/transcriptions/solve_ivp/solve_ivp.py
@@ -460,7 +460,7 @@ def configure_path_constraints(self, phase):
         """
         pass
 
-    def _configure_boundary_constraints(self, phase):
+    def configure_boundary_constraints(self, phase):
         """
         Not used in SolveIVP.
 
diff --git a/dymos/transcriptions/transcription_base.py b/dymos/transcriptions/transcription_base.py
index e994afb22..93d2d57f7 100644
--- a/dymos/transcriptions/transcription_base.py
+++ b/dymos/transcriptions/transcription_base.py
@@ -401,7 +401,7 @@ def configure_timeseries_outputs(self, phase):
                     phase.connect(src_name=src, tgt_name=f'{timeseries_name}.input_values:{name}',
                                   src_indices=ts_output['src_idxs'])
 
-    def _configure_boundary_constraints(self, phase):
+    def configure_boundary_constraints(self, phase):
         """
         Configures the boundary constraints.
 
@@ -448,7 +448,7 @@ def _get_constraint_kwargs(self, constraint_type, options, phase):
         con_name = constraint_kwargs.pop('constraint_name')
 
         # Determine the path to the variable which we will be constraining
-        var = options['name']
+        var = con_name if options['is_expr'] else options['name']
         var_type = phase.classify_var(var)
 
         # These are the flat indices at a single point in time used
@@ -503,6 +503,7 @@ def _get_constraint_kwargs(self, constraint_type, options, phase):
         con_path = constraint_kwargs.pop('constraint_path')
         constraint_kwargs.pop('shape')
         constraint_kwargs['flat_indices'] = True
+        constraint_kwargs.pop('is_expr')
 
         return con_path, constraint_kwargs
 
diff --git a/dymos/utils/introspection.py b/dymos/utils/introspection.py
index 0a35e89c6..3331546c9 100644
--- a/dymos/utils/introspection.py
+++ b/dymos/utils/introspection.py
@@ -177,7 +177,7 @@ def _configure_constraint_introspection(phase):
             time_units = phase.time_options['units']
 
             # Determine the path to the variable which we will be constraining
-            var = con['name']
+            var = con['constraint_name'] if con['is_expr'] else con['name']
             var_type = phase.classify_var(var)
 
             if var != con['constraint_name'] is not None and var_type != 'ode':
@@ -269,6 +269,11 @@ def _configure_constraint_introspection(phase):
                     if con['units'] is None else con['units']
                 con['constraint_path'] = f'timeseries.polynomial_control_rates:{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)
@@ -736,7 +741,7 @@ def configure_timeseries_output_introspection(phase):
 
         for output_name, output_options in ts_opts['outputs'].items():
             if output_options['is_expr']:
-                output_meta = configure_timeseries_expr_introspection(phase, ts_name, output_options)
+                output_meta = phase.timeseries_ec_vars[ts_name][output_name]['meta_data']
             else:
                 try:
                     output_meta = transcription._get_timeseries_var_source(output_options['name'],
@@ -763,7 +768,7 @@ def configure_timeseries_output_introspection(phase):
             ts_opts['outputs'].pop(s)
 
 
-def configure_timeseries_expr_introspection(phase, timeseries_name, expr_options):
+def configure_timeseries_expr_introspection(phase):
     """
     Introspect the timeseries expressions, add expressions to the exec comp, and make connections to the variables.
 
@@ -771,111 +776,110 @@ def configure_timeseries_expr_introspection(phase, timeseries_name, expr_options
     ----------
     phase : Phase
         The phase object whose timeseries outputs are to be introspected.
-    timeseries_name : str
-        The name of the timeseries that the expr is to be added to.
-    expr_options : dict
-        Options for the expression to be added to the timeseries exec comp.
-
-    Returns
-    -------
-    dict
-        A dictionary containing the metadata for the source. This consists of shape, units, and tags.
     """
-    timeseries_ec = phase._get_subsystem(f'{timeseries_name}.timeseries_exec_comp')
     transcription = phase.options['transcription']
     num_output_nodes = transcription._get_num_timeseries_nodes()
 
-    expr = expr_options['name']
-    expr_lhs = expr.split('=')[0].strip()
-    if '.' in expr_lhs or ':' in expr_lhs:
-        raise ValueError(f'{expr_lhs} is an invalid name for the timeseries LHS. Must use a valid variable name')
-    expr_reduced = expr
-
-    units = expr_options['units'] if expr_options['units'] is not _unspecified else None
-    shape = expr_options['shape'] if expr_options['shape'] is not _unspecified else (1,)
-
-    abs_names = [x.strip() for x in re.findall(re.compile(r'([_a-zA-Z]\w*[ ]*\(?:?[.]?)'), expr)
-                 if not x.endswith('(') and not x.endswith(':')]
-    for name in abs_names:
-        if name.endswith('.'):
-            idx = abs_names.index(name)
-            abs_names[idx + 1] = name + abs_names[idx + 1]
-            abs_names.remove(name)
-    for name in abs_names:
-        var_name = name.split('.')[-1]
-        if var_name not in phase.timeseries_ec_vars.keys():
-            phase.timeseries_ec_vars[var_name] = {'added_to_ec': False, 'abs_name': name}
-    expr_vars = {}
-    expr_kwargs = expr_options['expr_kwargs']
-    kwargs_rel_shape = {}
-    meta = {'units': units, 'shape': shape}
-
-    for var in phase.timeseries_ec_vars:
-        if var in ('output_name', 'units', 'shape', 'timeseries'):
-            raise RuntimeError(f'{phase.pathname}: Timeseries expression `{expr}` uses the following disallowed '
-                               f'variable name in its definition: \'{var}\'.')
-
-        if phase.timeseries_ec_vars[var]['added_to_ec']:
-            continue
-        elif var == expr_lhs:
-            var_units = units
-            var_shape = shape
-            meta['src'] = f'{timeseries_name}.timeseries_exec_comp.{var}'
-            phase.timeseries_ec_vars[var]['added_to_ec'] = True
-        else:
-            try:
-                expr_vars[var] = transcription._get_timeseries_var_source(var, output_name=var, phase=phase)
-            except ValueError as e:
-                expr_vars[var] = transcription._get_timeseries_var_source(phase.timeseries_ec_vars[var]['abs_name'],
-                                                                          output_name=var, phase=phase)
-
-            var_units = expr_vars[var]['units']
-            var_shape = expr_vars[var]['shape']
-
-        if var not in expr_kwargs:
-            expr_kwargs[var] = {}
-
-        if 'units' not in expr_kwargs[var]:
-            expr_kwargs[var]['units'] = var_units
-        elif var == expr_lhs:
-            meta['units'] = expr_kwargs[var]['units']
-
-        if 'shape' in expr_kwargs[var] and var != expr_lhs:
-            om.issue_warning(f'{phase.pathname}: User-provided shape for timeseries expression input `{var}` ignored.\n'
-                             f'Using automatically determined shape {(num_output_nodes,) + var_shape}.',
-                             category=om.UnusedOptionWarning)
-        elif 'shape' in expr_kwargs[var] and var == expr_lhs:
-            if expr_options['shape'] is not _unspecified:
-                om.issue_warning(f'{phase.pathname}: User-provided shape for timeseries expression output `{var}`'
-                                 f'using both the\n`shape` keyword argument and the `{var}` keyword argument '
-                                 f'dictionary.\nThe latter will take precedence.',
-                                 category=om.SetupWarning)
-            var_shape = expr_kwargs[var]['shape']
-            meta['shape'] = var_shape
-
-        expr_kwargs[var]['shape'] = var_shape
-        kwargs_rel_shape[var] = {'units': expr_kwargs[var]['units'],
-                                 'shape': (num_output_nodes,) + expr_kwargs[var]['shape']}
-
-    for input_var in expr_vars:
-        abs_name = phase.timeseries_ec_vars[input_var]['abs_name']
-        if '.' in abs_name:
-            expr_reduced = expr_reduced.replace(abs_name, input_var)
-
-    timeseries_ec.add_expr(expr_reduced, **kwargs_rel_shape)
-
-    for var, options in expr_vars.items():
-        if phase.timeseries_ec_vars[var]['added_to_ec']:
-            continue
-        else:
-            phase.connect(src_name=options['src'], tgt_name=f'{timeseries_name}.timeseries_exec_comp.{var}',
-                          src_indices=options['src_idxs'])
-
-        phase.timeseries_ec_vars[var]['added_to_ec'] = True
+    for ts_name, ts_opts in phase._timeseries.items():
+        timeseries_ec = phase._get_subsystem(f'{ts_name}.timeseries_exec_comp')
+        phase.timeseries_ec_vars[ts_name] = {}
 
-    meta['src_idxs'] = None
+        for output_name, output_options in ts_opts['outputs'].items():
+            if output_options['is_expr']:
 
-    return meta
+                expr = output_options['name']
+                expr_lhs = expr.split('=')[0].strip()
+                if '.' in expr_lhs or ':' in expr_lhs:
+                    raise ValueError(
+                        f'{expr_lhs} is an invalid name for the timeseries LHS. Must use a valid variable name')
+                expr_reduced = expr
+
+                units = output_options['units'] if output_options['units'] is not _unspecified else None
+                shape = output_options['shape'] if output_options['shape'] not in {_unspecified, None} else (1,)
+
+                abs_names = [x.strip() for x in re.findall(re.compile(r'([_a-zA-Z]\w*[ ]*\(?:?[.]?)'), expr)
+                             if not x.endswith('(') and not x.endswith(':')]
+                for name in abs_names:
+                    if name.endswith('.'):
+                        idx = abs_names.index(name)
+                        abs_names[idx + 1] = name + abs_names[idx + 1]
+                        abs_names.remove(name)
+                for name in abs_names:
+                    var_name = name.split('.')[-1]
+                    if var_name not in phase.timeseries_ec_vars[ts_name].keys():
+                        phase.timeseries_ec_vars[ts_name][var_name] = {'added_to_ec': False, 'abs_name': name}
+                expr_vars = {}
+                expr_kwargs = output_options['expr_kwargs']
+                kwargs_rel_shape = {}
+                meta = {'units': units, 'shape': shape}
+
+                for var in phase.timeseries_ec_vars[ts_name]:
+                    if var in ('output_name', 'units', 'shape', 'timeseries'):
+                        raise RuntimeError(f'{phase.pathname}: Timeseries expression `{expr}` uses the following'
+                                           f' disallowed variable name in its definition: \'{var}\'.')
+
+                    if phase.timeseries_ec_vars[ts_name][var]['added_to_ec']:
+                        continue
+                    elif var == expr_lhs:
+                        var_units = units
+                        var_shape = shape
+                        meta['src'] = f'{ts_name}.timeseries_exec_comp.{var}'
+                        phase.timeseries_ec_vars[ts_name][var]['added_to_ec'] = True
+                    else:
+                        try:
+                            expr_vars[var] = transcription._get_timeseries_var_source(var, output_name=var, phase=phase)
+                        except ValueError as e:
+                            expr_vars[var] = transcription._get_timeseries_var_source(
+                                phase.timeseries_ec_vars[ts_name][var]['abs_name'],
+                                output_name=var, phase=phase)
+
+                        var_units = expr_vars[var]['units']
+                        var_shape = expr_vars[var]['shape']
+
+                    if var not in expr_kwargs:
+                        expr_kwargs[var] = {}
+
+                    if 'units' not in expr_kwargs[var]:
+                        expr_kwargs[var]['units'] = var_units
+                    elif var == expr_lhs:
+                        meta['units'] = expr_kwargs[var]['units']
+
+                    if 'shape' in expr_kwargs[var] and var != expr_lhs:
+                        om.issue_warning(f'{phase.pathname}: User-provided shape for timeseries expression input'
+                                         f' `{var}` ignored.\nUsing automatically determined shape'
+                                         f' {(num_output_nodes,) + var_shape}.',
+                                         category=om.UnusedOptionWarning)
+                    elif 'shape' in expr_kwargs[var] and var == expr_lhs:
+                        if output_options['shape'] is not _unspecified:
+                            om.issue_warning(f'{phase.pathname}: User-provided shape for timeseries expression output'
+                                             f' `{var}` using both the\n`shape` keyword argument and the `{var}`'
+                                             f' keyword argument dictionary.\nThe latter will take precedence.',
+                                             category=om.SetupWarning)
+                        var_shape = expr_kwargs[var]['shape']
+                        meta['shape'] = var_shape
+
+                    expr_kwargs[var]['shape'] = var_shape
+                    kwargs_rel_shape[var] = {'units': expr_kwargs[var]['units'],
+                                             'shape': (num_output_nodes,) + expr_kwargs[var]['shape']}
+
+                for input_var in expr_vars:
+                    abs_name = phase.timeseries_ec_vars[ts_name][input_var]['abs_name']
+                    if '.' in abs_name:
+                        expr_reduced = expr_reduced.replace(abs_name, input_var)
+
+                timeseries_ec.add_expr(expr_reduced, **kwargs_rel_shape)
+
+                for var, options in expr_vars.items():
+                    if phase.timeseries_ec_vars[ts_name][var]['added_to_ec']:
+                        continue
+                    else:
+                        phase.connect(src_name=options['src'], tgt_name=f'{ts_name}.timeseries_exec_comp.{var}',
+                                      src_indices=options['src_idxs'])
+
+                    phase.timeseries_ec_vars[ts_name][var]['added_to_ec'] = True
+
+                meta['src_idxs'] = None
+                phase.timeseries_ec_vars[ts_name][expr_lhs]['meta_data'] = meta
 
 
 def filter_outputs(patterns, sys):

From 5e64e744b14aabd463ffe42d671493dd1b09acd6 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Tue, 29 Nov 2022 16:01:38 -0600
Subject: [PATCH 24/34] Modified add_path_constraint to accept expressions and
 added a test problem

---
 .../brachistochrone/test/test_brachistochrone_path_constraint.py  | 0
 1 file changed, 0 insertions(+), 0 deletions(-)
 create mode 100644 dymos/examples/brachistochrone/test/test_brachistochrone_path_constraint.py

diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_path_constraint.py b/dymos/examples/brachistochrone/test/test_brachistochrone_path_constraint.py
new file mode 100644
index 000000000..e69de29bb

From 775b52080f6e2366840e584d44141b058fa98903 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Tue, 29 Nov 2022 16:01:56 -0600
Subject: [PATCH 25/34] Modified add_path_constraint to accept expressions and
 added a test problem

---
 .../test_brachistochrone_path_constraint.py   | 75 +++++++++++++++++++
 dymos/phase/phase.py                          | 15 +++-
 2 files changed, 89 insertions(+), 1 deletion(-)

diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_path_constraint.py b/dymos/examples/brachistochrone/test/test_brachistochrone_path_constraint.py
index e69de29bb..bcf52c708 100644
--- a/dymos/examples/brachistochrone/test/test_brachistochrone_path_constraint.py
+++ b/dymos/examples/brachistochrone/test/test_brachistochrone_path_constraint.py
@@ -0,0 +1,75 @@
+import os
+import unittest
+
+import matplotlib
+import openmdao.api as om
+import matplotlib.pyplot as plt
+import dymos as dm
+
+from openmdao.utils.testing_utils import use_tempdirs
+from openmdao.utils.assert_utils import assert_near_equal
+from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE
+
+matplotlib.use('Agg')
+plt.style.use('ggplot')
+
+
+@use_tempdirs
+class TestBrachistochroneStaticGravity(unittest.TestCase):
+
+    def _make_problem(self, tx):
+        p = om.Problem(model=om.Group())
+
+        p.driver = om.ScipyOptimizeDriver()
+        p.driver.declare_coloring()
+
+        phase = dm.Phase(ode_class=BrachistochroneODE,
+                         transcription=tx)
+
+        p.model.add_subsystem('phase0', phase)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+
+        phase.set_state_options('x', fix_initial=True)
+        phase.set_state_options('y', fix_initial=True)
+        phase.set_state_options('v', fix_initial=True)
+
+        phase.add_control('theta', continuity=True, rate_continuity=True, opt=True,
+                          units='deg', lower=0.01, upper=179.9, ref=1, ref0=0)
+
+        phase.add_boundary_constraint('x', loc='final', equals=10.0)
+        phase.add_path_constraint('pc = y-x/2-1', lower=0.0)
+
+        phase.add_parameter('g', opt=False, units='m/s**2', val=9.80665, include_timeseries=True)
+
+        phase.add_objective('time_phase', loc='final', scaler=10)
+
+        p.model.options['assembled_jac_type'] = 'csc'
+        p.model.linear_solver = om.DirectSolver()
+        p.setup(check=True)
+
+        p['phase0.t_initial'] = 0.0
+        p['phase0.t_duration'] = 2.0
+
+        p['phase0.states:x'] = phase.interp('x', [0, 10])
+        p['phase0.states:y'] = phase.interp('y', [10, 5])
+        p['phase0.states:v'] = phase.interp('v', [0, 9.9])
+        p[f'phase0.controls:theta'] = phase.interp('theta', [5, 100])
+        p['phase0.parameters:g'] = 9.80665
+        return p
+
+    def tearDown(self):
+        for filename in ['total_coloring.pkl', 'SLSQP.out', 'SNOPT_print.out', 'SNOPT_summary.out']:
+            if os.path.exists(filename):
+                os.remove(filename)
+
+    def test_brachistochrone_path_constraint(self):
+        prob = self._make_problem(tx=dm.Radau(num_segments=5, order=3, compressed=True))
+        prob.run_driver()
+        yf = prob.get_val('phase0.timeseries.states:y')[-1]
+
+        assert_near_equal(yf, 6)
+
+
+if __name__ == '__main__':  # pragma: no cover
+    unittest.main()
diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index bf62c58e7..d995da879 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -1176,7 +1176,19 @@ def add_path_constraint(self, name, constraint_name=None, units=None, shape=None
             If True, treat indices as flattened C-ordered indices of elements to constrain at each given point in time.
             Otherwise, indices should be a tuple or list giving the elements to constrain at each point in time.
         """
-        if constraint_name is None:
+        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
+
+        if is_expr:
+            constraint_name = name.split('=')[0].strip()
+        elif constraint_name is None:
             constraint_name = name.rpartition('.')[-1]
 
         existing_pc = [pc for pc in self._path_constraints
@@ -1203,6 +1215,7 @@ def add_path_constraint(self, name, constraint_name=None, units=None, shape=None
         pc['linear'] = linear
         pc['units'] = units
         pc['flat_indices'] = flat_indices
+        pc['is_expr'] = is_expr
 
         # Automatically add the requested variable to the timeseries outputs if it's an ODE output.
         var_type = self.classify_var(name)

From 3124f3864cef62e01c95f217dd895b601d35ab54 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Wed, 30 Nov 2022 09:26:52 -0600
Subject: [PATCH 26/34] Fixed a bug where incorrect name was used for non
 expression constraints

---
 dymos/trajectory/trajectory.py | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/dymos/trajectory/trajectory.py b/dymos/trajectory/trajectory.py
index 0789e0acc..fb8eadaf8 100644
--- a/dymos/trajectory/trajectory.py
+++ b/dymos/trajectory/trajectory.py
@@ -1112,8 +1112,11 @@ def _print_constraints(phase, outstream):
             for loc, d in ds.items():
                 str_loc = f'[{loc}]'
                 for options in d:
-                    expr = options['constraint_name']
-                    _, shape, units, linear = tx._get_objective_src(expr, loc, phase, ode_outputs=ode_outputs)
+                    if options['is_expr']:
+                        name = options['constraint_name']
+                    else:
+                        name = options['name']
+                    _, shape, units, linear = tx._get_objective_src(name, loc, phase, ode_outputs=ode_outputs)
 
                     equals = options['equals']
                     lower = options['lower']
@@ -1146,11 +1149,11 @@ def _print_constraints(phase, outstream):
                         str_upper = ''
 
                     if equals is not None:
-                        printer(f'{2 * indent}{str_loc:<10s}{str_equals} == {expr} [{str_units}]',
+                        printer(f'{2 * indent}{str_loc:<10s}{str_equals} == {name} [{str_units}]',
                                 file=outstream)
                     else:
                         printer(
-                            f'{2 * indent}{str_loc:<10s}{str_lower} {expr} {str_upper} [{str_units}]',
+                            f'{2 * indent}{str_loc:<10s}{str_lower} {name} {str_upper} [{str_units}]',
                             file=outstream)
 
         for phase_name, phs in self._phases.items():

From 236c974dfcab1bf3b7a5ec933f9e444667f7cf61 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Wed, 30 Nov 2022 11:02:51 -0600
Subject: [PATCH 27/34] Added a check for the case in solve_ivp where
 output_nodes_per_seg is not defined

---
 dymos/transcriptions/solve_ivp/solve_ivp.py | 6 +++++-
 1 file changed, 5 insertions(+), 1 deletion(-)

diff --git a/dymos/transcriptions/solve_ivp/solve_ivp.py b/dymos/transcriptions/solve_ivp/solve_ivp.py
index bfab091fa..3b9adf165 100644
--- a/dymos/transcriptions/solve_ivp/solve_ivp.py
+++ b/dymos/transcriptions/solve_ivp/solve_ivp.py
@@ -686,4 +686,8 @@ def _get_num_timeseries_nodes(self):
         int
             The number of nodes in the default timeseries for this transcription.
         """
-        return self.grid_data.num_segments * self.options['output_nodes_per_seg']
+        if self.options['output_nodes_per_seg'] is None:
+            output_nodes = self.grid_data.num_segments
+        else:
+            output_nodes = self.grid_data.num_segments * self.options['output_nodes_per_seg']
+        return output_nodes

From 4c37c9cdb25ab8db3ba5e5dc39d9bb0fab33b715 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Wed, 30 Nov 2022 11:21:11 -0600
Subject: [PATCH 28/34] changed cannonball example to use scipy driver instead
 of pyoptsparse

---
 .../examples/cannonball/test/test_cannonball_expr_constraint.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/dymos/examples/cannonball/test/test_cannonball_expr_constraint.py b/dymos/examples/cannonball/test/test_cannonball_expr_constraint.py
index cdf0015f4..716e9fb26 100644
--- a/dymos/examples/cannonball/test/test_cannonball_expr_constraint.py
+++ b/dymos/examples/cannonball/test/test_cannonball_expr_constraint.py
@@ -20,7 +20,7 @@ def test_cannonball_design_boundary_constraint_expression(self):
 
         p = om.Problem(model=om.Group())
 
-        p.driver = om.pyOptSparseDriver()
+        p.driver = om.ScipyOptimizeDriver()
         p.driver.options['optimizer'] = 'SLSQP'
         p.driver.declare_coloring()
 

From e2a7bd8c7c15f86796b47f08c05bc7bd44f12c2a Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Thu, 1 Dec 2022 07:58:26 -0600
Subject: [PATCH 29/34] Renamed tests to make it more clear what they were
 testing and added matrix multiplication as an oerator to look for

---
 .../test/test_brachistochrone_path_constraint.py              | 4 ++--
 dymos/phase/phase.py                                          | 4 ++--
 2 files changed, 4 insertions(+), 4 deletions(-)

diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_path_constraint.py b/dymos/examples/brachistochrone/test/test_brachistochrone_path_constraint.py
index bcf52c708..2116aba7d 100644
--- a/dymos/examples/brachistochrone/test/test_brachistochrone_path_constraint.py
+++ b/dymos/examples/brachistochrone/test/test_brachistochrone_path_constraint.py
@@ -15,7 +15,7 @@
 
 
 @use_tempdirs
-class TestBrachistochroneStaticGravity(unittest.TestCase):
+class TestBrachistochroneExprPathConstraint(unittest.TestCase):
 
     def _make_problem(self, tx):
         p = om.Problem(model=om.Group())
@@ -63,7 +63,7 @@ def tearDown(self):
             if os.path.exists(filename):
                 os.remove(filename)
 
-    def test_brachistochrone_path_constraint(self):
+    def test_brachistochrone_expr_path_constraint(self):
         prob = self._make_problem(tx=dm.Radau(num_segments=5, order=3, compressed=True))
         prob.run_driver()
         yf = prob.get_val('phase0.timeseries.states:y')[-1]
diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index 215d3ec29..7f588fff0 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -1126,7 +1126,7 @@ def add_boundary_constraint(self, name, loc, constraint_name=None, units=None,
             raise ValueError(f'Invalid boundary constraint location "{loc}". Must be '
                              '"initial" or "final".')
 
-        expr_operators = ['(', '+', '-', '/', '*', '&', '%']
+        expr_operators = ['(', '+', '-', '/', '*', '&', '%', '@']
         if '=' in name:
             is_expr = True
         elif '=' not in name and any(opr in name for opr in expr_operators):
@@ -1231,7 +1231,7 @@ def add_path_constraint(self, name, constraint_name=None, units=None, shape=None
             If True, treat indices as flattened C-ordered indices of elements to constrain at each given point in time.
             Otherwise, indices should be a tuple or list giving the elements to constrain at each point in time.
         """
-        expr_operators = ['(', '+', '-', '/', '*', '&', '%']
+        expr_operators = ['(', '+', '-', '/', '*', '&', '%', '@']
         if '=' in name:
             is_expr = True
         elif '=' not in name and any(opr in name for opr in expr_operators):

From b1836e03572b77603b350eabe7870baf82f7c205 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Mon, 5 Dec 2022 10:16:22 -0600
Subject: [PATCH 30/34] Moved regex compilation outside loop

---
 dymos/utils/introspection.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/dymos/utils/introspection.py b/dymos/utils/introspection.py
index eacbed14b..c6a7c79e6 100644
--- a/dymos/utils/introspection.py
+++ b/dymos/utils/introspection.py
@@ -790,6 +790,7 @@ def configure_timeseries_expr_introspection(phase):
     """
     transcription = phase.options['transcription']
     num_output_nodes = transcription._get_num_timeseries_nodes()
+    var_names_regex = re.compile(r'([_a-zA-Z]\w*[ ]*\(?:?[.]?)')
 
     for ts_name, ts_opts in phase._timeseries.items():
         timeseries_ec = phase._get_subsystem(f'{ts_name}.timeseries_exec_comp')
@@ -808,7 +809,7 @@ def configure_timeseries_expr_introspection(phase):
                 units = output_options['units'] if output_options['units'] is not _unspecified else None
                 shape = output_options['shape'] if output_options['shape'] not in {_unspecified, None} else (1,)
 
-                abs_names = [x.strip() for x in re.findall(re.compile(r'([_a-zA-Z]\w*[ ]*\(?:?[.]?)'), expr)
+                abs_names = [x.strip() for x in re.findall(var_names_regex, expr)
                              if not x.endswith('(') and not x.endswith(':')]
                 for name in abs_names:
                     if name.endswith('.'):

From 5f36c14abfaa0c9d0e9d48089aaf7dcef63f543b Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Tue, 6 Dec 2022 13:57:05 -0600
Subject: [PATCH 31/34] Added tests to check error messages

---
 .../test_brachistochrone_path_constraint.py   |   5 -
 dymos/phase/phase.py                          |  50 ++--
 .../test/test_add_boundary_constraint.py      | 254 ++++++++++++++++++
 .../test_add_boundary_constraint_array.py     | 118 --------
 dymos/phase/test/test_add_path_constraint.py  | 190 +++++++++++++
 5 files changed, 475 insertions(+), 142 deletions(-)
 create mode 100644 dymos/phase/test/test_add_boundary_constraint.py
 delete mode 100644 dymos/phase/test/test_add_boundary_constraint_array.py
 create mode 100644 dymos/phase/test/test_add_path_constraint.py

diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_path_constraint.py b/dymos/examples/brachistochrone/test/test_brachistochrone_path_constraint.py
index 2116aba7d..2d51bc212 100644
--- a/dymos/examples/brachistochrone/test/test_brachistochrone_path_constraint.py
+++ b/dymos/examples/brachistochrone/test/test_brachistochrone_path_constraint.py
@@ -58,11 +58,6 @@ def _make_problem(self, tx):
         p['phase0.parameters:g'] = 9.80665
         return p
 
-    def tearDown(self):
-        for filename in ['total_coloring.pkl', 'SLSQP.out', 'SNOPT_print.out', 'SNOPT_summary.out']:
-            if os.path.exists(filename):
-                os.remove(filename)
-
     def test_brachistochrone_expr_path_constraint(self):
         prob = self._make_problem(tx=dm.Radau(num_segments=5, order=3, compressed=True))
         prob.run_driver()
diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index 7f588fff0..a81e2f8e8 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -1078,14 +1078,15 @@ def add_boundary_constraint(self, name, loc, constraint_name=None, units=None,
         Parameters
         ----------
         name : str
-            Name of the variable to constrain.  If name is not a state, control, or 'time',
+            Name of the variable to constrain. May also provide an expression to be evaluated and constrained.
+            If a single variable and the name is not a state, control, or 'time',
             then this is assumed to be the path of the variable to be constrained in the ODE.
+            If an expression, it must be provided in the form of an equation with a left- and right-hand side
         loc : str
             The location of the boundary constraint ('initial' or 'final').
         constraint_name : str or None
-            The name of the boundary constraint. By default, this is the left-hand side of
-            the given expression or "var_constraint" if var is a single variable. The user may
-            override the constraint name if desired.
+            The name of the boundary constraint. By default, this is 'var_constraint' if name is a single variable,
+             or the left-hand side of the equation if name is an expression.
         units : str or None
             The units in which the boundary constraint is to be applied.  If None, use the
             units associated with the constrained output.  If provided, must be compatible with
@@ -1130,9 +1131,9 @@ def add_boundary_constraint(self, name, loc, constraint_name=None, units=None,
         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)"')
+            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
 
@@ -1149,12 +1150,13 @@ def add_boundary_constraint(self, name, loc, constraint_name=None, units=None,
             raise ValueError(f'Cannot add new {loc} boundary constraint for variable `{name}` and indices {indices}. '
                              f'One already exists.')
 
-        existing_bc = [bc for bc in bc_list if bc['name'] == constraint_name and
-                       bc['indices'] is None and indices is None]
+        existing_bc_name = [bc for bc in bc_list if bc['name'] == constraint_name and
+                            bc['indices'] is None and indices is None]
 
-        if existing_bc:
-            raise ValueError(f'Cannot add new {loc} boundary constraint named `{constraint_name}` and indices{indices}.'
-                             f' `{constraint_name}` is already in use as a {loc} boundary constraint')
+        if existing_bc_name:
+            raise ValueError(f'Cannot add new {loc} boundary constraint named `{constraint_name}`'
+                             f' and indices {indices}. The name `{constraint_name}` is already in use'
+                             f' as a {loc} boundary constraint')
 
         bc = ConstraintOptionsDictionary()
         bc_list.append(bc)
@@ -1190,11 +1192,13 @@ def add_path_constraint(self, name, constraint_name=None, units=None, shape=None
         Parameters
         ----------
         name : str
-            Name of the response variable in the system.
+            Name of the variable to constrain. May also provide an expression to be evaluated and constrained.
+            If a single variable and the name is not a state, control, or 'time',
+            then this is assumed to be the path of the variable to be constrained in the ODE.
+            If an expression, it must be provided in the form of an equation with a left- and right-hand side
         constraint_name : str or None
-            The name of the variable as provided to the boundary constraint comp.  By
-            default this is the last element in `name` when split by dots.  The user may
-            override the constraint name if splitting the path causes name collisions.
+            The name of the path constraint. By default, this is 'var_constraint' if name is a single variable,
+             or the left-hand side of the equation if name is an expression.
         units : str or None
             The units in which the boundary constraint is to be applied.  If None, use the
             units associated with the constrained output.  If provided, must be compatible with
@@ -1235,9 +1239,9 @@ def add_path_constraint(self, name, constraint_name=None, units=None, shape=None
         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)"')
+            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
 
@@ -1253,6 +1257,14 @@ def add_path_constraint(self, name, constraint_name=None, units=None, shape=None
             raise ValueError(f'Cannot add new path constraint for variable `{name}` and indices {indices}. '
                              f'One already exists.')
 
+        existing_bc_name = [pc for pc in self._path_constraints
+                            if pc['name'] == constraint_name and pc['indices'] == indices
+                            and pc['flat_indices'] == flat_indices]
+
+        if existing_bc_name:
+            raise ValueError(f'Cannot add new path constraint named `{constraint_name}` and indices {indices}.'
+                             f' The name `{constraint_name}` is already in use as a path constraint')
+
         pc = ConstraintOptionsDictionary()
         self._path_constraints.append(pc)
 
diff --git a/dymos/phase/test/test_add_boundary_constraint.py b/dymos/phase/test/test_add_boundary_constraint.py
new file mode 100644
index 000000000..4e59a5e10
--- /dev/null
+++ b/dymos/phase/test/test_add_boundary_constraint.py
@@ -0,0 +1,254 @@
+import unittest
+
+import numpy as np
+
+import openmdao.api as om
+from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
+from openmdao.utils.assert_utils import assert_near_equal
+
+import dymos as dm
+
+
+class BrachistochroneVectorStatesODE(om.ExplicitComponent):
+    def initialize(self):
+        self.options.declare('num_nodes', types=int)
+
+    def setup(self):
+        nn = self.options['num_nodes']
+
+        self.add_input('v', val=np.zeros(nn), desc='velocity', units='m/s')
+        self.add_input('g', val=9.80665 * np.ones(nn), desc='grav. acceleration', units='m/s/s')
+        self.add_input('theta', val=np.zeros(nn), desc='angle of wire', units='rad')
+
+        self.add_output('pos_dot', val=np.zeros((nn, 2)), desc='velocity components', units='m/s')
+        self.add_output('vdot', val=np.zeros(nn), desc='acceleration magnitude', units='m/s**2')
+        self.add_output('check', val=np.zeros(nn), desc='check solution: v/sin(theta) = constant', units='m/s')
+
+        # Setup partials
+        arange = np.arange(self.options['num_nodes'])
+
+        self.declare_partials(of='vdot', wrt='g', rows=arange, cols=arange)
+        self.declare_partials(of='vdot', wrt='theta', rows=arange, cols=arange)
+
+        rs = np.arange(2*nn, dtype=int)
+        cs = np.repeat(np.arange(nn, dtype=int), 2)
+
+        self.declare_partials('*', '*', method='cs')
+        self.declare_coloring(wrt='*', method='cs', show_sparsity=True)
+
+    def compute(self, inputs, outputs):
+        theta = inputs['theta']
+        cos_theta = np.cos(theta)
+        sin_theta = np.sin(theta)
+        g = inputs['g']
+        v = inputs['v']
+
+        outputs['vdot'] = g * cos_theta
+        outputs['pos_dot'][:, 0] = v * sin_theta
+        outputs['pos_dot'][:, 1] = -v * cos_theta
+
+        outputs['check'] = v / sin_theta
+
+
+@use_tempdirs
+class TestAddBoundaryConstraint(unittest.TestCase):
+
+    @require_pyoptsparse(optimizer='SLSQP')
+    def test_simple_no_exception(self):
+        p = om.Problem(model=om.Group())
+
+        p.driver = om.pyOptSparseDriver()
+        p.driver.options['optimizer'] = 'SLSQP'
+
+        p.driver.declare_coloring()
+
+        transcription = dm.GaussLobatto(num_segments=3,
+                                        order=3,
+                                        compressed=True)
+
+        traj = dm.Trajectory()
+        phase = dm.Phase(ode_class=BrachistochroneVectorStatesODE,
+                         transcription=transcription)
+        traj.add_phase('phase0', phase)
+
+        p.model.add_subsystem('traj0', traj)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10), units='s')
+
+        # can't fix final position if you're solving the segments
+        phase.add_state('pos',
+                        rate_source='pos_dot', units='m',
+                        fix_initial=True)
+
+        # test add_boundary_constraint with arrays:
+        expected = np.array([10, 5])
+        phase.add_boundary_constraint(name='pos', loc='final', equals=expected)
+
+        phase.add_state('v',
+                        rate_source='vdot', units='m/s',
+                        fix_initial=True, fix_final=False)
+
+        phase.add_control('theta',
+                          continuity=True, rate_continuity=True,
+                          units='deg', lower=0.01, upper=179.9)
+
+        phase.add_parameter('g', units='m/s**2', val=9.80665, opt=False)
+
+        # Minimize time at the end of the phase
+        phase.add_objective('time', loc='final', scaler=10)
+
+        p.model.linear_solver = om.DirectSolver()
+        p.setup(check=True, force_alloc_complex=True)
+        p.final_setup()
+
+        p['traj0.phase0.t_initial'] = 0.0
+        p['traj0.phase0.t_duration'] = 1.8016
+
+        pos0 = [0, 10]
+        posf = [10, 5]
+
+        p['traj0.phase0.states:pos'] = phase.interp('pos', [pos0, posf])
+        p['traj0.phase0.states:v'] = phase.interp('v', [0, 9.9])
+        p['traj0.phase0.controls:theta'] = phase.interp('theta', [5, 100])
+        p['traj0.phase0.parameters:g'] = 9.80665
+
+        p.run_driver()
+
+        assert_near_equal(p.get_val('traj0.phase0.timeseries.states:pos')[-1, ...],
+                          [10, 5], tolerance=1.0E-5)
+
+    def test_invalid_expression(self):
+        p = om.Problem(model=om.Group())
+
+        p.driver = om.pyOptSparseDriver()
+        p.driver.options['optimizer'] = 'SLSQP'
+
+        p.driver.declare_coloring()
+
+        transcription = dm.GaussLobatto(num_segments=3,
+                                        order=3,
+                                        compressed=True)
+
+        traj = dm.Trajectory()
+        phase = dm.Phase(ode_class=BrachistochroneVectorStatesODE,
+                         transcription=transcription)
+        traj.add_phase('phase0', phase)
+
+        p.model.add_subsystem('traj0', traj)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10), units='s')
+
+        # can't fix final position if you're solving the segments
+        phase.add_state('pos',
+                        rate_source='pos_dot', units='m',
+                        fix_initial=True)
+
+        phase.add_state('v',
+                        rate_source='vdot', units='m/s',
+                        fix_initial=True, fix_final=False)
+
+        phase.add_control('theta',
+                          continuity=True, rate_continuity=True,
+                          units='deg', lower=0.01, upper=179.9)
+
+        phase.add_parameter('g', units='m/s**2', val=9.80665, opt=False)
+
+        # test add_boundary_constraint with arrays:
+        expected = 'The expression provided `pos**2` has invalid format. ' \
+                   'Expression may be a single variable or an equation ' \
+                   'of the form `constraint_name = func(vars)`'
+        with self.assertRaises(ValueError) as e:
+            phase.add_boundary_constraint(name='pos**2', loc='final', equals=np.array([10, 5]))
+
+        self.assertEqual(expected, str(e.exception))
+
+    def test_duplicate_name(self):
+        p = om.Problem(model=om.Group())
+
+        p.driver = om.pyOptSparseDriver()
+        p.driver.options['optimizer'] = 'SLSQP'
+
+        p.driver.declare_coloring()
+
+        transcription = dm.GaussLobatto(num_segments=3,
+                                        order=3,
+                                        compressed=True)
+
+        traj = dm.Trajectory()
+        phase = dm.Phase(ode_class=BrachistochroneVectorStatesODE,
+                         transcription=transcription)
+        traj.add_phase('phase0', phase)
+
+        p.model.add_subsystem('traj0', traj)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10), units='s')
+
+        # can't fix final position if you're solving the segments
+        phase.add_state('pos',
+                        rate_source='pos_dot', units='m',
+                        fix_initial=True)
+
+        phase.add_state('v',
+                        rate_source='vdot', units='m/s',
+                        fix_initial=True, fix_final=False)
+
+        phase.add_control('theta',
+                          continuity=True, rate_continuity=True,
+                          units='deg', lower=0.01, upper=179.9)
+
+        phase.add_parameter('g', units='m/s**2', val=9.80665, opt=False)
+
+        phase.add_boundary_constraint(name='pos', loc='final', equals=np.array([10, 5]))
+
+        # test add_boundary_constraint with arrays:
+        with self.assertRaises(ValueError) as e:
+            phase.add_boundary_constraint(name='pos=v**2', loc='final', equals=np.array([10, 5]))
+
+        expected = 'Cannot add new final boundary constraint named `pos` and indices None.' \
+                   ' The name `pos` is already in use as a final boundary constraint'
+        self.assertEqual(expected, str(e.exception))
+
+    def test_duplicate_constraint(self):
+        p = om.Problem(model=om.Group())
+
+        p.driver = om.pyOptSparseDriver()
+        p.driver.options['optimizer'] = 'SLSQP'
+
+        p.driver.declare_coloring()
+
+        transcription = dm.GaussLobatto(num_segments=3,
+                                        order=3,
+                                        compressed=True)
+
+        traj = dm.Trajectory()
+        phase = dm.Phase(ode_class=BrachistochroneVectorStatesODE,
+                         transcription=transcription)
+        traj.add_phase('phase0', phase)
+
+        p.model.add_subsystem('traj0', traj)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10), units='s')
+
+        # can't fix final position if you're solving the segments
+        phase.add_state('pos',
+                        rate_source='pos_dot', units='m',
+                        fix_initial=True)
+
+        phase.add_state('v',
+                        rate_source='vdot', units='m/s',
+                        fix_initial=True, fix_final=False)
+
+        phase.add_control('theta',
+                          continuity=True, rate_continuity=True,
+                          units='deg', lower=0.01, upper=179.9)
+
+        phase.add_parameter('g', units='m/s**2', val=9.80665, opt=False)
+
+        phase.add_boundary_constraint(name='pos', loc='final', equals=np.array([10, 5]))
+
+        # test add_boundary_constraint with arrays:
+        with self.assertRaises(ValueError) as e:
+            phase.add_boundary_constraint(name='pos', loc='final', equals=np.array([10, 5]))
+
+        expected = 'Cannot add new final boundary constraint for variable `pos` and indices None. One already exists.'
+        self.assertEqual(expected, str(e.exception))
diff --git a/dymos/phase/test/test_add_boundary_constraint_array.py b/dymos/phase/test/test_add_boundary_constraint_array.py
deleted file mode 100644
index b524c4f06..000000000
--- a/dymos/phase/test/test_add_boundary_constraint_array.py
+++ /dev/null
@@ -1,118 +0,0 @@
-import unittest
-
-import numpy as np
-
-import openmdao.api as om
-from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
-from openmdao.utils.assert_utils import assert_near_equal
-
-import dymos as dm
-
-
-class BrachistochroneVectorStatesODE(om.ExplicitComponent):
-    def initialize(self):
-        self.options.declare('num_nodes', types=int)
-
-    def setup(self):
-        nn = self.options['num_nodes']
-
-        self.add_input('v', val=np.zeros(nn), desc='velocity', units='m/s')
-        self.add_input('g', val=9.80665 * np.ones(nn), desc='grav. acceleration', units='m/s/s')
-        self.add_input('theta', val=np.zeros(nn), desc='angle of wire', units='rad')
-
-        self.add_output('pos_dot', val=np.zeros((nn, 2)), desc='velocity components', units='m/s')
-        self.add_output('vdot', val=np.zeros(nn), desc='acceleration magnitude', units='m/s**2')
-        self.add_output('check', val=np.zeros(nn), desc='check solution: v/sin(theta) = constant', units='m/s')
-
-        # Setup partials
-        arange = np.arange(self.options['num_nodes'])
-
-        self.declare_partials(of='vdot', wrt='g', rows=arange, cols=arange)
-        self.declare_partials(of='vdot', wrt='theta', rows=arange, cols=arange)
-
-        rs = np.arange(2*nn, dtype=int)
-        cs = np.repeat(np.arange(nn, dtype=int), 2)
-
-        self.declare_partials('*', '*', method='cs')
-        self.declare_coloring(wrt='*', method='cs', show_sparsity=True)
-
-    def compute(self, inputs, outputs):
-        theta = inputs['theta']
-        cos_theta = np.cos(theta)
-        sin_theta = np.sin(theta)
-        g = inputs['g']
-        v = inputs['v']
-
-        outputs['vdot'] = g * cos_theta
-        outputs['pos_dot'][:, 0] = v * sin_theta
-        outputs['pos_dot'][:, 1] = -v * cos_theta
-
-        outputs['check'] = v / sin_theta
-
-
-@use_tempdirs
-class TestAddBoundaryConstraint(unittest.TestCase):
-
-    @require_pyoptsparse(optimizer='SLSQP')
-    def test_simple_no_exception(self):
-        p = om.Problem(model=om.Group())
-
-        p.driver = om.pyOptSparseDriver()
-        p.driver.options['optimizer'] = 'SLSQP'
-
-        p.driver.declare_coloring()
-
-        transcription = dm.GaussLobatto(num_segments=3,
-                                        order=3,
-                                        compressed=True)
-
-        traj = dm.Trajectory()
-        phase = dm.Phase(ode_class=BrachistochroneVectorStatesODE,
-                         transcription=transcription)
-        traj.add_phase('phase0', phase)
-
-        p.model.add_subsystem('traj0', traj)
-
-        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10), units='s')
-
-        # can't fix final position if you're solving the segments
-        phase.add_state('pos',
-                        rate_source='pos_dot', units='m',
-                        fix_initial=True)
-
-        # test add_boundary_constraint with arrays:
-        expected = np.array([10, 5])
-        phase.add_boundary_constraint(name='pos', loc='final', equals=expected)
-
-        phase.add_state('v',
-                        rate_source='vdot', units='m/s',
-                        fix_initial=True, fix_final=False)
-
-        phase.add_control('theta',
-                          continuity=True, rate_continuity=True,
-                          units='deg', lower=0.01, upper=179.9)
-
-        phase.add_parameter('g', units='m/s**2', val=9.80665, opt=False)
-
-        # Minimize time at the end of the phase
-        phase.add_objective('time', loc='final', scaler=10)
-
-        p.model.linear_solver = om.DirectSolver()
-        p.setup(check=True, force_alloc_complex=True)
-        p.final_setup()
-
-        p['traj0.phase0.t_initial'] = 0.0
-        p['traj0.phase0.t_duration'] = 1.8016
-
-        pos0 = [0, 10]
-        posf = [10, 5]
-
-        p['traj0.phase0.states:pos'] = phase.interp('pos', [pos0, posf])
-        p['traj0.phase0.states:v'] = phase.interp('v', [0, 9.9])
-        p['traj0.phase0.controls:theta'] = phase.interp('theta', [5, 100])
-        p['traj0.phase0.parameters:g'] = 9.80665
-
-        p.run_driver()
-
-        assert_near_equal(p.get_val('traj0.phase0.timeseries.states:pos')[-1, ...],
-                          [10, 5], tolerance=1.0E-5)
diff --git a/dymos/phase/test/test_add_path_constraint.py b/dymos/phase/test/test_add_path_constraint.py
new file mode 100644
index 000000000..ed9847e08
--- /dev/null
+++ b/dymos/phase/test/test_add_path_constraint.py
@@ -0,0 +1,190 @@
+import unittest
+
+import numpy as np
+
+import openmdao.api as om
+from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
+from openmdao.utils.assert_utils import assert_near_equal
+
+import dymos as dm
+
+
+class BrachistochroneVectorStatesODE(om.ExplicitComponent):
+    def initialize(self):
+        self.options.declare('num_nodes', types=int)
+
+    def setup(self):
+        nn = self.options['num_nodes']
+
+        self.add_input('v', val=np.zeros(nn), desc='velocity', units='m/s')
+        self.add_input('g', val=9.80665 * np.ones(nn), desc='grav. acceleration', units='m/s/s')
+        self.add_input('theta', val=np.zeros(nn), desc='angle of wire', units='rad')
+
+        self.add_output('pos_dot', val=np.zeros((nn, 2)), desc='velocity components', units='m/s')
+        self.add_output('vdot', val=np.zeros(nn), desc='acceleration magnitude', units='m/s**2')
+        self.add_output('check', val=np.zeros(nn), desc='check solution: v/sin(theta) = constant', units='m/s')
+
+        # Setup partials
+        arange = np.arange(self.options['num_nodes'])
+
+        self.declare_partials(of='vdot', wrt='g', rows=arange, cols=arange)
+        self.declare_partials(of='vdot', wrt='theta', rows=arange, cols=arange)
+
+        rs = np.arange(2*nn, dtype=int)
+        cs = np.repeat(np.arange(nn, dtype=int), 2)
+
+        self.declare_partials('*', '*', method='cs')
+        self.declare_coloring(wrt='*', method='cs', show_sparsity=True)
+
+    def compute(self, inputs, outputs):
+        theta = inputs['theta']
+        cos_theta = np.cos(theta)
+        sin_theta = np.sin(theta)
+        g = inputs['g']
+        v = inputs['v']
+
+        outputs['vdot'] = g * cos_theta
+        outputs['pos_dot'][:, 0] = v * sin_theta
+        outputs['pos_dot'][:, 1] = -v * cos_theta
+
+        outputs['check'] = v / sin_theta
+
+
+@use_tempdirs
+class TestAddPathConstraint(unittest.TestCase):
+
+    def test_invalid_expression(self):
+        p = om.Problem(model=om.Group())
+
+        p.driver = om.pyOptSparseDriver()
+        p.driver.options['optimizer'] = 'SLSQP'
+
+        p.driver.declare_coloring()
+
+        transcription = dm.GaussLobatto(num_segments=3,
+                                        order=3,
+                                        compressed=True)
+
+        traj = dm.Trajectory()
+        phase = dm.Phase(ode_class=BrachistochroneVectorStatesODE,
+                         transcription=transcription)
+        traj.add_phase('phase0', phase)
+
+        p.model.add_subsystem('traj0', traj)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10), units='s')
+
+        # can't fix final position if you're solving the segments
+        phase.add_state('pos',
+                        rate_source='pos_dot', units='m',
+                        fix_initial=True)
+
+        phase.add_state('v',
+                        rate_source='vdot', units='m/s',
+                        fix_initial=True, fix_final=False)
+
+        phase.add_control('theta',
+                          continuity=True, rate_continuity=True,
+                          units='deg', lower=0.01, upper=179.9)
+
+        phase.add_parameter('g', units='m/s**2', val=9.80665, opt=False)
+
+        # test add_boundary_constraint with arrays:
+        expected = 'The expression provided `pos**2` has invalid format. ' \
+                   'Expression may be a single variable or an equation ' \
+                   'of the form `constraint_name = func(vars)`'
+        with self.assertRaises(ValueError) as e:
+            phase.add_path_constraint(name='pos**2', equals=np.array([10, 5]))
+
+        self.assertEqual(expected, str(e.exception))
+
+    def test_duplicate_name(self):
+        p = om.Problem(model=om.Group())
+
+        p.driver = om.pyOptSparseDriver()
+        p.driver.options['optimizer'] = 'SLSQP'
+
+        p.driver.declare_coloring()
+
+        transcription = dm.GaussLobatto(num_segments=3,
+                                        order=3,
+                                        compressed=True)
+
+        traj = dm.Trajectory()
+        phase = dm.Phase(ode_class=BrachistochroneVectorStatesODE,
+                         transcription=transcription)
+        traj.add_phase('phase0', phase)
+
+        p.model.add_subsystem('traj0', traj)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10), units='s')
+
+        # can't fix final position if you're solving the segments
+        phase.add_state('pos',
+                        rate_source='pos_dot', units='m',
+                        fix_initial=True)
+
+        phase.add_state('v',
+                        rate_source='vdot', units='m/s',
+                        fix_initial=True, fix_final=False)
+
+        phase.add_control('theta',
+                          continuity=True, rate_continuity=True,
+                          units='deg', lower=0.01, upper=179.9)
+
+        phase.add_parameter('g', units='m/s**2', val=9.80665, opt=False)
+
+        phase.add_path_constraint(name='pos', equals=np.array([10, 5]))
+
+        # test add_boundary_constraint with arrays:
+        with self.assertRaises(ValueError) as e:
+            phase.add_path_constraint(name='pos=v**2', equals=np.array([10, 5]))
+
+        expected = 'Cannot add new path constraint named `pos` and indices None.' \
+                   ' The name `pos` is already in use as a path constraint'
+        self.assertEqual(expected, str(e.exception))
+
+    def test_duplicate_constraint(self):
+        p = om.Problem(model=om.Group())
+
+        p.driver = om.pyOptSparseDriver()
+        p.driver.options['optimizer'] = 'SLSQP'
+
+        p.driver.declare_coloring()
+
+        transcription = dm.GaussLobatto(num_segments=3,
+                                        order=3,
+                                        compressed=True)
+
+        traj = dm.Trajectory()
+        phase = dm.Phase(ode_class=BrachistochroneVectorStatesODE,
+                         transcription=transcription)
+        traj.add_phase('phase0', phase)
+
+        p.model.add_subsystem('traj0', traj)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10), units='s')
+
+        # can't fix final position if you're solving the segments
+        phase.add_state('pos',
+                        rate_source='pos_dot', units='m',
+                        fix_initial=True)
+
+        phase.add_state('v',
+                        rate_source='vdot', units='m/s',
+                        fix_initial=True, fix_final=False)
+
+        phase.add_control('theta',
+                          continuity=True, rate_continuity=True,
+                          units='deg', lower=0.01, upper=179.9)
+
+        phase.add_parameter('g', units='m/s**2', val=9.80665, opt=False)
+
+        phase.add_path_constraint(name='pos', equals=np.array([10, 5]))
+
+        # test add_boundary_constraint with arrays:
+        with self.assertRaises(ValueError) as e:
+            phase.add_path_constraint(name='pos', equals=np.array([10, 5]))
+
+        expected = 'Cannot add new path constraint for variable `pos` and indices None. One already exists.'
+        self.assertEqual(expected, str(e.exception))

From e354af2a229a09608c79c7af18249deb5b2c3e08 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Tue, 6 Dec 2022 15:21:36 -0600
Subject: [PATCH 32/34] Fixed docstrings and formatting

---
 dymos/phase/phase.py                             | 8 ++++----
 dymos/phase/test/test_add_boundary_constraint.py | 6 +++---
 dymos/phase/test/test_add_path_constraint.py     | 6 +++---
 3 files changed, 10 insertions(+), 10 deletions(-)

diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index a81e2f8e8..a7613e161 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -1081,7 +1081,7 @@ def add_boundary_constraint(self, name, loc, constraint_name=None, units=None,
             Name of the variable to constrain. May also provide an expression to be evaluated and constrained.
             If a single variable and the name is not a state, control, or 'time',
             then this is assumed to be the path of the variable to be constrained in the ODE.
-            If an expression, it must be provided in the form of an equation with a left- and right-hand side
+            If an expression, it must be provided in the form of an equation with a left- and right-hand side.
         loc : str
             The location of the boundary constraint ('initial' or 'final').
         constraint_name : str or None
@@ -1195,7 +1195,7 @@ def add_path_constraint(self, name, constraint_name=None, units=None, shape=None
             Name of the variable to constrain. May also provide an expression to be evaluated and constrained.
             If a single variable and the name is not a state, control, or 'time',
             then this is assumed to be the path of the variable to be constrained in the ODE.
-            If an expression, it must be provided in the form of an equation with a left- and right-hand side
+            If an expression, it must be provided in the form of an equation with a left- and right-hand side.
         constraint_name : str or None
             The name of the path constraint. By default, this is 'var_constraint' if name is a single variable,
              or the left-hand side of the equation if name is an expression.
@@ -1258,8 +1258,8 @@ def add_path_constraint(self, name, constraint_name=None, units=None, shape=None
                              f'One already exists.')
 
         existing_bc_name = [pc for pc in self._path_constraints
-                            if pc['name'] == constraint_name and pc['indices'] == indices
-                            and pc['flat_indices'] == flat_indices]
+                            if pc['name'] == constraint_name and pc['indices'] == indices and
+                            pc['flat_indices'] == flat_indices]
 
         if existing_bc_name:
             raise ValueError(f'Cannot add new path constraint named `{constraint_name}` and indices {indices}.'
diff --git a/dymos/phase/test/test_add_boundary_constraint.py b/dymos/phase/test/test_add_boundary_constraint.py
index 4e59a5e10..1876b71e9 100644
--- a/dymos/phase/test/test_add_boundary_constraint.py
+++ b/dymos/phase/test/test_add_boundary_constraint.py
@@ -57,7 +57,7 @@ class TestAddBoundaryConstraint(unittest.TestCase):
     def test_simple_no_exception(self):
         p = om.Problem(model=om.Group())
 
-        p.driver = om.pyOptSparseDriver()
+        p.driver = om.ScipyOptimizeDriver()
         p.driver.options['optimizer'] = 'SLSQP'
 
         p.driver.declare_coloring()
@@ -120,7 +120,7 @@ def test_simple_no_exception(self):
     def test_invalid_expression(self):
         p = om.Problem(model=om.Group())
 
-        p.driver = om.pyOptSparseDriver()
+        p.driver = om.ScipyOptimizeDriver()
         p.driver.options['optimizer'] = 'SLSQP'
 
         p.driver.declare_coloring()
@@ -165,7 +165,7 @@ def test_invalid_expression(self):
     def test_duplicate_name(self):
         p = om.Problem(model=om.Group())
 
-        p.driver = om.pyOptSparseDriver()
+        p.driver = om.ScipyOptimizeDriver()
         p.driver.options['optimizer'] = 'SLSQP'
 
         p.driver.declare_coloring()
diff --git a/dymos/phase/test/test_add_path_constraint.py b/dymos/phase/test/test_add_path_constraint.py
index ed9847e08..83f2f8622 100644
--- a/dymos/phase/test/test_add_path_constraint.py
+++ b/dymos/phase/test/test_add_path_constraint.py
@@ -56,7 +56,7 @@ class TestAddPathConstraint(unittest.TestCase):
     def test_invalid_expression(self):
         p = om.Problem(model=om.Group())
 
-        p.driver = om.pyOptSparseDriver()
+        p.driver = om.ScipyOptimizeDriver()
         p.driver.options['optimizer'] = 'SLSQP'
 
         p.driver.declare_coloring()
@@ -101,7 +101,7 @@ def test_invalid_expression(self):
     def test_duplicate_name(self):
         p = om.Problem(model=om.Group())
 
-        p.driver = om.pyOptSparseDriver()
+        p.driver = om.ScipyOptimizeDriver()
         p.driver.options['optimizer'] = 'SLSQP'
 
         p.driver.declare_coloring()
@@ -147,7 +147,7 @@ def test_duplicate_name(self):
     def test_duplicate_constraint(self):
         p = om.Problem(model=om.Group())
 
-        p.driver = om.pyOptSparseDriver()
+        p.driver = om.ScipyOptimizeDriver()
         p.driver.options['optimizer'] = 'SLSQP'
 
         p.driver.declare_coloring()

From 2fabf73d784c8ad5f5cd19c3cce60aaf544a0fd5 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Wed, 7 Dec 2022 07:54:01 -0600
Subject: [PATCH 33/34] Changed tests to use scipy driver

---
 dymos/phase/test/test_add_boundary_constraint.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/dymos/phase/test/test_add_boundary_constraint.py b/dymos/phase/test/test_add_boundary_constraint.py
index 1876b71e9..03d631491 100644
--- a/dymos/phase/test/test_add_boundary_constraint.py
+++ b/dymos/phase/test/test_add_boundary_constraint.py
@@ -211,7 +211,7 @@ def test_duplicate_name(self):
     def test_duplicate_constraint(self):
         p = om.Problem(model=om.Group())
 
-        p.driver = om.pyOptSparseDriver()
+        p.driver = om.ScipyOptimizeDriver()
         p.driver.options['optimizer'] = 'SLSQP'
 
         p.driver.declare_coloring()

From 255d6aa84f9c8c4a62e69bb65f33d909e04160d9 Mon Sep 17 00:00:00 2001
From: Kaushik Ponnapalli <kaushik.s.ponnapalli@nasa.gov>
Date: Wed, 7 Dec 2022 11:22:11 -0600
Subject: [PATCH 34/34] Added documentation

---
 .../features/phases/constraints.ipynb         | 19 ++++++++++++++++++-
 1 file changed, 18 insertions(+), 1 deletion(-)

diff --git a/docs/dymos_book/features/phases/constraints.ipynb b/docs/dymos_book/features/phases/constraints.ipynb
index f1c7a0dfa..fb5ce38d9 100644
--- a/docs/dymos_book/features/phases/constraints.ipynb
+++ b/docs/dymos_book/features/phases/constraints.ipynb
@@ -143,6 +143,23 @@
     "\n",
     "The downside of path constraints is that they add a considerable number of constraints to the NLP problem and thus may negatively impact performance, although this is generally minor for many problems.\n",
     "\n",
+    "##  Constraining Expressions\n",
+    "\n",
+    "Constraints may be defined using mathematical expressions of the form `y=f(x)` to be evaluated. Here `x` may be vector combination of time, states, controls, parameters, or any outputs of the ODE. The variable `y` is added to the timeseries outputs of the phase and the desired constraint is applied to it.\n",
+    "\n",
+    "Consider, again, the example of maximizing the range flown by a cannonball. But now, rather than a constraint on the initial velocity, we wish to apply a constraint to the initial normalized kinetic energy.\n",
+    "\n",
+    "\\begin{align}\n",
+    "    ke &= 0.5 * v^2 \\\\\n",
+    "    ke_0 &= 5000 \\, \\frac{\\mathrm{m^2}}{\\mathrm{s^2}}\n",
+    "\\end{align}\n",
+    "\n",
+    "The first way to achieve this is to add kinetic energy as a state in the model. This state may then be constrained either using `fix_initial=True` or `add_boundary_constraint(‘ke’, loc=’initial’ , equals=5000)`.\n",
+    "\n",
+    "The second method to add this constraint is to add the constraint using an expression. This may be done as `add_boundary_constraint(‘ke=0.5*v**2, loc=’initial’, equals=5000)`.\n",
+    "\n",
+    "The advantage to the latter method is that no changes are required to the previous model, simply modifying the `add_boundary_constraint` statement is sufficient. The disadvantage is that the derivatives of the constraints are evaluated using complex step rather than analytical expressions. This may negatively impact performance, but the effect should be minor for most problems.\n",
+    "\n",
     "## Implementation Detail - Constraint Aliasing\n",
     "\n",
     "As of OpenMDAO Version 3.17.0, multiple constraints are allowed to exist on the same output as long as they use different indices and are provided different aliases.\n",
@@ -157,7 +174,7 @@
     "\n",
     "The use of the `->` in this case is intended to remind the user that this is not the actual path to the variable being constrained.\n",
     "\n",
-    "## Constraint Linearity\n",
+    "##  Constraint Linearity\n",
     "\n",
     "OpenMDAO will treat all boundary and path consraints as nonlinear unless the user provides the argument `linear=True` to `add_boundary_constraint` or `add_path_constraint`.\n",
     "Note that it is the responsibility of the user to understand the inner workings of Dymos and their model well enough to know if the constraint may be treated as linear.\n",