Allow double linkage between phases (#1121)

* Changed linkage constraint to use the name of the time variable specified instead of default. Added a test

* Changed behavior when a cyclical link is detected across phases. Now raises a warning instead of an error

* modified how the warning was issued to be non-interrupting

* fixed indentation

* removed warning for cycles

* removed unused import

dymos/examples/cannonball/test/test_two_phase_cannonball_renamed_time.py

@@ -0,0 +1,112 @@
+import unittest
+
+import openmdao.api as om
+import dymos as dm
+
+from openmdao.utils.assert_utils import assert_near_equal
+from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
+
+from dymos.examples.cannonball.size_comp import CannonballSizeComp
+from dymos.examples.cannonball.cannonball_ode import CannonballODE
+
+
+@use_tempdirs
+class TestTwoPhaseCannonballRenamedTime(unittest.TestCase):
+
+    @require_pyoptsparse(optimizer='SLSQP')
+    def test_rename_time(self):
+        p = om.Problem(model=om.Group())
+
+        p.driver = om.pyOptSparseDriver()
+        p.driver.options['optimizer'] = 'SLSQP'
+        p.driver.declare_coloring()
+
+        p.model.add_subsystem('size_comp', CannonballSizeComp(),
+                              promotes_inputs=['radius', 'dens'])
+        p.model.set_input_defaults('dens', val=7.87, units='g/cm**3')
+        p.model.add_design_var('radius', lower=0.01, upper=0.10,
+                               ref0=0.01, ref=0.10, units='m')
+
+        traj = p.model.add_subsystem('traj', dm.Trajectory())
+
+        transcription = dm.Radau(num_segments=5, order=3, compressed=True)
+        ascent = dm.Phase(ode_class=CannonballODE, transcription=transcription)
+
+        ascent = traj.add_phase('ascent', ascent)
+
+        ascent.set_time_options(fix_initial=True, duration_bounds=(1, 100),
+                                duration_ref=100, units='s', name='t')
+
+        ascent.set_state_options('r', rate_source='r_dot', fix_initial=True, fix_final=False)
+        ascent.set_state_options('h', rate_source='h_dot', fix_initial=True, fix_final=False)
+        ascent.set_state_options('gam', rate_source='gam_dot', fix_initial=False, fix_final=True)
+        ascent.set_state_options('v', rate_source='v_dot', fix_initial=False, fix_final=False)
+
+        ascent.add_parameter('S', units='m**2', static_target=True)
+        ascent.add_parameter('m', units='kg', static_target=True)
+
+        ascent.add_boundary_constraint('ke', loc='initial',
+                                       upper=400000, lower=0, ref=100000)
+
+        # Second Phase (descent)
+        transcription = dm.GaussLobatto(num_segments=5, order=3, compressed=True)
+        descent = dm.Phase(ode_class=CannonballODE, transcription=transcription)
+
+        traj.add_phase('descent', descent)
+
+        descent.set_time_options(initial_bounds=(.5, 100), duration_bounds=(.5, 100),
+                                 duration_ref=100, units='s', name='t')
+
+        descent.add_state('r', rate_source='r_dot')
+        descent.add_state('h', rate_source='h_dot', fix_initial=False, fix_final=True)
+        descent.add_state('gam', rate_source='gam_dot', fix_initial=False, fix_final=False)
+        descent.add_state('v', rate_source='v_dot', fix_initial=False, fix_final=False)
+
+        descent.add_parameter('S', units='m**2', static_target=True)
+        descent.add_parameter('m', units='kg', static_target=True)
+
+        descent.add_objective('r', loc='final', scaler=-1.0)
+
+        traj.add_parameter('CD',
+                           targets={'ascent': ['CD'], 'descent': ['CD']},
+                           val=0.5, units=None, opt=False, static_target=True)
+        traj.add_parameter('m', units='kg', val=1.0,
+                           targets={'ascent': 'm', 'descent': 'm'}, static_target=True)
+        traj.add_parameter('S', units='m**2', val=0.005, static_target=True)
+
+        traj.link_phases(phases=['ascent', 'descent'], vars=['*'])
+
+        p.model.connect('size_comp.mass', 'traj.parameters:m')
+        p.model.connect('size_comp.S', 'traj.parameters:S')
+
+        p.model.linear_solver = om.DirectSolver()
+
+        p.setup()
+
+        p.set_val('radius', 0.05, units='m')
+        p.set_val('dens', 7.87, units='g/cm**3')
+
+        traj.set_parameter_val('CD', 0.5)
+
+        ascent.set_time_val(initial=0.0, duration=10.0)
+
+        ascent.set_state_val('r', [0, 100])
+        ascent.set_state_val('h', [0, 100])
+        ascent.set_state_val('v', [200, 150])
+        ascent.set_state_val('gam', [25, 0], units='deg')
+
+        descent.set_time_val(initial=10.0, duration=10.0)
+
+        descent.set_state_val('r', [100, 200])
+        descent.set_state_val('h', [100, 0])
+        descent.set_state_val('v', [150, 200])
+        descent.set_state_val('gam', [0, -45], units='deg')
+
+        dm.run_problem(p, simulate=True)
+
+        assert_near_equal(p.get_val('traj.descent.states:r')[-1],
+                          3183.25, tolerance=1.0E-2)
+
+
+if __name__ == '__main__':  # pragma: no cover
+    unittest.main()

dymos/examples/finite_burn_orbit_raise/test/test_two_phase_orbit_double_linked.py

@@ -0,0 +1,91 @@
+import unittest
+import numpy as np
+import openmdao.api as om
+import dymos as dm
+
+from dymos.examples.finite_burn_orbit_raise.finite_burn_eom import FiniteBurnODE
+
+from openmdao.utils.testing_utils import use_tempdirs, require_pyoptsparse
+from openmdao.utils.assert_utils import assert_near_equal
+
+
+@require_pyoptsparse(optimizer='SLSQP')
+@use_tempdirs
+class TestTwoPhaseOrbitDoubleLinked(unittest.TestCase):
+
+    def test_two_phase_orbit_double_linked(self):
+        traj = dm.Trajectory()
+        p = om.Problem(model=om.Group())
+        p.model.add_subsystem('traj', traj)
+
+        p.driver = om.pyOptSparseDriver()
+        p.driver.options['optimizer'] = 'SLSQP'
+
+        p.driver.declare_coloring()
+
+        traj.add_parameter('c', opt=False, val=0.0, units='DU/TU',
+                           targets={'phase0': ['c']})
+
+        traj.add_parameter('accel', opt=False, val=0.0, units='DU/TU**2',
+                           targets={'phase0': ['accel']})
+
+        traj.add_parameter('u1', opt=False, val=0.0, units='deg',
+                           targets={'phase0': ['u1']})
+
+        # First half of the orbit
+        phase0 = dm.Phase(ode_class=FiniteBurnODE,
+                          transcription=dm.GaussLobatto(num_segments=10, order=3, compressed=True))
+
+        phase0 = traj.add_phase('phase0', phase0)
+
+        phase0.set_time_options(fix_initial=True, fix_duration=True, units='TU')
+        phase0.add_state('r', fix_initial=False, fix_final=False,
+                         rate_source='r_dot', targets=['r'], units='DU')
+        phase0.add_state('theta', fix_initial=True, fix_final=False,
+                         rate_source='theta_dot', targets=['theta'], units='rad')
+        phase0.add_state('vr', fix_initial=False, fix_final=False,
+                         rate_source='vr_dot', targets=['vr'], units='DU/TU')
+        phase0.add_state('vt', fix_initial=False, fix_final=False,
+                         rate_source='vt_dot', targets=['vt'], units='DU/TU')
+
+        phase0.add_objective('time', loc='final')
+
+        # Second half of the orbit
+        phase1 = dm.Phase(ode_class=FiniteBurnODE,
+                          transcription=dm.GaussLobatto(num_segments=10, order=3, compressed=True))
+
+        phase1 = traj.add_phase('phase1', phase1)
+
+        phase1.set_time_options(fix_initial=True, fix_duration=True, units='TU')
+        phase1.add_state('r', fix_initial=False, fix_final=False,
+                         rate_source='r_dot', targets=['r'], units='DU')
+        phase1.add_state('theta', fix_initial=False, fix_final=False,
+                         rate_source='theta_dot', targets=['theta'], units='rad')
+        phase1.add_state('vr', fix_initial=False, fix_final=False,
+                         rate_source='vr_dot', targets=['vr'], units='DU/TU')
+        phase1.add_state('vt', fix_initial=False, fix_final=False,
+                         rate_source='vt_dot', targets=['vt'], units='DU/TU')
+
+        traj.link_phases(phases=['phase0', 'phase1'], vars=['r', 'theta', 'vr', 'vt'], connected=False)
+        traj.link_phases(phases=['phase1', 'phase0'], vars=['r', 'vr', 'vt'], connected=False)
+
+        p.setup()
+        phase0.set_time_val(initial=0.0, duration=np.pi)
+        phase0.set_state_val('r', vals=[2, 1])
+        phase0.set_state_val('theta', vals=[0, np.pi])
+        phase0.set_state_val('vr', vals=[0, 0])
+        phase0.set_state_val('vt', vals=[2, 1])
+
+        phase1.set_time_val(initial=np.pi, duration=np.pi)
+        phase0.set_state_val('r', vals=[1, 1])
+        phase0.set_state_val('theta', vals=[0, np.pi])
+        phase0.set_state_val('vr', vals=[0, 0])
+        phase0.set_state_val('vt', vals=[1, 1])
+
+        dm.run_problem(p, make_plots=True)
+
+        r_1 = p.get_val('traj.phase0.timeseries.r')
+        theta_2 = p.get_val('traj.phase1.timeseries.theta')
+
+        assert_near_equal(r_1[0], 1.0, tolerance=1e-9)
+        assert_near_equal(theta_2[-1], 2*np.pi, tolerance=1e-9)

dymos/trajectory/test/test_t_initial_bounds.py

@@ -5,6 +5,7 @@
 import dymos as dm
 
 from openmdao.utils.testing_utils import use_tempdirs
+from openmdao.utils.assert_utils import assert_warning
 
 
 def make_phases(traj, time_option_kwargs):
@@ -98,30 +99,6 @@ def try_model(self, probname, kwargs, conns):
 
         dm.run_problem(prob, run_driver=False)
 
-    def test_phase_link_cycle(self):
-        nphases = 3
-
-        kwargs = {}
-        conns = []
-        for i in range(nphases):
-            pname = f'phase{i}'
-            kwargs[pname] = dct = {}
-            dct['fix_initial'] = True
-            dct['fix_duration'] = False
-            dct['initial_val'] = 5. * i if i < nphases - 1 else 5. * (i - 1)
-
-            if i > 0:
-                conns.append((f'phase{i-1}', [pname]))
-
-        conns.append((pname, ['phase0']))
-
-        with self.assertRaises(Exception) as cm:
-            self.try_model('phase_link_cycle', kwargs, conns)
-
-        msg = ("'traj' <class Trajectory>: The following cycles were found in the phase "
-               "linkage graph: [['phase0', 'phase1', 'phase2']].")
-        self.assertEqual(cm.exception.args[0], msg)
-
     def test_pair_fixed_t_initial_below(self):
         kwargs = {
             'phase0': {

dymos/trajectory/trajectory.py

@@ -5,7 +5,6 @@
 import sys
 
 from openmdao.utils.om_warnings import warn_deprecation
-from openmdao.utils.graph_utils import get_sccs_topo
 from openmdao.utils.units import unit_conversion
 
 import numpy as np
@@ -711,8 +710,9 @@ def _expand_star_linkage_configure(self):
             var_pair = ('*', '*')
             if var_pair in var_dict:
                 options = var_dict[var_pair]
+                time_name = phase_b.time_options['name']
                 self.add_linkage_constraint(phase_name_a, phase_name_b,
-                                            var_a='time', var_b='time',
+                                            var_a=time_name, var_b=time_name,
                                             loc_a=options['loc_a'], loc_b=options['loc_b'],
                                             mult_a=options['mult_a'], mult_b=options['mult_b'],
                                             connected=options['connected'])
@@ -981,13 +981,6 @@ def get_final_tbounds(phase_name, node_data, old_tmin=-INF_BOUND, old_tmax=INF_B
 
         phase_graph = self._phase_graph
 
-        # since we have a graph, do a quick check that we have no cycles
-        sccs = get_sccs_topo(phase_graph)
-        cycles = sorted([s for s in sccs if len(s) > 1], key=lambda x: len(x))
-        if cycles:
-            raise RuntimeError(f"{self.msginfo}: The following cycles were found in the phase "
-                               f"linkage graph: {[sorted(c) for c in cycles]}.")
-
         node_data = phase_graph.nodes(data=True)
 
         # only keep the part of the graph where 't' connects phases