Repository: OpenMDAO/dymos

Issue:
Two phases can only be linked to each other once
### Description

Hi everyone, 

Lateley I have been experimenting with multi phase trajectories. I am modelling a racing car that drives a lap around a circuit, similar to the racing car example in the docs. I am splitting the lap up into two parts to be able to force the car to keep accelerating after the final corner. Let's call the first phase A and the second B. I need the end of A to be matched to the start of B. This works without issues. However, I also need the end of B to be matched to the start of A to ensure the lap is repeatable (i.e. the car doesn't start with a ridiculously high speed and starts and ends at the same position on track). This final linkage constraint seems impossible to enforce together with the first one.

### Example

Phase B (or `phase1`) is a duplicate of my 'normal' phase.

```
phase1 = phase.duplicate(transcription=dm.Radau(num_segments=15, order=6, compressed=True))
traj.add_phase('phase1',phase1)


# Link the two parts of the lap together
traj.link_phases(phases=['phase0', 'phase1'],
                vars=['u','v','t','I_FC','SoC','n','xi','U_bus'])
# Link the states at the start and end of the phase in order to ensure a continous lap
traj.link_phases(phases=['phase1','phase0'], vars = ['u','v','r','n','I_FC'],
                 locs=['final','initial'])

```
When I use this code the latter linkage constraints are not used (except for the variables that are 'new' in this case `'r'`) and their residuals do not show up when the results are printed either. So to me it seems like I can only define one linkage between these two phases. This might be expected behaviour, but as one phase is able to be linked to itself I assumed two phases could be linked to eachother multiple times (at different ends). 

Also I can't use `vars = ['*']`  in the declaration of the constraint as my integration variable is not `time` but `s` so I get an error if I try that, not a big deal but might be good to know for people that change their integration variable names.



### Dymos Version

1.10.1.dev0

### Relevant environment information

_No response_

Discussion:
@robfalck: My suspicion is that this is a bug and that the redeclaration of a `phase0 -> phase1` linkage is overwriting some information, but I'm looking into it.

Solution patch:
diff --git a/dymos/examples/oscillator/oscillator_ode.py b/dymos/examples/oscillator/oscillator_ode.py
index dac0fcd7dd..a46aac6474 100644
--- a/dymos/examples/oscillator/oscillator_ode.py
+++ b/dymos/examples/oscillator/oscillator_ode.py
@@ -36,3 +36,37 @@ def compute(self, inputs, outputs):
         f_damper = -c * v
 
         outputs['v_dot'] = (f_spring + f_damper) / m
+
+
+class OscillatorVectorODE(om.ExplicitComponent):
+    """
+    A Dymos ODE for a damped harmonic oscillator with vectorized states.
+    """
+
+    def initialize(self):
+        self.options.declare('num_nodes', types=int)
+        self.options.declare('static_params', types=bool)
+
+    def setup(self):
+        nn = self.options['num_nodes']
+        static = self.options['static_params']
+
+        # Inputs
+        self.add_input('x', val=np.ones((nn, 2)))
+        if static:
+            self.add_input('A', val=np.ones((2, 2)))
+        else:
+            self.add_input('A', val=np.ones((nn, 2, 2)))
+
+        # Output
+        self.add_output('x_dot', val=np.zeros((nn, 2)))
+
+        self.declare_partials(of='*', wrt='*', method='fd')
+
+    def compute(self, inputs, outputs):
+        A = inputs['A']
+        x = inputs['x']
+
+        static = self.options['static_params']
+        ein_sequence = 'jk, ik->ij' if static else 'ijk, ik->ij'
+        outputs['x_dot'] = np.einsum(ein_sequence, A, x)
diff --git a/dymos/trajectory/trajectory.py b/dymos/trajectory/trajectory.py
index e687c088db..99a3e9bc38 100644
--- a/dymos/trajectory/trajectory.py
+++ b/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
diff --git a/dymos/transcriptions/explicit_shooting/explicit_shooting.py b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
index a83572864e..a194f9e57f 100644
--- a/dymos/transcriptions/explicit_shooting/explicit_shooting.py
+++ b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
@@ -283,10 +283,6 @@ def configure_states(self, phase):
                                      ref0=options['ref0'],
                                      ref=options['ref'])
 
-    def _get_ode(self, phase):
-        ode = phase._get_subsystem('ode')
-        return ode
-
     def setup_ode(self, phase):
         """
         Setup the ode for this transcription.
diff --git a/dymos/transcriptions/explicit_shooting/ode_evaluation_group.py b/dymos/transcriptions/explicit_shooting/ode_evaluation_group.py
index a7caef3805..2aa64ad0ee 100644
--- a/dymos/transcriptions/explicit_shooting/ode_evaluation_group.py
+++ b/dymos/transcriptions/explicit_shooting/ode_evaluation_group.py
@@ -218,7 +218,6 @@ def _configure_states(self):
             self.add_constraint(f'state_rate_collector.state_rates:{name}_rate')
 
     def _configure_params(self):
-        vec_size = self._vec_size
         ode_inputs = get_promoted_vars(self.ode, iotypes='input', metadata_keys=['shape', 'units', 'val', 'tags'])
 
         for name, options in self._parameter_options.items():
@@ -242,13 +241,8 @@ def _configure_params(self):
             # Promote targets from the ODE
             for tgt in targets:
                 if tgt in options['static_targets']:
-                    src_idxs = None
                     shape = None
-                else:
-                    src_rows = np.zeros(vec_size, dtype=int)
-                    src_idxs = om.slicer[src_rows, ...]
-
-                self.promotes('ode', inputs=[(tgt, var_name)], src_indices=src_idxs,
+                self.promotes('ode', inputs=[(tgt, var_name)],
                               src_shape=shape)
             if targets:
                 self.set_input_defaults(name=var_name,


Test patch:
diff --git a/dymos/examples/cannonball/test/test_two_phase_cannonball_renamed_time.py b/dymos/examples/cannonball/test/test_two_phase_cannonball_renamed_time.py
new file mode 100644
index 0000000000..7756e636d0
--- /dev/null
+++ b/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()
diff --git a/dymos/examples/finite_burn_orbit_raise/test/test_two_phase_orbit_double_linked.py b/dymos/examples/finite_burn_orbit_raise/test/test_two_phase_orbit_double_linked.py
new file mode 100644
index 0000000000..9d8a5b4a46
--- /dev/null
+++ b/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)
diff --git a/dymos/examples/oscillator/test/test_oscillator_vector_states.py b/dymos/examples/oscillator/test/test_oscillator_vector_states.py
new file mode 100644
index 0000000000..924d19afad
--- /dev/null
+++ b/dymos/examples/oscillator/test/test_oscillator_vector_states.py
@@ -0,0 +1,99 @@
+import unittest
+from openmdao.utils.testing_utils import use_tempdirs
+from openmdao.utils.assert_utils import assert_near_equal
+
+import openmdao.api as om
+import dymos as dm
+import numpy as np
+
+
+# @use_tempdirs
+class TestDocOscillator(unittest.TestCase):
+
+    def test_matrix_param(self):
+
+        from dymos.examples.oscillator.oscillator_ode import OscillatorVectorODE
+
+        # Instantiate an OpenMDAO Problem instance.
+        prob = om.Problem()
+        prob.driver = om.ScipyOptimizeDriver()
+        prob.driver.options["optimizer"] = 'SLSQP'
+
+        static_params = False
+
+        t = dm.Radau(num_segments=2, order=3)
+        phase = dm.Phase(ode_class=OscillatorVectorODE, transcription=t,
+                         ode_init_kwargs={'static_params': static_params})
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(1, 2), duration_ref=1)
+        phase.add_state("x", fix_initial=True, rate_source="x_dot")
+
+        A_mat = np.array(
+            [
+                [0, 1],
+                [-1, 0]
+            ]
+        )
+
+        # argument "dynamic" doesn't seem to help
+        phase.add_parameter("A", val=A_mat, targets=["A"], static_target=static_params)
+        phase.add_objective("time", loc="final", scaler=1)
+
+        traj = dm.Trajectory()
+        traj.add_phase("phase0", phase)
+
+        prob.model.add_subsystem("traj", traj)
+
+        prob.driver.declare_coloring()
+        prob.setup(force_alloc_complex=True)
+        phase.set_state_val('x', vals=[[1, 0], [1, 0]])
+
+        dm.run_problem(prob, run_driver=True, simulate=True, make_plots=True)
+        t_f = prob.get_val('traj.phase0.timeseries.time')[-1]
+        final_state = prob.get_val('traj.phase0.timeseries.x')[-1, :]
+        assert_near_equal(final_state, np.array([np.cos(t_f), -np.sin(t_f)]).ravel(),
+                          tolerance=1e-5)
+
+    def test_matrix_static_param(self):
+
+        from dymos.examples.oscillator.oscillator_ode import OscillatorVectorODE
+
+        # Instantiate an OpenMDAO Problem instance.
+        prob = om.Problem()
+        prob.driver = om.ScipyOptimizeDriver()
+        prob.driver.options["optimizer"] = 'SLSQP'
+
+        static_params = True
+
+        t = dm.Radau(num_segments=2, order=3)
+        phase = dm.Phase(ode_class=OscillatorVectorODE, transcription=t,
+                         ode_init_kwargs={'static_params': static_params})
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(1, 2), duration_ref=1)
+        phase.add_state("x", fix_initial=True, rate_source="x_dot")
+
+        A_mat = np.array(
+            [
+                [0, 1],
+                [-1, 0]
+            ]
+        )
+
+        # argument "dynamic" doesn't seem to help
+        phase.add_parameter("A", val=A_mat, targets=["A"], static_target=static_params)
+        phase.add_objective("time", loc="final", scaler=1)
+
+        traj = dm.Trajectory()
+        traj.add_phase("phase0", phase)
+
+        prob.model.add_subsystem("traj", traj)
+
+        prob.driver.declare_coloring()
+        prob.setup(force_alloc_complex=True)
+        phase.set_state_val('x', vals=[[1, 0], [1, 0]])
+
+        dm.run_problem(prob, run_driver=True, simulate=True, make_plots=True)
+        t_f = prob.get_val('traj.phase0.timeseries.time')[-1]
+        final_state = prob.get_val('traj.phase0.timeseries.x')[-1, :]
+        assert_near_equal(final_state, np.array([np.cos(t_f), -np.sin(t_f)]).ravel(),
+                          tolerance=1e-5)
diff --git a/dymos/trajectory/test/test_t_initial_bounds.py b/dymos/trajectory/test/test_t_initial_bounds.py
index 31db8ab87c..970fef9527 100644
--- a/dymos/trajectory/test/test_t_initial_bounds.py
+++ b/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': {