Deprecated SolveIVP transcription. simulate method now uses the ExplicitShooting transcription without derivatives.

.github/workflows/dymos_docs_workflow.yml

@@ -72,7 +72,7 @@ jobs:
           echo "$SSH_KNOWN_HOSTS" > ~/.ssh/known_hosts
 
       - name: Checkout code
-        uses: actions/checkout@v2
+        uses: actions/checkout@v3
 
       - name: Fetch tags
         run: |

.github/workflows/dymos_tests_workflow.yml

@@ -125,7 +125,7 @@ jobs:
 
       - name: Checkout code
         if: env.RUN_BUILD
-        uses: actions/checkout@v2
+        uses: actions/checkout@v3
 
       - name: Fetch tags
         if: env.RUN_BUILD

docs/dymos_book/examples/robertson_problem/robertson_problem.ipynb

@@ -240,10 +240,19 @@
    ]
   },
   {
+   "attachments": {},
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "## Building and running the problem"
+    "## Building and running the problem\n",
+    "\n",
+    "Here we're using the ExplicitShooting transcription in Dymos.\n",
+    "The ExplicitShooting transcription explicit integrates the given ODE using the `solve_ivp` method of [scipy](https://scipy.org/).\n",
+    "\n",
+    "Since this is purely an integration with no controls to be determined, a single call to `run_model` will propagate the solution for us. There's no need to run a driver. Even the typical follow-up call to `traj.simulate` is unnecessary.\n",
+    "\n",
+    "Technically, we could even solve this using a single segment since segment spacing in the explicit shooting transcription determines the density of the control nodes, and there are no controls for this simulation.\n",
+    "Providing more segments in this case (or a higher segment order) increases the number of nodes at which the outputs are provided."
    ]
   },
   {
@@ -265,11 +274,11 @@
     "    # Create a trajectory and add a phase to it\n",
     "    #\n",
     "    traj = p.model.add_subsystem('traj', dm.Trajectory())\n",
+    "    \n",
+    "    tx = dm.ExplicitShooting(num_segments=10, method='LSODA')\n",
     "\n",
     "    phase = traj.add_phase('phase0',\n",
-    "                           dm.Phase(ode_class=RobertsonODE,\n",
-    "                                    transcription=dm.GaussLobatto(num_segments=50)\n",
-    "                                   ))\n",
+    "                           dm.Phase(ode_class=RobertsonODE, transcription=tx))\n",
     "\n",
     "    #\n",
     "    # Set the variables\n",
@@ -311,9 +320,7 @@
     "# just set up the problem, test it elsewhere\n",
     "p = robertson_problem(t_final=40)\n",
     "\n",
-    "p.run_model()\n",
-    "\n",
-    "p_sim = p.model.traj.simulate(method='LSODA')"
+    "p.run_model()"
    ]
   },
   {
@@ -329,7 +336,7 @@
    "source": [
     "from openmdao.utils.assert_utils import assert_near_equal\n",
     "\n",
-    "assert_near_equal(p_sim.get_val('traj.phase0.timeseries.states:x0')[-1], 0.71583161, tolerance=1E-5)"
+    "assert_near_equal(p.get_val('traj.phase0.timeseries.states:x0')[-1], 0.71583161, tolerance=1E-4)"
    ]
   },
   {
@@ -345,7 +352,7 @@
    "source": [
     "from openmdao.utils.assert_utils import assert_near_equal\n",
     "\n",
-    "assert_near_equal(p_sim.get_val('traj.phase0.timeseries.states:y0')[-1], 9.18571144e-06, tolerance=1E-1)"
+    "assert_near_equal(p.get_val('traj.phase0.timeseries.states:y0')[-1], 9.18571144e-06, tolerance=1E-4)"
    ]
   },
   {
@@ -361,7 +368,7 @@
    "source": [
     "from openmdao.utils.assert_utils import assert_near_equal\n",
     "\n",
-    "assert_near_equal(p_sim.get_val('traj.phase0.timeseries.states:z0')[-1], 0.2841592, tolerance=1E-5)"
+    "assert_near_equal(p.get_val('traj.phase0.timeseries.states:z0')[-1], 0.2841592, tolerance=1E-4)"
    ]
   },
   {
@@ -372,12 +379,12 @@
    "source": [
     "import matplotlib.pyplot as plt\n",
     "\n",
-    "t_sim = p_sim.get_val('traj.phase0.timeseries.time')\n",
+    "t = p.get_val('traj.phase0.timeseries.time')\n",
     "\n",
     "states = ['x0', 'y0', 'z0']\n",
     "fig, axes = plt.subplots(len(states), 1)\n",
     "for i, state in enumerate(states):\n",
-    "    axes[i].plot(t_sim, p_sim.get_val(f'traj.phase0.timeseries.states:{state}'), '-')\n",
+    "    axes[i].plot(t, p.get_val(f'traj.phase0.timeseries.states:{state}'), 'o')\n",
     "    axes[i].set_ylabel(state[0])\n",
     "axes[-1].set_xlabel('time (s)')\n",
     "plt.tight_layout()\n",
@@ -396,7 +403,7 @@
    }
   },
   "kernelspec": {
-   "display_name": "Python 3",
+   "display_name": "Python 3 (ipykernel)",
    "language": "python",
    "name": "python3"
   },
@@ -410,7 +417,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.7.3"
+   "version": "3.11.0"
   }
  },
  "nbformat": 4,

dymos/examples/aircraft_steady_flight/test/test_aircraft_cruise.py

@@ -83,7 +83,7 @@ def test_cruise_results_gl(self):
         p.setup()
 
         p['phase0.t_initial'] = 0.0
-        p['phase0.t_duration'] = 1.515132 * 3600.0
+        p['phase0.t_duration'] = 3600.0
         p['phase0.states:range'] = phase.interp('range', (0, 1296.4))
         p['phase0.states:mass_fuel'] = phase.interp('mass_fuel', (12236.594555, 1))
         p['phase0.states:alt'] = 5.0
@@ -181,7 +181,7 @@ def test_cruise_results_radau(self):
         p.setup()
 
         p['phase0.t_initial'] = 0.0
-        p['phase0.t_duration'] = 1.515132 * 3600.0
+        p['phase0.t_duration'] = 1.0 * 3600.0
         p['phase0.states:range'] = phase.interp('range', (0, 1296.4))
         p['phase0.states:mass_fuel'] = phase.interp('mass_fuel', (12236.594555, 0))
         p['phase0.states:alt'] = 5.0

dymos/examples/brachistochrone/test/ex_brachistochrone.py

@@ -1,5 +1,4 @@
 import matplotlib
-import matplotlib.pyplot as plt
 
 import openmdao.api as om
 from openmdao.utils.testing_utils import require_pyoptsparse
@@ -30,6 +29,17 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran
         t = dm.Radau(num_segments=num_segments,
                      order=transcription_order,
                      compressed=compressed)
+    elif transcription == 'shooting-gauss-lobatto':
+        grid = dm.GaussLobattoGrid(num_segments=num_segments,
+                                   nodes_per_seg=transcription_order,
+                                   compressed=compressed)
+        t = dm.ExplicitShooting(grid=grid)
+    elif transcription == 'shooting-radau':
+        grid = dm.RadauGrid(num_segments=num_segments,
+                            nodes_per_seg=transcription_order + 1,
+                            compressed=compressed)
+        t = dm.ExplicitShooting(grid=grid)
+
     traj = dm.Trajectory()
     phase = dm.Phase(ode_class=BrachistochroneODE, transcription=t)
     p.model.add_subsystem('traj0', traj)
@@ -81,5 +91,5 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran
 
 if __name__ == '__main__':
     p = brachistochrone_min_time(transcription='gauss-lobatto', num_segments=5, run_driver=True,
-                                 transcription_order=5, compressed=False, optimizer='SNOPT',
+                                 transcription_order=5, compressed=False, optimizer='IPOPT',
                                  solve_segments=False, force_alloc_complex=True)

dymos/examples/brachistochrone/test/test_brachistochrone_subprob_reports.py

@@ -166,8 +166,8 @@ def test_make_sim_reports(self):
             report_subdirs = sorted([e for e in pathlib.Path(get_reports_dir()).iterdir() if e.is_dir()])
 
             # Test that a report subdir was made
-            # # There is the nominal problem, the simulation problem, and a subproblem for each segment in the simulation.
-            self.assertEqual(len(report_subdirs), 12)
+            # There is the nominal problem, the simulation problem, and a subproblem for the simulation.
+            self.assertEqual(len(report_subdirs), 3)
 
             for subdir in report_subdirs:
                 path = pathlib.Path(subdir).joinpath(self.n2_filename)

dymos/examples/brachistochrone/test/test_ex_brachistochrone.py

@@ -55,32 +55,52 @@ def test_ex_brachistochrone_radau_compressed(self):
                                                         compressed=True)
         self.run_asserts(p)
         self.tearDown()
-        if os.path.exists('ex_brach_radau_compressed.db'):
-            os.remove('ex_brach_radau_compressed.db')
 
     def test_ex_brachistochrone_radau_uncompressed(self):
         ex_brachistochrone.SHOW_PLOTS = True
         p = ex_brachistochrone.brachistochrone_min_time(transcription='radau-ps',
                                                         compressed=False)
         self.run_asserts(p)
         self.tearDown()
-        if os.path.exists('ex_brach_radau_uncompressed.db'):
-            os.remove('ex_brach_radau_uncompressed.db')
 
     def test_ex_brachistochrone_gl_compressed(self):
         ex_brachistochrone.SHOW_PLOTS = True
         p = ex_brachistochrone.brachistochrone_min_time(transcription='gauss-lobatto',
                                                         compressed=True)
         self.run_asserts(p)
         self.tearDown()
-        if os.path.exists('ex_brach_gl_compressed.db'):
-            os.remove('ex_brach_gl_compressed.db')
 
     def test_ex_brachistochrone_gl_uncompressed(self):
         ex_brachistochrone.SHOW_PLOTS = True
         p = ex_brachistochrone.brachistochrone_min_time(transcription='gauss-lobatto',
                                                         compressed=False)
         self.run_asserts(p)
         self.tearDown()
-        if os.path.exists('ex_brach_gl_uncompressed.db'):
-            os.remove('ex_brach_gl_uncompressed.db')
+
+    def test_ex_brachistochrone_shooting_gl_compressed(self):
+        ex_brachistochrone.SHOW_PLOTS = True
+        p = ex_brachistochrone.brachistochrone_min_time(transcription='shooting-gauss-lobatto',
+                                                        compressed=True)
+        self.run_asserts(p)
+        self.tearDown()
+
+    def test_ex_brachistochrone_shooting_gl_uncompressed(self):
+        ex_brachistochrone.SHOW_PLOTS = True
+        p = ex_brachistochrone.brachistochrone_min_time(transcription='shooting-gauss-lobatto',
+                                                        compressed=False)
+        self.run_asserts(p)
+        self.tearDown()
+
+    def test_ex_brachistochrone_shooting_radau_compressed(self):
+        ex_brachistochrone.SHOW_PLOTS = True
+        p = ex_brachistochrone.brachistochrone_min_time(transcription='shooting-radau',
+                                                        compressed=True)
+        self.run_asserts(p)
+        self.tearDown()
+
+    def test_ex_brachistochrone_shooting_radau_uncompressed(self):
+        ex_brachistochrone.SHOW_PLOTS = True
+        p = ex_brachistochrone.brachistochrone_min_time(transcription='shooting-radau',
+                                                        compressed=False)
+        self.run_asserts(p)
+        self.tearDown()

dymos/examples/brachistochrone/test/test_state_rate_introspection.py

@@ -239,7 +239,7 @@ def _test_integrate_control_rate2(self, transcription):
         # Define a Dymos Phase object with GaussLobatto Transcription
         #
         phase = dm.Phase(ode_class=BrachistochroneODE,
-                         transcription=transcription(num_segments=10, order=3))
+                         transcription=transcription(num_segments=15, order=3))
 
         traj.add_phase(name='phase0', phase=phase)
 
@@ -304,7 +304,7 @@ def _test_integrate_control_rate2(self, transcription):
         p.set_val('traj.phase0.controls:theta', phase.interp('theta', [0, 100]), units='deg')
 
         # Run the driver to solve the problem
-        dm.run_problem(p, simulate=True, make_plots=False)
+        dm.run_problem(p, simulate=True, make_plots=False, simulate_kwargs={'atol': 1.0E-9, 'rtol': 1.0E-9})
 
         sol = om.CaseReader('dymos_solution.db').get_case('final')
         sim = om.CaseReader('dymos_simulation.db').get_case('final')

dymos/examples/min_time_climb/test/test_min_time_climb_simulate_failure.py

@@ -93,10 +93,10 @@ def test_simulate_raises_analysis_error(self):
 
         with self.assertRaises(om.AnalysisError) as e:
             dm.run_problem(p, run_driver=False, simulate=True)
-        expected_pattern = r"'traj.phases.phase0.segments.segment_[0-9]' <class SegmentSimulationComp>: Error calling " \
-                           r"compute\(\), solve_ivp failed: Required step size is less than spacing " \
-                           r"between numbers. Dynamics changing too dramatically"
-        self.assertRegex(str(e.exception), expected_pattern)
+        expected_pattern = r"'traj.phases.phase0.integrator' <class ODEIntegrationComp>: Error calling " \
+                           r"compute(), solve_ivp failed: Required step size is less than spacing " \
+                           r"between numbers."
+        self.assertEqual(str(e.exception), expected_pattern)
 
 
 if __name__ == '__main__':  # pragma: no cover

dymos/examples/robertson_problem/doc/test_doc_robertson_problem.py

@@ -7,7 +7,7 @@
 from dymos.examples.robertson_problem.doc.robertson_ode import RobertsonODE
 
 
-matplotlib.use('Agg')
+# matplotlib.use('Agg')
 plt.style.use('ggplot')
 
 
@@ -31,9 +31,8 @@ def robertson_problem(self, t_final=1.0):
 
         phase = traj.add_phase('phase0',
                                dm.Phase(ode_class=RobertsonODE,
-                                        transcription=dm.GaussLobatto(num_segments=50)
+                                        transcription=dm.ExplicitShooting(num_segments=10, method='LSODA')
                                         ))
-
         #
         # Set the variables
         #
@@ -93,18 +92,16 @@ def test_robertson_problem_for_docs(self):
 
         p.run_model()
 
-        p_sim = p.model.traj.simulate(method='LSODA')
-
-        assert_near_equal(p_sim.get_val('traj.phase0.timeseries.states:x0')[-1], 0.71583161, tolerance=1E-5)
-        assert_near_equal(p_sim.get_val('traj.phase0.timeseries.states:y0')[-1], 9.18571144e-06, tolerance=1E-1)
-        assert_near_equal(p_sim.get_val('traj.phase0.timeseries.states:z0')[-1], 0.2841592, tolerance=1E-5)
+        assert_near_equal(p.get_val('traj.phase0.timeseries.states:x0')[-1], 0.71583161, tolerance=1E-4)
+        assert_near_equal(p.get_val('traj.phase0.timeseries.states:y0')[-1], 9.18571144e-06, tolerance=1E-4)
+        assert_near_equal(p.get_val('traj.phase0.timeseries.states:z0')[-1], 0.2841592, tolerance=1E-4)
 
-        t_sim = p_sim.get_val('traj.phase0.timeseries.time')
+        t_sim = p.get_val('traj.phase0.timeseries.time')
 
         states = ['x0', 'y0', 'z0']
         fig, axes = plt.subplots(len(states), 1)
         for i, state in enumerate(states):
-            axes[i].plot(t_sim, p_sim.get_val(f'traj.phase0.timeseries.states:{state}'), '-')
+            axes[i].plot(t_sim, p.get_val(f'traj.phase0.timeseries.states:{state}'), 'o')
             axes[i].set_ylabel(state[0])
         axes[-1].set_xlabel('time (s)')
         plt.tight_layout()

dymos/phase/__init__.py

@@ -1,2 +1,3 @@
 from .analytic_phase import AnalyticPhase
 from .phase import Phase
+from .simulation_phase import SimulationPhase

dymos/phase/options.py

@@ -611,7 +611,7 @@ class SimulateOptionsDictionary(om.OptionsDictionary):
     def __init__(self, read_only=False):
         super(SimulateOptionsDictionary, self).__init__(read_only)
 
-        self.declare('method', values=('RK23', 'RK45', 'DOP853', 'BDF', 'Radau', 'LSODA'), default='RK45',
+        self.declare('method', values=('RK23', 'RK45', 'DOP853', 'BDF', 'Radau', 'LSODA'), default='DOP853',
                      desc='The method used by simulate to propagate the ODE.')
 
         self.declare(name='atol', types=(float, np.array), default=1.0E-6,