Added timeseries rates for pseudospectral transcriptions (#781)

* interim

* got rid of unnecessary unit conversions (when scaler==1. and adder=0)

* changed timeseries['outputs'] to a dict keyed on output_name

* started adding differentiation matrix stuff

* getting closer, but single segment order=21 case doesn't converge

* moved some unrelated fixes to another branch

* updated timeseries rate output test to more accurately test rates

* added an add_timeseries_output_rate method

* added time_units

Co-authored-by: Rob Falck <[email protected]>

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",

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",

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",

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."
    ]

dymos/examples/brachistochrone/test/test_duplicate_constraints.py

@@ -614,4 +614,4 @@ def test_duplicate_timeseries_output(self):
             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()])

dymos/examples/finite_burn_orbit_raise/test/test_ex_two_burn_orbit_raise_bokeh_plots.py

@@ -216,7 +216,7 @@ def test_mpl_plots(self):
                           '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')))
 

dymos/examples/min_time_climb/test/test_ex_min_time_climb.py

@@ -1,15 +1,15 @@
 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
 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()

dymos/phase/options.py

@@ -685,3 +685,6 @@ def __init__(self, read_only=False):
 
         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.')