Birkhoff Transcription Take 2 (#1008)

* Added Birkhoff transcription and modified pseudospectral base to accomodate it

* modified defects to match new paper

* Added a birkhoff transcription for GL grid

* Fixed a bug in time conversion

* attempt to compute defects using the matrix form

* typo

* changed collocation comp to use vectorized state

* removed duplicate Birhkoff transcription for GL grid

* pre-computed A and C matrics

* modified birhkoff collocation to be standalone transcription independent of pseudospectral base

* Demonstration of a functional birkhoff iter group under LGL and solve_segments=forward

* better testing

* added backwards propagation test

* updated to newest algorithm from Ross

* Added analytic derivatives

* Updated Birkhoff.get_parameter_connections to accommodate dymos changes.

* Fixed indexing on collocaition constraint test

* Changed the way state rates are calculated in BirkhoffCollocationComp to better match the referenced.  Updated partials.

* more partials cleanup

* Minor fixes to collocation comp

* adding birkhoff radau capability.

* fixed a unit mismatch in collocation comp

* Added CGL collocation

* Modified birkhoff collocation comp to allow for matrix states

* modified birkhoff collocation comp to work with matrix states

* made SimulationPhase and ExplicitShooting transcription work with a birkhoff grid

* Minor fixes to the birkhoff collocation comp test

* Fixed an issue with the cgl weights and added tests for the cgl nodes and birkhoff matrices

* Added birkhoff transcription and grid to init

* birkhoff transcription now uses the birkhoff_iter_group

* Added moon landing problem from Ross's paper

* added an endpoint ODE and started connecting it to boundary constraints, but have a derivative issue somewhere

* Added second version of double integrator

* Fixed an issue in the CGL partials

* Double integrator problem now solves

* exec comp for boundary constraints and objective now added to boundary eval group when necessary

* Vandermonde control interpolation is inaccurate as number of grid points increases.

* Added birkhoff test for matrix state cannonball. Fixed a bug with indexing the final boundary constraints

* Fixed an issue with indexing connections to endpoint ODE

* added birkhoff transcription tests for basic examples

* added birkhoff test case to the double integrator

* birkhoff state rates can now come from sources other than the ODE. Double integrator works but does not quite hit the correct value for final x

* Added more example problems

* Fixed configure_io in birkhoff_iter_group and a promotion issue with the boundary_ode

* fixed accuracy issue in double integrator problem

* Fixed accuracy issue of double integrator and minor bugs in the rate_source functions

* tests for hypersensitive with birkhoff

* Fixed shuttle reentry problems to work with Birkhoff transcription

* Fixed accuracy issue in double integrator problem and modified breakwell double integrator to use inputs as rate sources

* decoupled x_a and x_b from x_ab

* Decoupled initial_states:x (x_a), final_states:x (x_b), and state_segment_ends:x (x_ab = [x_a x_b].
The last is always a design variable, but it is now constrained to match x_a at the start and x_b at
the end.  This will facilitate state jumps at the end of phases.

x_ab will also be more easily turned into a segment-by-segment variable, with a size of (num_segs, 2, shape).

* Boundary ODE evaluation is now done with initial_states:x and final_states:x, instead of the endpoint values from states:x.

* cleanup

* fixed the number of segments and nodes in the aircraft cruise porblem

* works with multiple segments and scalar states, matrix cannonball test case not working yet.

* BirkhoffCollocationComp works with segments and has correct derivatives

* issue in the residual birkhoff iter group dealing with explicitly named residuals

* added a two phaseexample for birkhoff and fixed a bug in paramter connections

* Added correct inputs/outputs/resids to the birkhoff state resids comp, but currently only works with fd/cs derivs and a modified version of openmdao

* Solve segments now works once OpenMDAO issue [#3030](OpenMDAO/OpenMDAO#3030) is addressed. Sparse partials are still problematic but it works with cs/fd

* fixed an issue where input initial was still adding variables as design variables

* Added goddard rocket problem and low thrust spiral problem to examples. Only rocket problem currently works

* fixed test problems to work with IPOPT

* Fixed collocation comp tests

* formatting fixes to comply with pycodestyle

* Fixed a bug where array valued ref and ref0 were throwing an error

* fixed test_birkhoff_collocation

* added nonlinear solver to iter group. removed mux group and boundary group from check partials:

* changed optimizer back to IPOPT for examples

* Fixed moon landing problem and docstring linting

* Fixed a bug where ref and ref0 where being set to same value

* increased number of nodes on reentry problem

* fixes for new timeseries comp updates

* Removing coverage from the latest entry due to poor performance.

* More cleanup from the Timeseries component shuffle

* workflow fixes

* minimum version of openmdao bumped to 3.27

* remove coverage submission from latest

---------

Co-authored-by: Kaushik Ponnapalli <[email protected]>
Co-authored-by: kaushikponnapalli <[email protected]>

.github/workflows/dymos_tests_workflow.yml

@@ -88,7 +88,7 @@ jobs:
             PETSc: 3.13
             PYOPTSPARSE: 'v2.6.1'
             SNOPT: 7.2
-            OPENMDAO: 3.26.0
+            OPENMDAO: 3.27.0
             OPTIONAL: '[test]'
 
     steps:
@@ -312,10 +312,14 @@ jobs:
           testflo -n 1 joss/test --pre_announce
           testflo -b benchmark --pre_announce
           cd $HOME
-          testflo dymos -n 2 --show_skipped --coverage --durations 20 --coverpkg dymos
+          if [[ "${{ matrix.NAME }}" != "latest" ]]; then
+            testflo dymos -n 2 --show_skipped --durations 20 --coverage --coverpkg dymos
+          else
+            testflo dymos -n 2 --show_skipped --durations 20
+          fi
 
       - name: Submit coverage
-        if: (github.event_name != 'workflow_dispatch')
+        if: ((github.event_name != 'workflow_dispatch') && (matrix.NAME != 'latest'))
         shell: bash -l {0}
         env:
           COVERALLS_REPO_TOKEN: ${{ secrets.GITHUB_TOKEN }}

dymos/__init__.py

@@ -2,8 +2,8 @@
 
 
 from .phase import Phase, AnalyticPhase
-from .transcriptions import GaussLobatto, Radau, ExplicitShooting, Analytic
-from .transcriptions.grid_data import GaussLobattoGrid, RadauGrid, UniformGrid
+from .transcriptions import GaussLobatto, Radau, ExplicitShooting, Analytic, Birkhoff
+from .transcriptions.grid_data import GaussLobattoGrid, RadauGrid, UniformGrid, BirkhoffGrid
 from .trajectory.trajectory import Trajectory
 from .run_problem import run_problem
 from .load_case import load_case

dymos/examples/aircraft_steady_flight/test/test_aircraft_cruise.py

@@ -208,6 +208,106 @@ def test_cruise_results_radau(self):
 
         assert_near_equal(range, tas*time, tolerance=1.0E-4)
 
+    def test_cruise_results_birkhoff(self):
+        p = om.Problem(model=om.Group())
+        if optimizer == 'SNOPT':
+            p.driver = om.pyOptSparseDriver()
+            p.driver.options['optimizer'] = optimizer
+            p.driver.declare_coloring()
+            p.driver.opt_settings['Major iterations limit'] = 100
+            p.driver.opt_settings['Major step limit'] = 0.05
+            p.driver.opt_settings['Major feasibility tolerance'] = 1.0E-6
+            p.driver.opt_settings['Major optimality tolerance'] = 1.0E-6
+            p.driver.opt_settings["Linesearch tolerance"] = 0.10
+            p.driver.opt_settings['iSumm'] = 6
+            p.driver.opt_settings['Verify level'] = 3
+        else:
+            p.driver = om.ScipyOptimizeDriver()
+            p.driver.declare_coloring()
+
+        # Pass Reference Area from an external source
+        assumptions = p.model.add_subsystem('assumptions', om.IndepVarComp())
+        assumptions.add_output('S', val=427.8, units='m**2')
+        assumptions.add_output('mass_empty', val=1.0, units='kg')
+        assumptions.add_output('mass_payload', val=1.0, units='kg')
+
+        grid = dm.BirkhoffGrid(num_segments=1, nodes_per_seg=10, grid_type='lgl')
+
+        transcription = dm.Birkhoff(grid=grid)
+        phase = dm.Phase(ode_class=AircraftODE, transcription=transcription)
+        p.model.add_subsystem('phase0', phase)
+
+        phase.set_time_options(fix_initial=True,
+                               duration_bounds=(3600, 3600),
+                               duration_ref=3600)
+
+        phase.set_time_options(fix_initial=True,
+                               duration_bounds=(3600, 3600),
+                               duration_ref=3600)
+
+        phase.add_state('range', units='km', fix_initial=True, fix_final=False, scaler=0.01,
+                        rate_source='range_rate_comp.dXdt:range',
+                        defect_scaler=0.01)
+        phase.add_state('mass_fuel', units='kg', fix_final=True, upper=20000.0, lower=0.0,
+                        rate_source='propulsion.dXdt:mass_fuel',
+                        scaler=1.0E-4, defect_scaler=1.0E-2)
+        phase.add_state('alt',
+                        rate_source='climb_rate',
+                        units='km', fix_initial=True)
+
+        phase.add_control('mach', targets=['tas_comp.mach', 'aero.mach'], opt=False)
+
+        phase.add_control('climb_rate', targets=['gam_comp.climb_rate'], units='m/s', opt=False)
+
+        phase.add_parameter('S',
+                            targets=['aero.S', 'flight_equilibrium.S', 'propulsion.S'],
+                            units='m**2')
+
+        phase.add_parameter('mass_empty', targets=['mass_comp.mass_empty'], units='kg')
+        phase.add_parameter('mass_payload', targets=['mass_comp.mass_payload'], units='kg')
+
+        phase.add_path_constraint('propulsion.tau', lower=0.01, upper=1.0)
+
+        phase.add_timeseries_output('tas_comp.TAS')
+
+        p.model.connect('assumptions.S', 'phase0.parameters:S')
+        p.model.connect('assumptions.mass_empty', 'phase0.parameters:mass_empty')
+        p.model.connect('assumptions.mass_payload', 'phase0.parameters:mass_payload')
+
+        phase.add_objective('time', loc='final', ref=3600)
+
+        p.model.linear_solver = om.DirectSolver()
+
+        p.setup()
+
+        p['phase0.t_initial'] = 0.0
+        p['phase0.t_duration'] = 1.0 * 3600.0
+        p['phase0.initial_states:range'] = 0.0
+        p['phase0.final_states:mass_fuel'] = 0.0
+        p['phase0.initial_states:alt'] = 5.0
+        p['phase0.controls:mach'] = 0.8
+        p['phase0.controls:climb_rate'] = 0.0
+
+        p['assumptions.S'] = 427.8
+        p['assumptions.mass_empty'] = 0.15E6
+        p['assumptions.mass_payload'] = 84.02869 * 400
+
+        dm.run_problem(p)
+
+        time = p.get_val('phase0.timeseries.time')
+        tas = p.get_val('phase0.timeseries.TAS', units='km/s')
+        range = p.get_val('phase0.timeseries.range')
+
+        assert_near_equal(range, tas*time, tolerance=1.0E-4)
+
+        exp_out = phase.simulate()
+
+        time = exp_out.get_val('phase0.timeseries.time')
+        tas = exp_out.get_val('phase0.timeseries.TAS', units='km/s')
+        range = exp_out.get_val('phase0.timeseries.range')
+
+        assert_near_equal(range, tas*time, tolerance=1.0E-4)
+
 
 if __name__ == '__main__':  # pragma: no cover
     unittest.main()

dymos/examples/brachistochrone/test/ex_brachistochrone.py

@@ -1,4 +1,5 @@
 import matplotlib
+import numpy as np
 
 import openmdao.api as om
 from openmdao.utils.testing_utils import require_pyoptsparse
@@ -39,13 +40,19 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran
                             nodes_per_seg=transcription_order + 1,
                             compressed=compressed)
         t = dm.ExplicitShooting(grid=grid)
+    elif transcription == 'birkhoff':
+        grid = dm.BirkhoffGrid(num_segments=num_segments,
+                               nodes_per_seg=transcription_order + 1,
+                               compressed=compressed, grid_type='cgl')
+        t = dm.Birkhoff(grid=grid)
+        # phase = dm.ImplicitPhase(ode_class=BrachistochroneODE, num_nodes=11)
 
     traj = dm.Trajectory()
     phase = dm.Phase(ode_class=BrachistochroneODE, transcription=t)
     p.model.add_subsystem('traj0', traj)
     traj.add_phase('phase0', phase)
 
-    phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+    phase.set_time_options(fix_initial=True, duration_bounds=(0.5, 10))
 
     phase.add_state('x', fix_initial=True, fix_final=False, solve_segments=solve_segments)
     phase.add_state('y', fix_initial=True, fix_final=False, solve_segments=solve_segments)
@@ -68,17 +75,20 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran
 
     phase.add_boundary_constraint('x', loc='final', equals=10)
     phase.add_boundary_constraint('y', loc='final', equals=5)
+    phase.add_path_constraint('y', lower=0, upper=20)
     # Minimize time at the end of the phase
     phase.add_objective('time_phase', loc='final', scaler=10)
 
+    phase.add_timeseries_output('check')
+
     p.setup(check=['unconnected_inputs'], force_alloc_complex=force_alloc_complex)
 
     phase.set_simulate_options(method='RK23')
 
     p['traj0.phase0.t_initial'] = 0.0
     p['traj0.phase0.t_duration'] = 2.0
 
-    if transcription.startswith('shooting'):
+    if transcription.startswith('shooting') or transcription == 'birkhoff':
         p['traj0.phase0.initial_states:x'] = 0
         p['traj0.phase0.initial_states:y'] = 10
         p['traj0.phase0.initial_states:v'] = 0
@@ -90,12 +100,15 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran
     p['traj0.phase0.controls:theta'] = phase.interp('theta', [5, 100])
     p['traj0.phase0.parameters:g'] = 9.80665
 
-    dm.run_problem(p, run_driver=run_driver, simulate=True, make_plots=True)
+    dm.run_problem(p, run_driver=run_driver, simulate=True, make_plots=True,
+                   simulate_kwargs={'times_per_seg': 100})
 
     return p
 
 
 if __name__ == '__main__':
-    p = brachistochrone_min_time(transcription='gauss-lobatto', num_segments=5, run_driver=True,
-                                 transcription_order=5, compressed=False, optimizer='IPOPT',
-                                 solve_segments=False, force_alloc_complex=True)
+
+    with dm.options.temporary(include_check_partials=True):
+        p = brachistochrone_min_time(transcription='birkhoff', num_segments=1, run_driver=True,
+                                     transcription_order=19, compressed=False, optimizer='SLSQP',
+                                     solve_segments=False, force_alloc_complex=True)

dymos/examples/brachistochrone/test/ex_brachistochrone_vector_states.py

@@ -5,7 +5,7 @@
 
 
 def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, transcription_order=3,
-                             compressed=True, optimizer='SLSQP',
+                             grid_type='lgl', compressed=True, optimizer='SLSQP',
                              dynamic_simul_derivs=True, force_alloc_complex=False,
                              solve_segments=False, run_driver=True):
     p = om.Problem(model=om.Group())
@@ -33,6 +33,11 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran
                                  order=transcription_order,
                                  compressed=compressed)
         fix_final = not solve_segments
+    elif transcription == 'birkhoff':
+        gd = dm.BirkhoffGrid(num_segments=num_segments, nodes_per_seg=transcription_order,
+                             grid_type=grid_type)
+        transcription = dm.Birkhoff(grid=gd)
+        fix_final = not solve_segments
 
     traj = dm.Trajectory()
     phase = dm.Phase(ode_class=BrachistochroneVectorStatesODE,
@@ -45,7 +50,8 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran
 
     # can't fix final position if you're solving the segments
 
-    phase.add_state('pos', fix_initial=True, fix_final=fix_final, solve_segments=solve_segments, ref=[1, 1])
+    phase.add_state('pos', fix_initial=True, fix_final=fix_final,
+                    solve_segments=solve_segments, ref=[1, 1])
     #
     phase.add_state('v', fix_initial=True, fix_final=False, solve_segments=solve_segments)
     #
@@ -74,6 +80,11 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran
     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
+    if isinstance(transcription, dm.Birkhoff):
+        p['traj0.phase0.initial_states:pos'] = pos0
+        p['traj0.phase0.initial_states:v'] = 0.0
+        p['traj0.phase0.final_states:pos'] = posf
+        p['traj0.phase0.final_states:v'] = 9.9
 
     p.run_model()
     if run_driver:
@@ -85,4 +96,5 @@ def brachistochrone_min_time(transcription='gauss-lobatto', num_segments=8, tran
 if __name__ == '__main__':
     p = brachistochrone_min_time(transcription='radau-ps', num_segments=5, run_driver=True,
                                  transcription_order=5, compressed=False, optimizer='SLSQP',
-                                 solve_segments=False, force_alloc_complex=True, dynamic_simul_derivs=True)
+                                 solve_segments=False, force_alloc_complex=True,
+                                 dynamic_simul_derivs=True)

dymos/examples/brachistochrone/test/test_brachistochrone_solve_segments.py

@@ -204,6 +204,33 @@ def test_ex_brachistochrone_vs_gl_single_segment(self):
                                                            transcription_order=11)
         self.assert_results(p)
 
+    def test_ex_brachistochrone_vs_birkhoff_lgl(self):
+        ex_brachistochrone_vs.SHOW_PLOTS = False
+        p = ex_brachistochrone_vs.brachistochrone_min_time(transcription='birkhoff',
+                                                           force_alloc_complex=True,
+                                                           solve_segments=False,
+                                                           grid_type='lgl',
+                                                           num_segments=1,
+                                                           transcription_order=15,
+                                                           optimizer='SLSQP')
+
+        cpd = p.check_partials(compact_print=True, method='cs')
+        from dymos.utils.testing_utils import assert_check_partials
+        assert_check_partials(cpd)
+
+        self.assert_results(p)
+
+    def test_ex_brachistochrone_vs_birkhoff_cgl(self):
+        ex_brachistochrone_vs.SHOW_PLOTS = False
+        p = ex_brachistochrone_vs.brachistochrone_min_time(transcription='birkhoff',
+                                                           compressed=False,
+                                                           force_alloc_complex=True,
+                                                           solve_segments='forward',
+                                                           grid_type='cgl',
+                                                           num_segments=1,
+                                                           transcription_order=11)
+        self.assert_results(p)
+
 
 @require_pyoptsparse(optimizer='IPOPT')
 @use_tempdirs

dymos/examples/cannonball/test/test_cannonball_matrix_state.py

@@ -1,5 +1,6 @@
 import openmdao.api as om
 import dymos as dm
+import numpy as np
 
 import unittest
 from openmdao.utils.assert_utils import assert_near_equal
@@ -56,6 +57,9 @@ def _make_problem(self, tx):
 
         p.set_val('traj.phase.t_initial', 0)
         p.set_val('traj.phase.t_duration', 5)
+        if isinstance(tx, dm.Birkhoff):
+            p.set_val('traj.phase.initial_states:z', np.array([[0, 0], [10, 10]]))
+            p.set_val('traj.phase.final_states:z', np.array([[10, 0], [10, -10]]))
         p.set_val('traj.phase.states:z', phase.interp('z', [[[0, 0], [10, 10]], [[10, 0], [10, -10]]]))
 
         return p
@@ -99,6 +103,42 @@ def test_cannonball_matrix_state_gl(self):
         assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
         assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
 
+    @require_pyoptsparse(optimizer='IPOPT')
+    def test_cannonball_matrix_state_birkhoff_lgl(self):
+        tx = dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=1, nodes_per_seg=20, grid_type='lgl'))
+
+        p = self._make_problem(tx)
+
+        dm.run_problem(p, simulate=True)
+
+        assert_near_equal(p.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
+        assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
+        assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
+
+        c = om.CaseReader('dymos_simulation.db').get_case('final')
+
+        assert_near_equal(c.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
+        assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
+        assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
+
+    @require_pyoptsparse(optimizer='IPOPT')
+    def test_cannonball_matrix_state_birkhoff_cgl(self):
+        tx = dm.Birkhoff(grid=dm.BirkhoffGrid(num_segments=1, nodes_per_seg=20, grid_type='cgl'))
+
+        p = self._make_problem(tx)
+
+        dm.run_problem(p, simulate=True)
+
+        assert_near_equal(p.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
+        assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
+        assert_near_equal(p.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
+
+        c = om.CaseReader('dymos_simulation.db').get_case('final')
+
+        assert_near_equal(c.get_val('traj.phase.timeseries.time')[-1], 2.03873598, tolerance=1E-5)
+        assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 1], 0.0, tolerance=1E-5)
+        assert_near_equal(c.get_val('traj.phase.timeseries.z')[-1, 0, 0], 20.3873598, tolerance=1E-5)
+
     @require_pyoptsparse(optimizer='IPOPT')
     def test_cannonball_matrix_state_radau_solve_segments(self):