Repository: OpenMDAO/dymos
Pull request #781 (2022-07-12 20:00:28)

Added timeseries rates for pseudospectral transcriptions

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

Coverage decreased (-0.06%) to 95.461% when pulling **270a14af89d45e67961bf138adc287babf793cfc on naylor-b:timeseries_rates** into **bb1cadc80a2cc1c2cdf2d922c8b69ea011d8bead on OpenMDAO:master**.


Patch:
diff --git a/docs/dymos_book/faq/downstream_analysis.ipynb b/docs/dymos_book/faq/downstream_analysis.ipynb
index ca515d5a..7c0dae32 100644
--- a/docs/dymos_book/faq/downstream_analysis.ipynb
+++ b/docs/dymos_book/faq/downstream_analysis.ipynb
@@ -47,7 +47,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# How do connect the outputs of Dymos to a downstream analysis?\n",
+    "# How do I connect the outputs of Dymos to a downstream analysis?\n",
     "\n",
     "One of the design goals of Dymos is to allow the trajectory to be a part of a larger multidisciplinary optimization problem.\n",
     "Sometimes, you may want to take the results from the Dymos trajectory and feed them to some downstream analysis.\n",
diff --git a/docs/dymos_book/faq/upstream_analysis.ipynb b/docs/dymos_book/faq/upstream_analysis.ipynb
index 9a50d954..204fe61c 100644
--- a/docs/dymos_book/faq/upstream_analysis.ipynb
+++ b/docs/dymos_book/faq/upstream_analysis.ipynb
@@ -47,7 +47,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "# How do connect the outputs of an upstream analysis as inputs to Dymos?\n",
+    "# How do I connect the outputs of an upstream analysis as inputs to Dymos?\n",
     "\n",
     "One of the design goals of Dymos is to allow the trajectory to be a part of a larger multidisciplinary optimization problem.\n",
     "Sometimes, you may want to take the results of some upstream design and use them as an input to a Dymos trajectory.\n",
diff --git a/docs/dymos_book/features/trajectories/trajectories.ipynb b/docs/dymos_book/features/trajectories/trajectories.ipynb
index 73add00b..5fd96cda 100644
--- a/docs/dymos_book/features/trajectories/trajectories.ipynb
+++ b/docs/dymos_book/features/trajectories/trajectories.ipynb
@@ -224,7 +224,7 @@
     "If we implemented these as parameters within each phase individually, we would need some constraints to ensure that they held the same value within each phase.\n",
     "To avoid this complexity, Dymos Trajectory objects support their own Parameters.\n",
     "\n",
-    "Like their Phase-based counterparts, Trajectory parameters produce may be design variables for the problem or used as inputs to the trajectory from external sources.\n",
+    "Like their Phase-based counterparts, Trajectory parameters may be design variables for the problem or they may be inputs to the trajectory from external sources.\n",
     "\n",
     "When using Trajectory parameters, their values are connected to each phase as an Input Parameter within the Phase.\n",
     "Because ODEs in different phases may have different names for parameters (e.g. 'mass', 'm', 'm_total', etc) Dymos allows the user to specify the targeted ODE parameters on a phase-by-phase basis using the `targets` and `target_params` option.\n",
diff --git a/docs/dymos_book/getting_started/intro_to_dymos/intro_segments.ipynb b/docs/dymos_book/getting_started/intro_to_dymos/intro_segments.ipynb
index bdf9ac31..d1c35b61 100644
--- a/docs/dymos_book/getting_started/intro_to_dymos/intro_segments.ipynb
+++ b/docs/dymos_book/getting_started/intro_to_dymos/intro_segments.ipynb
@@ -156,7 +156,7 @@
     "\n",
     "## How does the number and order of segments affect the solution?\n",
     "\n",
-    "Obviously, a single third-order polynomial won't be able to fit highly oscillator behavior.\n",
+    "Obviously, a single third-order polynomial won't be able to fit highly oscillatory behavior.\n",
     "In this case, our guess of using four segments (equally spaced) in the phase wasn't quite sufficient.\n",
     "Let's try increasing that number to ten third-order segments."
    ]
diff --git a/dymos/phase/options.py b/dymos/phase/options.py
index 7efab2bc..a7cbd373 100644
--- a/dymos/phase/options.py
+++ b/dymos/phase/options.py
@@ -685,3 +685,6 @@ class TimeseriesOutputOptionsDictionary(om.OptionsDictionary):
 
         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='is_rate', default=False, allow_none=False,
+                     desc='If True this is a rate.')
diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index 1f622b85..86f1eb5a 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -75,7 +75,7 @@ class Phase(om.Group):
             self._path_constraints = []
             self._timeseries = {'timeseries': {'transcription': None,
                                                'subset': 'all',
-                                               'outputs': []}}
+                                               'outputs': {}}}
             self._objectives = {}
         else:
             self.time_options.update(from_phase.time_options)
@@ -613,7 +613,7 @@ class Phase(om.Group):
         name : str
             Name of the controllable parameter in the ODE.
         order : int
-            The order of the interpolating polynomial used to represent the control valeu in
+            The order of the interpolating polynomial used to represent the control value in
             phase tau space.
         desc : str
             A description of the polynomial control.
@@ -1050,8 +1050,7 @@ class Phase(om.Group):
         # Automatically add the requested variable to the timeseries outputs if it's an ODE output.
         var_type = self.classify_var(name)
         if var_type == 'ode':
-            current_output_names = [op['output_name'] for op in self._timeseries['timeseries']['outputs']]
-            if constraint_name not in current_output_names:
+            if constraint_name not in self._timeseries['timeseries']['outputs']:
                 self.add_timeseries_output(name, output_name=constraint_name, units=units, shape=shape)
 
     def add_path_constraint(self, name, constraint_name=None, units=None, shape=None, indices=None,
@@ -1135,8 +1134,7 @@ class Phase(om.Group):
         # Automatically add the requested variable to the timeseries outputs if it's an ODE output.
         var_type = self.classify_var(name)
         if var_type == 'ode':
-            current_output_names = [op['output_name'] for op in self._timeseries['timeseries']['outputs']]
-            if constraint_name not in current_output_names:
+            if constraint_name not in self._timeseries['timeseries']['outputs']:
                 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,
@@ -1149,7 +1147,61 @@ class Phase(om.Group):
         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, or the path to an output variable in the ODE.
+            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
+            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.
+        units : str or None or _unspecified
+            The units to express the timeseries output.  If None, use the
+            units associated with the target.  If provided, must be compatible with
+            the target units.
+            If a list of names is provided, units can be a matching list or dictionary.
+        shape : tuple or _unspecified
+            The shape of the timeseries output variable.  This must be provided (if not scalar)
+            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.
+        """
+        if type(name) is list:
+            for i, name_i in enumerate(name):
+                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
+                    unit = units[i]
+                else:
+                    unit = units
+
+                oname = self._add_timeseries_output(name_i, output_name=output_name,
+                                                    units=unit,
+                                                    shape=shape,
+                                                    timeseries=timeseries,
+                                                    rate=False)
+
+                # Handle specific units for wildcard names.
+                if oname is not None and '*' in name_i:
+                    self._timeseries[timeseries]['outputs'][oname]['wildcard_units'] = units
+
+        else:
+            self._add_timeseries_output(name, output_name=output_name,
+                                        units=units,
+                                        shape=shape,
+                                        timeseries=timeseries,
+                                        rate=False)
+
+    def add_timeseries_rate_output(self, name, output_name=None, units=_unspecified, shape=_unspecified,
+                                   timeseries='timeseries'):
+        r"""
+        Add the rate of a variable to the timeseries outputs of the phase.
+
+        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
+            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
             default this is the last element in `name` when split by dots.  The user may
@@ -1174,27 +1226,30 @@ class Phase(om.Group):
                 else:
                     unit = units
 
-                self._add_timeseries_output(name_i, output_name=output_name,
-                                            units=unit,
-                                            shape=shape,
-                                            timeseries=timeseries)
+                oname = self._add_timeseries_output(name_i, output_name=output_name,
+                                                    units=unit,
+                                                    shape=shape,
+                                                    timeseries=timeseries,
+                                                    rate=True)
 
                 # Handle specific units for wildcard names.
-                if '*' in name_i:
-                    self._timeseries[timeseries]['outputs'][-1]['wildcard_units'] = units
+                if oname is not None and '*' in name_i:
+                    self._timeseries[timeseries]['outputs'][oname]['wildcard_units'] = units
 
         else:
             self._add_timeseries_output(name, output_name=output_name,
                                         units=units,
                                         shape=shape,
-                                        timeseries=timeseries)
+                                        timeseries=timeseries,
+                                        rate=True)
 
     def _add_timeseries_output(self, name, output_name=None, units=_unspecified, shape=_unspecified,
-                               timeseries='timeseries'):
+                               timeseries='timeseries', rate=False):
         r"""
-        Add a single variable to the timeseries outputs of the phase.
+        Add a single variable or rate to the timeseries outputs of the phase.
 
-        This is called by add_timeseries_output for each variable that is added.
+        This is called by add_timeseries_output or add_timeseries_rate_output for each variable or rate
+        that is added.
 
         Parameters
         ----------
@@ -1205,7 +1260,8 @@ class Phase(om.Group):
         output_name : str or None
             The name of the variable as listed in the phase timeseries outputs.  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.
+            override the constraint name if splitting the path causes name collisions.  If rate
+            is True, the rate name will be this name + _rate.
         units : str or None
             The units to express the timeseries output.  If None, use the
             units associated with the target.  If provided, must be compatible with
@@ -1215,27 +1271,39 @@ class Phase(om.Group):
             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.
-        """
-        if output_name is None:
-            output_name = name.split('.')[-1]
+        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.
 
+        Returns
+        -------
+        str or None
+           Name of output that was added to the timeseries or None if nothing was added.
+        """
         if timeseries not in self._timeseries:
             raise ValueError(f'Timeseries {timeseries} does not exist in phase {self.pathname}')
 
-        current_output_names = [op['output_name'] for op in self._timeseries[timeseries]['outputs']]
+        if output_name is None:
+            output_name = name.split('.')[-1]
+
+            if rate:
+                output_name = output_name + '_rate'
 
-        if output_name not in current_output_names:
+        if output_name in self._timeseries[timeseries]['outputs']:
+            om.issue_warning(f'Output name `{output_name}` is already in timeseries `{timeseries}`. '
+                             f'New output ignored.')
+        else:
             ts_output = TimeseriesOutputOptionsDictionary()
             ts_output['name'] = name
             ts_output['output_name'] = output_name
             ts_output['wildcard_units'] = {}
             ts_output['units'] = units
             ts_output['shape'] = shape
+            ts_output['is_rate'] = rate
 
-            self._timeseries[timeseries]['outputs'].append(ts_output)
-        else:
-            om.issue_warning(f'Output name `{output_name}` is already in timeseries `{timeseries}`. '
-                             f'New output ignored.')
+            self._timeseries[timeseries]['outputs'][output_name] = ts_output
+
+            return output_name
 
     def add_timeseries(self, name, transcription, subset='all'):
         r"""
@@ -1254,7 +1322,7 @@ class Phase(om.Group):
         """
         self._timeseries[name] = {'transcription': transcription,
                                   'subset': subset,
-                                  'outputs': []}
+                                  'outputs': {}}
 
     def add_objective(self, name, loc='final', index=None, shape=(1,), ref=None, ref0=None,
                       adder=None, scaler=None, parallel_deriv_color=None):
@@ -1585,13 +1653,11 @@ class Phase(om.Group):
         """
         transcription = self.options['transcription']
         state_options = self.state_options
-        out_meta = get_promoted_vars(transcription._get_ode(self), 'output')
+        out_meta = get_promoted_vars(transcription._get_ode(self), 'output', metadata_keys=('tags',))
 
-        for name, meta in out_meta.items():
-            tags = meta['tags']
-            prom_name = meta['prom_name']
+        for prom_name, meta in out_meta.items():
             state = None
-            for tag in sorted(tags):
+            for tag in sorted(meta['tags']):
 
                 # Declared as rate_source.
                 if tag.startswith('dymos.state_rate_source:') or tag.startswith('state_rate_source:'):
@@ -1607,7 +1673,7 @@ class Phase(om.Group):
                     if state_options[state]['rate_source'] is not None:
                         if state_options[state]['rate_source'] != prom_name:
                             raise ValueError(f"rate_source has been declared twice for state "
-                                             f"'{state}' which is tagged on '{name}'.")
+                                             f"'{state}' which is tagged on '{prom_name}'.")
 
                     state_options[state]['rate_source'] = prom_name
 
@@ -2275,7 +2341,6 @@ class Phase(om.Group):
         all_flat_idxs : set
             A C-order flattened set of indices that apply to the constraint.
         """
-        s = {'initial': 'initial boundary', 'final': 'final boundary', 'path': 'path'}
         cons = {'initial': self._initial_boundary_constraints,
                 'final': self._final_boundary_constraints,
                 'path': self._path_constraints}
@@ -2285,12 +2350,15 @@ class Phase(om.Group):
         for con in cons[loc]:
             if con['name'] != name:
                 continue
+
             flat_idxs = get_constraint_flat_idxs(con)
             duplicate_idxs = all_flat_idxs.intersection(flat_idxs)
             if duplicate_idxs:
+                s = {'initial': 'initial boundary', 'final': 'final boundary', 'path': 'path'}
                 raise ValueError(f'Duplicate constraint in phase {self.pathname}. '
                                  f'The following indices of `{name}` are used in '
                                  f'multiple {s[loc]} constraints:\n{duplicate_idxs}')
+
             all_flat_idxs.update(flat_idxs)
 
         return all_flat_idxs
diff --git a/dymos/transcriptions/common/timeseries_output_comp.py b/dymos/transcriptions/common/timeseries_output_comp.py
index 1afb2889..8a0fbdee 100644
--- a/dymos/transcriptions/common/timeseries_output_comp.py
+++ b/dymos/transcriptions/common/timeseries_output_comp.py
@@ -28,8 +28,6 @@ class TimeseriesOutputCompBase(om.ExplicitComponent):
         """
         Declare component options.
         """
-        self._timeseries_outputs = []
-
         self._vars = {}
 
         self.options.declare('input_grid_data',
diff --git a/dymos/transcriptions/explicit_shooting/explicit_shooting.py b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
index 804d0b85..eb5c462d 100644
--- a/dymos/transcriptions/explicit_shooting/explicit_shooting.py
+++ b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
@@ -570,7 +570,7 @@ class ExplicitShooting(TranscriptionBase):
             units = phase.state_options[var]['units']
             linear = loc == 'initial'
             obj_path = f'integrator.states_out:{var}'
-        elif var_type in 'indep_control':
+        elif var_type == 'indep_control':
             shape = phase.control_options[var]['shape']
             units = phase.control_options[var]['units']
             linear = True
@@ -580,7 +580,7 @@ class ExplicitShooting(TranscriptionBase):
             units = phase.control_options[var]['units']
             linear = False
             obj_path = f'control_values:{var}'
-        elif var_type in 'indep_polynomial_control':
+        elif var_type == 'indep_polynomial_control':
             shape = phase.polynomial_control_options[var]['shape']
             units = phase.polynomial_control_options[var]['units']
             linear = True
diff --git a/dymos/transcriptions/explicit_shooting/explicit_timeseries_comp.py b/dymos/transcriptions/explicit_shooting/explicit_timeseries_comp.py
index 666dd05c..778fc0d6 100644
--- a/dymos/transcriptions/explicit_shooting/explicit_timeseries_comp.py
+++ b/dymos/transcriptions/explicit_shooting/explicit_timeseries_comp.py
@@ -31,7 +31,6 @@ class ExplicitTimeseriesComp(TimeseriesOutputCompBase):
         # from another variable has potentially different units
         self._units = {}
         self._conversion_factors = {}
-        self._vars = {}
 
         self._no_check_partials = not dymos_options['include_check_partials']
 
@@ -50,7 +49,7 @@ class ExplicitTimeseriesComp(TimeseriesOutputCompBase):
         self.input_num_nodes = igd.subset_num_nodes['segment_ends']
         self.output_num_nodes = self.input_num_nodes
 
-    def _add_output_configure(self, name, units, shape, desc='', src=None):
+    def _add_output_configure(self, name, units, shape, desc='', src=None, rate=False):
         """
         Add a single timeseries output.
 
@@ -70,12 +69,18 @@ class ExplicitTimeseriesComp(TimeseriesOutputCompBase):
             description of the timeseries output variable.
         src : str
             The src path of the variables input, used to prevent redundant inputs.
+        rate : bool
+            If True, timeseries output is a rate.
 
         Returns
         -------
         bool
             True if a new input was added for the output, or False if it reuses an existing input.
         """
+        if rate:
+            raise NotImplementedError("Timeseries output rates are not currently supported for "
+                                      "ExplicitShooting transcriptions.")
+
         input_num_nodes = self.input_num_nodes
         output_num_nodes = self.output_num_nodes
         added_source = False
@@ -102,7 +107,7 @@ class ExplicitTimeseriesComp(TimeseriesOutputCompBase):
                         shape=(output_num_nodes,) + shape,
                         units=units, desc=desc)
 
-        self._vars[name] = (input_name, output_name, shape)
+        self._vars[name] = (input_name, output_name, shape, rate)
 
         size = np.prod(shape)
         rs = cs = np.arange(output_num_nodes * size, dtype=int)
@@ -114,7 +119,7 @@ class ExplicitTimeseriesComp(TimeseriesOutputCompBase):
             offset = 0
         else:
             scale, offset = unit_conversion(input_units, units)
-        self._conversion_factors[output_name] = scale, offset
+            self._conversion_factors[output_name] = scale, offset
 
         self.declare_partials(of=output_name,
                               wrt=input_name,
@@ -133,6 +138,9 @@ class ExplicitTimeseriesComp(TimeseriesOutputCompBase):
         outputs : `Vector`
             `Vector` containing outputs.
         """
-        for (input_name, output_name, _) in self._vars.values():
-            scale, offset = self._conversion_factors[output_name]
-            outputs[output_name] = scale * (inputs[input_name] + offset)
+        for (input_name, output_name, _, _) in self._vars.values():
+            if output_name in self._conversion_factors:
+                scale, offset = self._conversion_factors[output_name]
+                outputs[output_name] = scale * (inputs[input_name] + offset)
+            else:
+                outputs[output_name] = inputs[input_name]
diff --git a/dymos/transcriptions/explicit_shooting/rk_integration_comp.py b/dymos/transcriptions/explicit_shooting/rk_integration_comp.py
index 099c5e6e..9ef7a749 100644
--- a/dymos/transcriptions/explicit_shooting/rk_integration_comp.py
+++ b/dymos/transcriptions/explicit_shooting/rk_integration_comp.py
@@ -577,22 +577,22 @@ class RKIntegrationComp(om.ExplicitComponent):
         self._timeseries_idxs_in_y = {}
         self._filtered_timeseries_outputs = {}
 
-        for ts_name, ts_opts in self.timeseries_options.items():
-            patterns = [output['name'] for output in ts_opts['outputs']]
+        for ts_opts in self.timeseries_options.values():
+            patterns = [output['name'] for output in ts_opts['outputs'].values()]
             matching_outputs = filter_outputs(patterns, ode_outputs)
 
             explicit_requests = set([var for var in patterns if '*' not in var])
 
-            unmatched_requests = sorted(list(set(explicit_requests) - set(matching_outputs.keys())))
+            unmatched_requests = explicit_requests.difference(matching_outputs)
 
             if unmatched_requests:
                 om.issue_warning(msg='The following timeseries outputs were requested but '
-                                     f'not found in the ODE: {", ".join(unmatched_requests)}',
+                                     f'not found in the ODE: {", ".join(sorted(unmatched_requests))}',
                                  category=om.OpenMDAOWarning)
 
             for var, var_meta in matching_outputs.items():
                 if var in explicit_requests:
-                    ts_output = next((output for output in ts_opts['outputs'] if output['name'] == var))
+                    ts_output = next((output for output in ts_opts['outputs'].values() if output['name'] == var))
                     # var explicitly matched
                     output_name = ts_output['output_name'] if ts_output['output_name'] else ts_output['name'].split('.')[-1]
                     units = ts_output['units'] or var_meta.get('units', None)
@@ -1545,12 +1545,12 @@ class RKIntegrationComp(om.ExplicitComponent):
         outputs['time_phase'] = self._t[idxs, ...] - inputs['t_initial']
 
         # Extract the state values
-        for state_name, options in self.state_options.items():
+        for state_name in self.state_options:
             of = self._state_output_names[state_name]
             outputs[of] = self._x[idxs, self.state_idxs[state_name]]
 
         # Extract the control values and rates
-        for control_name, options in self.control_options.items():
+        for control_name in self.control_options:
             oname = self._control_output_names[control_name]
             rate_name = self._control_rate_names[control_name]
             rate2_name = self._control_rate2_names[control_name]
@@ -1559,7 +1559,7 @@ class RKIntegrationComp(om.ExplicitComponent):
             outputs[rate2_name] = self._y[idxs, self._control_rate2_idxs_in_y[control_name]]
 
         # Extract the control values and rates
-        for control_name, options in self.polynomial_control_options.items():
+        for control_name in self.polynomial_control_options:
             oname = self._polynomial_control_output_names[control_name]
             rate_name = self._polynomial_control_rate_names[control_name]
             rate2_name = self._polynomial_control_rate2_names[control_name]
@@ -1568,7 +1568,7 @@ class RKIntegrationComp(om.ExplicitComponent):
             outputs[rate2_name] = self._y[idxs, self._polynomial_control_rate2_idxs_in_y[control_name]]
 
         # Extract the timeseries outputs
-        for name, options in self._filtered_timeseries_outputs.items():
+        for name in self._filtered_timeseries_outputs:
             oname = self._timeseries_output_names[name]
             outputs[oname] = self._y[idxs, self._timeseries_idxs_in_y[name]]
 
diff --git a/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py b/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py
index 77b40702..782efcde 100644
--- a/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py
+++ b/dymos/transcriptions/pseudospectral/components/pseudospectral_timeseries_output_comp.py
@@ -40,6 +40,14 @@ class PseudospectralTimeseriesOutputComp(TimeseriesOutputCompBase):
         # Flag to set if no multiplication by the interpolation matrix is necessary
         self._no_interp = False
 
+    def _declare_options(self):
+        """
+        Declare options before kwargs are processed in the init method.
+        """
+        super()._declare_options()
+        self.options.declare('time_units', default=None, allow_none=True, types=str,
+                             desc='Units of time')
+
     def setup(self):
         """
         Define the independent variables as output variables.
@@ -62,7 +70,8 @@ class PseudospectralTimeseriesOutputComp(TimeseriesOutputCompBase):
         # To do this, find the nodes in the output grid which fall in each segment of the input
         # grid.  Then build a Lagrange interpolating polynomial for that segment
         L_blocks = []
-        output_nodes_ptau = ogd.node_ptau[ogd.subset_node_indices[output_subset]].tolist()
+        D_blocks = []
+        output_nodes_ptau = ogd.node_ptau[ogd.subset_node_indices[output_subset]]
 
         for iseg in range(igd.num_segments):
             i1, i2 = igd.segment_indices[iseg]
@@ -75,32 +84,31 @@ class PseudospectralTimeseriesOutputComp(TimeseriesOutputCompBase):
             else:
                 ptau_hi = igd.segment_ends[iseg+1]
                 if iseg < igd.num_segments - 1:
-                    idxs_in_iseg = np.where(output_nodes_ptau <= ptau_hi)[0]
+                    optau_segi = output_nodes_ptau[output_nodes_ptau <= ptau_hi]
                 else:
-                    idxs_in_iseg = np.arange(len(output_nodes_ptau))
-                optau_segi = np.asarray(output_nodes_ptau)[idxs_in_iseg]
+                    optau_segi = output_nodes_ptau
+
                 # Remove the captured nodes so we don't accidentally include them again
-                output_nodes_ptau = output_nodes_ptau[len(idxs_in_iseg):]
+                output_nodes_ptau = output_nodes_ptau[len(optau_segi):]
 
             # # Now get the output nodes which fall in iseg in iseg's segment tau space.
             ostau_segi = 2.0 * (optau_segi - iptau_segi[0]) / (iptau_segi[-1] - iptau_segi[0]) - 1
 
             # Create the interpolation matrix and add it to the blocks
-            L, _ = lagrange_matrices(istau_segi, ostau_segi)
+            L, D = lagrange_matrices(istau_segi, ostau_segi)
             L_blocks.append(L)
+            D_blocks.append(D)
 
         if _USE_SPARSE:
             self.interpolation_matrix = sp.block_diag(L_blocks, format='csr')
+            self.differentiation_matrix = sp.block_diag(D_blocks, format='csr')
         else:
             self.interpolation_matrix = block_diag(*L_blocks)
+            self.differentiation_matrix = block_diag(*D_blocks)
 
-        for (name, kwargs) in self._timeseries_outputs:
-            units = kwargs['units']
-            desc = kwargs['units']
-            shape = kwargs['shape']
-            self._add_output_configure(name, units, shape, desc)
+        self.add_input('dt_dstau', shape=(self.input_num_nodes,), units=self.options['time_units'])
 
-    def _add_output_configure(self, name, units, shape, desc='', src=None):
+    def _add_output_configure(self, name, units, shape, desc='', src=None, rate=False):
         """
         Add a single timeseries output.
 
@@ -120,6 +128,8 @@ class PseudospectralTimeseriesOutputComp(TimeseriesOutputCompBase):
             description of the timeseries output variable.
         src : str
             The src path of the variables input, used to prevent redundant inputs.
+        rate : bool
+            If True, timeseries output is a rate.
 
         Returns
         -------
@@ -148,40 +158,40 @@ class PseudospectralTimeseriesOutputComp(TimeseriesOutputCompBase):
             added_source = True
 
         output_name = name
-        self.add_output(output_name,
-                        shape=(output_num_nodes,) + shape,
-                        units=units, desc=desc)
+        self.add_output(output_name, shape=(output_num_nodes,) + shape, units=units, desc=desc)
 
-        self._vars[name] = (input_name, output_name, shape)
+        self._vars[name] = (input_name, output_name, shape, rate)
 
         size = np.prod(shape)
-        val_jac = np.zeros((output_num_nodes, size, input_num_nodes, size))
+        jac = np.zeros((output_num_nodes, size, input_num_nodes, size))
+
+        if rate:
+            mat = self.differentiation_matrix
+        else:
+            mat = self.interpolation_matrix
 
         for i in range(size):
             if _USE_SPARSE:
-                val_jac[:, i, :, i] = self.interpolation_matrix.toarray()
+                jac[:, i, :, i] = mat.toarray()
             else:
-                val_jac[:, i, :, i] = self.interpolation_matrix
-
-        val_jac = val_jac.reshape((output_num_nodes * size, input_num_nodes * size),
-                                  order='C')
+                jac[:, i, :, i] = mat
 
-        val_jac_rows, val_jac_cols = np.where(val_jac != 0)
+        jac = jac.reshape((output_num_nodes * size, input_num_nodes * size), order='C')
 
-        rs, cs = val_jac_rows, val_jac_cols
+        jac_rows, jac_cols = np.nonzero(jac)
 
         # There's a chance that the input for this output was pulled from another variable with
         # different units, so account for that with a conversion.
-        if None in {input_units, units}:
-            scale = 1.0
-            offset = 0
+        if input_units is None or units is None:
+            self.declare_partials(of=output_name, wrt=input_name,
+                                  rows=jac_rows, cols=jac_cols, val=jac[jac_rows, jac_cols])
         else:
             scale, offset = unit_conversion(input_units, units)
-        self._conversion_factors[output_name] = scale, offset
+            self._conversion_factors[output_name] = scale, offset
 
-        self.declare_partials(of=output_name,
-                              wrt=input_name,
-                              rows=rs, cols=cs, val=scale * val_jac[rs, cs])
+            self.declare_partials(of=output_name, wrt=input_name,
+                                  rows=jac_rows, cols=jac_cols,
+                                  val=scale * jac[jac_rows, jac_cols])
 
         return added_source
 
@@ -196,9 +206,10 @@ class PseudospectralTimeseriesOutputComp(TimeseriesOutputCompBase):
         outputs : `Vector`
             `Vector` containing outputs.
         """
-        for (input_name, output_name, _) in self._vars.values():
-            scale, offset = self._conversion_factors[output_name]
+        # convert dt_dstau to a column vector
+        dt_dstau = inputs['dt_dstau'][:, np.newaxis]
 
+        for (input_name, output_name, _, is_rate) in self._vars.values():
             if self._no_interp:
                 interp_vals = inputs[input_name]
             else:
@@ -208,4 +219,12 @@ class PseudospectralTimeseriesOutputComp(TimeseriesOutputCompBase):
                     inp = inp.swapaxes(0, 1)
 
                 interp_vals = self.interpolation_matrix.dot(inp)
-            outputs[output_name] = scale * (interp_vals + offset)
+
+            if is_rate:
+                interp_vals = self.differentiation_matrix.dot(interp_vals) / dt_dstau
+
+            if output_name in self._conversion_factors:
+                scale, offset = self._conversion_factors[output_name]
+                outputs[output_name] = scale * (interp_vals + offset)
+            else:
+                outputs[output_name] = interp_vals
diff --git a/dymos/transcriptions/pseudospectral/gauss_lobatto.py b/dymos/transcriptions/pseudospectral/gauss_lobatto.py
index 58aa9dfa..bc250ab9 100644
--- a/dymos/transcriptions/pseudospectral/gauss_lobatto.py
+++ b/dymos/transcriptions/pseudospectral/gauss_lobatto.py
@@ -3,7 +3,7 @@ import warnings
 
 import numpy as np
 import openmdao.api as om
-from openmdao.utils.general_utils import simple_warning
+from openmdao.utils.om_warnings import issue_warning
 
 from .pseudospectral_base import PseudospectralBase
 from .components import GaussLobattoInterleaveComp
@@ -518,12 +518,12 @@ class GaussLobatto(PseudospectralBase):
                     # phase.promotes(timeseries_name, inputs=[(tgt_name, prom_name)],
                     #                src_indices=(src_idxs,), flat_src_indices=True)
 
-            for ts_output in phase._timeseries[timeseries_name]['outputs']:
+            for output_name, ts_output in phase._timeseries[timeseries_name]['outputs'].items():
                 var = ts_output['name']
-                output_name = ts_output['output_name']
                 units = ts_output['units']
                 wildcard_units = ts_output['wildcard_units']
                 shape = ts_output['shape']
+                is_rate = ts_output['is_rate']
 
                 if '*' in var:  # match outputs from the ODE
                     matches = filter(list(ode_outputs.keys()), var)
@@ -543,10 +543,10 @@ class GaussLobatto(PseudospectralBase):
                     for output_name, var_list in output_name_groups.items():
                         if len(var_list) > 1:
                             var_list_as_string = ', '.join(var_list)
-                            simple_warning(f"The timeseries variable name {output_name} is "
-                                           f"duplicated in these variables: {var_list_as_string}. "
-                                           "Disambiguate by using the add_timeseries_output "
-                                           "output_name option.")
+                            issue_warning(f"The timeseries variable name {output_name} is "
+                                          f"duplicated in these variables: {var_list_as_string}. "
+                                          "Disambiguate by using the add_timeseries_output "
+                                          "output_name option.")
                 else:
                     matches = [var]
 
@@ -585,7 +585,8 @@ class GaussLobatto(PseudospectralBase):
 
                     timeseries_input_added = timeseries_comp._add_output_configure(output_name, units, shape,
                                                                                    src=f'interleave_comp.'
-                                                                                       f'all_values:{output_name}')
+                                                                                       f'all_values:{output_name}',
+                                                                                   rate=is_rate)
 
                     interleave_comp = phase._get_subsystem('interleave_comp')
                     src_added = interleave_comp.add_var(output_name, shape, units,
diff --git a/dymos/transcriptions/pseudospectral/pseudospectral_base.py b/dymos/transcriptions/pseudospectral/pseudospectral_base.py
index 32a48711..fad215da 100644
--- a/dymos/transcriptions/pseudospectral/pseudospectral_base.py
+++ b/dymos/transcriptions/pseudospectral/pseudospectral_base.py
@@ -45,13 +45,13 @@ class PseudospectralBase(TranscriptionBase):
         phase : dymos.Phase
             The phase object to which this transcription instance applies.
         """
-        time_units = phase.time_options['units']
         grid_data = self.grid_data
 
         super(PseudospectralBase, self).setup_time(phase)
 
         time_comp = TimeComp(num_nodes=grid_data.num_nodes, node_ptau=grid_data.node_ptau,
-                             node_dptau_dstau=grid_data.node_dptau_dstau, units=time_units,
+                             node_dptau_dstau=grid_data.node_dptau_dstau,
+                             units=phase.time_options['units'],
                              initial_val=phase.time_options['initial_val'],
                              duration_val=phase.time_options['duration_val'])
 
@@ -252,14 +252,11 @@ class PseudospectralBase(TranscriptionBase):
             The phase object to which this transcription instance applies.
         """
         grid_data = self.grid_data
-        transcription = grid_data.transcription
-        time_units = phase.time_options['units']
-
         phase.add_subsystem('state_interp',
                             subsys=StateInterpComp(grid_data=grid_data,
                                                    state_options=phase.state_options,
-                                                   time_units=time_units,
-                                                   transcription=transcription))
+                                                   time_units=phase.time_options['units'],
+                                                   transcription=grid_data.transcription))
 
     def configure_ode(self, phase):
         """
@@ -292,14 +289,10 @@ class PseudospectralBase(TranscriptionBase):
         phase : dymos.Phase
             The phase object to which this transcription instance applies.
         """
-        grid_data = self.grid_data
-
-        time_units = phase.time_options['units']
-
         phase.add_subsystem('collocation_constraint',
-                            CollocationComp(grid_data=grid_data,
+                            CollocationComp(grid_data=self.grid_data,
                                             state_options=phase.state_options,
-                                            time_units=time_units))
+                                            time_units=phase.time_options['units']))
 
     def _configure_solve_segments(self, state_name, options, phase):
         """
@@ -498,9 +491,12 @@ class PseudospectralBase(TranscriptionBase):
 
             timeseries_comp = PseudospectralTimeseriesOutputComp(input_grid_data=gd,
                                                                  output_grid_data=ogd,
-                                                                 output_subset=options['subset'])
+                                                                 output_subset=options['subset'],
+                                                                 time_units=phase.time_options['units'])
             phase.add_subsystem(name, subsys=timeseries_comp)
 
+            phase.connect('dt_dstau', (f'{name}.dt_dstau'), flat_src_indices=True)
+
     def _get_objective_src(self, var, loc, phase, ode_outputs=None):
         """
         Return the path to the variable that will be used as the objective.
@@ -557,7 +553,7 @@ class PseudospectralBase(TranscriptionBase):
             else:
                 linear = False
             constraint_path = f'states:{var}'
-        elif var_type in 'indep_control':
+        elif var_type == 'indep_control':
             shape = phase.control_options[var]['shape']
             units = phase.control_options[var]['units']
             linear = True
@@ -567,7 +563,7 @@ class PseudospectralBase(TranscriptionBase):
             units = phase.control_options[var]['units']
             linear = False
             constraint_path = f'control_values:{var}'
-        elif var_type in 'indep_polynomial_control':
+        elif var_type == 'indep_polynomial_control':
             shape = phase.polynomial_control_options[var]['shape']
             units = phase.polynomial_control_options[var]['units']
             linear = True
diff --git a/dymos/transcriptions/pseudospectral/radau_pseudospectral.py b/dymos/transcriptions/pseudospectral/radau_pseudospectral.py
index baca78a7..4f601280 100644
--- a/dymos/transcriptions/pseudospectral/radau_pseudospectral.py
+++ b/dymos/transcriptions/pseudospectral/radau_pseudospectral.py
@@ -4,7 +4,7 @@ import warnings
 
 import numpy as np
 import openmdao.api as om
-from openmdao.utils.general_utils import simple_warning
+from openmdao.utils.om_warnings import issue_warning
 
 from .pseudospectral_base import PseudospectralBase
 from ..common import RadauPSContinuityComp
@@ -388,12 +388,12 @@ class Radau(PseudospectralBase):
                                       src_indices=(src_idxs,),
                                       flat_src_indices=True)
 
-            for ts_output in phase._timeseries[timeseries_name]['outputs']:
+            for output_name, ts_output in phase._timeseries[timeseries_name]['outputs'].items():
                 var = ts_output['name']
-                output_name = ts_output['output_name']
                 units = ts_output['units']
                 wildcard_units = ts_output['wildcard_units']
                 shape = ts_output['shape']
+                is_rate = ts_output['is_rate']
 
                 if '*' in var:  # match outputs from the ODE
                     matches = filter(list(ode_outputs.keys()), var)
@@ -413,10 +413,10 @@ class Radau(PseudospectralBase):
                     for output_name, var_list in output_name_groups.items():
                         if len(var_list) > 1:
                             var_list_as_string = ', '.join(var_list)
-                            simple_warning(f"The timeseries variable name {output_name} is "
-                                           f"duplicated in these variables: {var_list_as_string}. "
-                                           "Disambiguate by using the add_timeseries_output "
-                                           "output_name option.")
+                            issue_warning(f"The timeseries variable name {output_name} is "
+                                          f"duplicated in these variables: {var_list_as_string}. "
+                                          "Disambiguate by using the add_timeseries_output "
+                                          "output_name option.")
                 else:
                     matches = [var]
 
@@ -453,7 +453,8 @@ class Radau(PseudospectralBase):
 
                     add_connection = timeseries_comp._add_output_configure(output_name, units,
                                                                            shape, desc='',
-                                                                           src=f'rhs_all.{v}')
+                                                                           src=f'rhs_all.{v}',
+                                                                           rate=is_rate)
 
                     if add_connection:
                         phase.connect(src_name=f'rhs_all.{v}',
diff --git a/dymos/transcriptions/solve_ivp/components/solve_ivp_timeseries_comp.py b/dymos/transcriptions/solve_ivp/components/solve_ivp_timeseries_comp.py
index 62afbc6a..3c5c9987 100644
--- a/dymos/transcriptions/solve_ivp/components/solve_ivp_timeseries_comp.py
+++ b/dymos/transcriptions/solve_ivp/components/solve_ivp_timeseries_comp.py
@@ -31,13 +31,7 @@ class SolveIVPTimeseriesOutputComp(TimeseriesOutputCompBase):
         else:
             self.num_nodes = grid_data.num_segments * self.options['output_nodes_per_seg']
 
-        for (name, kwargs) in self._timeseries_outputs:
-            units = kwargs['units']
-            desc = kwargs['units']
-            shape = kwargs['shape']
-            self._add_output_configure(name, units, shape, desc)
-
-    def _add_output_configure(self, name, units, shape, desc):
+    def _add_output_configure(self, name, units, shape, desc, rate=False):
         """
         Add a single timeseries output.
 
@@ -55,20 +49,20 @@ class SolveIVPTimeseriesOutputComp(TimeseriesOutputCompBase):
             Default is None, which means it has no units.
         desc : str
             description of the timeseries output variable.
+        rate : bool
+            If True, timeseries output is a rate.
         """
-        num_nodes = self.num_nodes
+        if rate:
+            raise NotImplementedError("Timeseries output rates are not currently supported for "
+                                      "SolveIVP transcriptions.")
 
+        nodeshape = (self.num_nodes,)
         input_name = f'all_values:{name}'
-        self.add_input(input_name,
-                       shape=(num_nodes,) + shape,
-                       units=units, desc=desc)
 
-        output_name = name
-        self.add_output(output_name,
-                        shape=(num_nodes,) + shape,
-                        units=units, desc=desc)
+        self.add_input(input_name, shape=nodeshape + shape,  units=units, desc=desc)
+        self.add_output(name, shape=nodeshape + shape, units=units, desc=desc)
 
-        self._vars[name] = (input_name, output_name, shape)
+        self._vars[name] = (input_name, name, shape, rate)
 
     def compute(self, inputs, outputs):
         """
@@ -81,5 +75,4 @@ class SolveIVPTimeseriesOutputComp(TimeseriesOutputCompBase):
         outputs : `Vector`
             `Vector` containing outputs.
         """
-        for (input_name, output_name, _) in self._vars.values():
-            outputs[output_name] = inputs[input_name]
+        outputs.set_val(inputs.asarray())
diff --git a/dymos/transcriptions/solve_ivp/solve_ivp.py b/dymos/transcriptions/solve_ivp/solve_ivp.py
index a890cfc4..be10e8f6 100644
--- a/dymos/transcriptions/solve_ivp/solve_ivp.py
+++ b/dymos/transcriptions/solve_ivp/solve_ivp.py
@@ -633,22 +633,18 @@ class SolveIVP(TranscriptionBase):
                     src_idxs_raw = np.zeros(num_seg * output_nodes_per_seg, dtype=int)
                 src_idxs = get_src_indices_by_row(src_idxs_raw, options['shape']).ravel()
 
-                # tgt_name = f'all_values:parameters:{name}'
-                # phase.promotes('timeseries', inputs=[(tgt_name, prom_name)],
-                #                src_indices=(src_idxs,), flat_src_indices=True)
-                # print(src_idxs)
                 phase.connect(f'parameter_vals:{name}', f'timeseries.all_values:parameters:{name}',
                               src_indices=src_idxs, flat_src_indices=True)
 
-        for ts_output in phase._timeseries['timeseries']['outputs']:
+        for output_name, ts_output in phase._timeseries['timeseries']['outputs'].items():
             var = ts_output['name']
-            output_name = ts_output['output_name']
             units = ts_output['units']
             wildcard_units = ts_output['wildcard_units']
             shape = ts_output['shape']
 
             if '*' in var:  # match outputs from the ODE
-                matches = filter(list(ode_outputs.keys()), var)
+                # TODO: is filter still case INSENSTIVE on windows?  If so, fix this
+                matches = filter(ode_outputs.keys(), var)
             else:
                 matches = [var]
 


Test patch:
diff --git a/dymos/examples/brachistochrone/test/test_duplicate_constraints.py b/dymos/examples/brachistochrone/test/test_duplicate_constraints.py
index 4e18cba9..4f8e6fcb 100644
--- a/dymos/examples/brachistochrone/test/test_duplicate_constraints.py
+++ b/dymos/examples/brachistochrone/test/test_duplicate_constraints.py
@@ -614,4 +614,4 @@ class TestDuplicateConstraints(unittest.TestCase):
             warnings.simplefilter("error")
             phase.add_path_constraint('check', constraint_name='bounded_check', upper=100, lower=-100)
 
-        self.assertIn('bounded_check', [op['output_name'] for op in phase._timeseries['timeseries']['outputs']])
+        self.assertIn('bounded_check', [op['output_name'] for op in phase._timeseries['timeseries']['outputs'].values()])
diff --git a/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_bokeh_plots.py b/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_bokeh_plots.py
index 6f1f4ca3..3e9c6b0f 100644
--- a/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_bokeh_plots.py
+++ b/dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_bokeh_plots.py
@@ -216,7 +216,7 @@ class TestExampleTwoBurnOrbitRaise(unittest.TestCase):
                           'states_accel.png', 'controls_u1.png', 'states_vr.png', 'pos_x.png',
                           'states_vt.png', 'pos_y.png', 'parameters_u1.png', 'states_theta.png',
                           'control_rates_u1_rate2.png', 'state_rates_vt.png', 'time_phase.png',
-                          'parameters_c.png', 'state_rates_theta.png', 'state_rates_vr.png'}
+                          'parameters_c.png', 'state_rates_theta.png', 'state_rates_vr.png', 'dt_dstau.png'}
 
         self.assertSetEqual(expected_files, set(os.listdir('plots')))
 
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 94dfe816..29d2eab3 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
@@ -1,5 +1,6 @@
 import unittest
-
+import numpy as np
+from numpy.polynomial import Polynomial as P
 import openmdao.api as om
 from openmdao.utils.assert_utils import assert_near_equal
 import dymos as dm
@@ -7,9 +8,8 @@ from dymos.examples.min_time_climb.min_time_climb_ode import MinTimeClimbODE
 from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
 
 
-@require_pyoptsparse(optimizer='SLSQP')
 def min_time_climb(optimizer='SLSQP', num_seg=3, transcription='gauss-lobatto',
-                   transcription_order=3, force_alloc_complex=False):
+                   transcription_order=3, force_alloc_complex=False, add_rate=False):
 
     p = om.Problem(model=om.Group())
 
@@ -25,6 +25,12 @@ def min_time_climb(optimizer='SLSQP', num_seg=3, transcription='gauss-lobatto',
         p.driver.opt_settings['Function precision'] = 1.0E-12
         p.driver.opt_settings['Linesearch tolerance'] = 0.1
         p.driver.opt_settings['Major step limit'] = 0.5
+    elif optimizer == 'IPOPT':
+        p.driver.opt_settings['tol'] = 1.0E-5
+        p.driver.opt_settings['print_level'] = 5
+        p.driver.opt_settings['mu_strategy'] = 'monotone'
+        p.driver.opt_settings['bound_mult_init_method'] = 'mu-based'
+        p.driver.opt_settings['mu_init'] = 0.01
 
     t = {'gauss-lobatto': dm.GaussLobatto(num_segments=num_seg, order=transcription_order),
          'radau-ps': dm.Radau(num_segments=num_seg, order=transcription_order)}
@@ -86,6 +92,10 @@ def min_time_climb(optimizer='SLSQP', num_seg=3, transcription='gauss-lobatto',
     phase.add_timeseries_output(['aero.*', 'prop.thrust', 'prop.m_dot'],
                                 units={'aero.f_lift': 'lbf', 'prop.thrust': 'lbf'})
 
+    # test adding rate as timeseries output
+    if add_rate:
+        phase.add_timeseries_rate_output('aero.mach')
+
     p.model.linear_solver = om.DirectSolver()
 
     p.setup(check=True, force_alloc_complex=force_alloc_complex)
@@ -108,48 +118,111 @@ def min_time_climb(optimizer='SLSQP', num_seg=3, transcription='gauss-lobatto',
 @use_tempdirs
 class TestMinTimeClimb(unittest.TestCase):
 
-    def test_results_gauss_lobatto(self):
-        p = min_time_climb(optimizer='SLSQP', num_seg=12, transcription_order=3,
-                           transcription='gauss-lobatto')
+    def _test_results(self, p):
+        """ Verify the results of the optimization. """
+        # Verify that ODE output mach is added to the timeseries
+        assert_near_equal(p.get_val('traj.phase0.timeseries.mach')[-1], 1.0, tolerance=1.0E-2)
 
         # Check that time matches to within 1% of an externally verified solution.
         assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 321.0, tolerance=0.02)
 
-        # Verify that ODE output mach is added to the timeseries
-        assert_near_equal(p.get_val('traj.phase0.timeseries.mach')[-1], 1.0, tolerance=1.0E-2)
-
-        # verify all wildcard timeseries exist
-        output_dict = dict(p.model.list_outputs(units=True))
-        ts = [k for k, v in output_dict.items() if 'timeseries' in k]
+    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))
+        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]))
-        # verify time series units
+
+    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
 
+    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))
+        ts = {k: v for k, v in output_dict.items() if 'timeseries.' in k}
+        self.assertTrue('traj.phases.phase0.timeseries.mach_rate' in ts)
+
+        time = p['traj.phases.phase0.timeseries.time'][:, 0]
+        mach = p['traj.phases.phase0.timeseries.mach'][:, 0]
+        mach_rate = p['traj.phases.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
+        seg_idxs = gd.subset_segment_indices['all']
+        num_seg = seg_idxs.shape[0]
+
+        if plot:
+            import matplotlib.pyplot as plt
+            import matplotlib.cm as cm
+            fig, axes = plt.subplots(2, 1, sharex=True)
+            color = iter(cm.viridis(np.linspace(0, 1, num_seg)))
+            axes[0].set_ylabel('Mach')
+            axes[1].set_ylabel('Mach rate (1/s)')
+            axes[1].set_xlabel('time (s)')
+
+        for i in range(num_seg):
+            time_seg = time[seg_idxs[i, 0]: seg_idxs[i, 1]]
+            mach_seg = mach[seg_idxs[i, 0]: seg_idxs[i, 1]]
+            order = len(time_seg) - 1
+            p = P.fit(time_seg, mach_seg, order)
+            deriv = p.deriv(1)
+
+            time_nodes = time[seg_idxs[i, 0]: seg_idxs[i, 1]]
+            mach_nodes = mach[seg_idxs[i, 0]: seg_idxs[i, 1]]
+            mach_rate_nodes = mach_rate[seg_idxs[i, 0]: seg_idxs[i, 1]]
+
+            if plot:
+                c = next(color)
+                t_plot = np.linspace(time_seg[0], time_seg[-1], 20)
+                m_plot = p(t_plot)
+                m_rate_plot = deriv(t_plot)
+                axes[0].plot(t_plot, m_plot, c=c)
+                axes[1].plot(t_plot, m_rate_plot, '--', c=c)
+
+                axes[0].plot(time[seg_idxs[i, 0]: seg_idxs[i, 1]], mach[seg_idxs[i, 0]: seg_idxs[i, 1]], 'o', c=c)
+                axes[1].plot(time[seg_idxs[i, 0]: seg_idxs[i, 1]], mach_rate[seg_idxs[i, 0]: seg_idxs[i, 1]], 'o', c=c)
+
+            assert_near_equal(mach_nodes, p(time_nodes), tolerance=1.0E-9)
+            assert_near_equal(mach_rate_nodes, deriv(time_nodes), tolerance=1.0E-9)
+
+        if plot:
+            plt.show()
+
+    @require_pyoptsparse(optimizer='SLSQP')
+    def test_results_gauss_lobatto(self):
+        NUM_SEG = 12
+        ORDER = 3
+        p = min_time_climb(optimizer='IPOPT', num_seg=NUM_SEG, transcription_order=ORDER,
+                           transcription='gauss-lobatto', add_rate=True)
+
+        self._test_results(p)
+
+        self._test_wilcard_outputs(p)
+
+        self._test_timeseries_units(p)
+
+        self._test_mach_rate(p)
+
+    @require_pyoptsparse(optimizer='SLSQP')
     def test_results_radau(self):
-        p = min_time_climb(optimizer='SLSQP', num_seg=12, transcription_order=3,
-                           transcription='radau-ps')
+        NUM_SEG = 12
+        ORDER = 3
+        p = min_time_climb(optimizer='SLSQP', num_seg=NUM_SEG, transcription_order=ORDER,
+                           transcription='radau-ps', add_rate=True)
 
-        # Check that time matches to within 1% of an externally verified solution.
-        assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 321.0, tolerance=0.02)
+        self._test_results(p)
 
-        # Verify that ODE output mach is added to the timeseries
-        assert_near_equal(p.get_val('traj.phase0.timeseries.mach')[-1], 1.0, tolerance=1.0E-2)
+        self._test_wilcard_outputs(p)
 
-        # verify all wildcard timeseries exist
-        output_dict = dict(p.model.list_outputs(units=True))
-        ts = [k 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]))
-        # verify time series units
-        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
+        self._test_timeseries_units(p)
 
+        self._test_mach_rate(p)
 
 if __name__ == '__main__':  # pragma: no cover
     unittest.main()