Repository: OpenMDAO/dymos (python)

Issue:
Overlapping constraints in Birkhoff transcription
### Description

Adding a boundary constraint and a path constraint for variables of any type other than state causes the error.

The lines of code when added to the brachistochrone problem throw the error

`RuntimeError: <model> <class Group>: Indices for aliases ['traj.phases.phase0->final_boundary_constraint->theta', 'traj.phases.phase0->path_constraint->theta'] are overlapping constraint/objective 'traj.phases.phase0.timeseries.timeseries_comp.theta'.`


### Example

phase.add_path_constraint('theta', lower=0.01, upper=179.9)
phase.add_boundary_constraint('theta', loc='final', lower=0.01, upper=179.9)


### Dymos Version

1.9.2.dev0

### Relevant environment information

_No response_

Pull request:
Boundary constraints in Birkhoff transcription changed to pull from boundary ODE
### Summary

Added additional introspection to birkhoff constraints so that boundary constraints on all variables now pull from the boundary ODE. Previously only states did so.
Added a test for controls and ODE outputs

### Related Issues

- Resolves #1017

### Backwards incompatibilities

None

### New Dependencies

None


Patch:
diff --git a/dymos/utils/introspection.py b/dymos/utils/introspection.py
index 73497d24..775fa29b 100644
--- a/dymos/utils/introspection.py
+++ b/dymos/utils/introspection.py
@@ -227,7 +227,10 @@ def _configure_constraint_introspection(phase):
 
                 con['shape'] = control_shape
                 con['units'] = control_units if con['units'] is None else con['units']
-                con['constraint_path'] = f'timeseries.{prefix}{var}'
+                if birkhoff and constraint_type in ('initial', 'final'):
+                    con['constraint_path'] = f'boundary_vals.{var}'
+                else:
+                    con['constraint_path'] = f'timeseries.{prefix}{var}'
 
             elif var_type in ['indep_polynomial_control', 'input_polynomial_control']:
                 prefix = 'polynomial_controls:' if phase.timeseries_options['use_prefix'] else ''
@@ -235,7 +238,10 @@ def _configure_constraint_introspection(phase):
                 control_units = phase.polynomial_control_options[var]['units']
                 con['shape'] = control_shape
                 con['units'] = control_units if con['units'] is None else con['units']
-                con['constraint_path'] = f'timeseries.{prefix}{var}'
+                if birkhoff and constraint_type in ('initial', 'final'):
+                    con['constraint_path'] = f'boundary_vals.{var}'
+                else:
+                    con['constraint_path'] = f'timeseries.{prefix}{var}'
 
             elif var_type == 'control_rate':
                 prefix = 'control_rates:' if phase.timeseries_options['use_prefix'] else ''
@@ -245,7 +251,10 @@ def _configure_constraint_introspection(phase):
                 con['shape'] = control_shape
                 con['units'] = get_rate_units(control_units, time_units, deriv=1) \
                     if con['units'] is None else con['units']
-                con['constraint_path'] = f'timeseries.{prefix}{var}'
+                if birkhoff and constraint_type in ('initial', 'final'):
+                    con['constraint_path'] = f'boundary_vals.{var}'
+                else:
+                    con['constraint_path'] = f'timeseries.{prefix}{var}'
 
             elif var_type == 'control_rate2':
                 prefix = 'control_rates:' if phase.timeseries_options['use_prefix'] else ''
@@ -255,7 +264,10 @@ def _configure_constraint_introspection(phase):
                 con['shape'] = control_shape
                 con['units'] = get_rate_units(control_units, time_units, deriv=2) \
                     if con['units'] is None else con['units']
-                con['constraint_path'] = f'timeseries.{prefix}{var}'
+                if birkhoff and constraint_type in ('initial', 'final'):
+                    con['constraint_path'] = f'boundary_vals.{var}'
+                else:
+                    con['constraint_path'] = f'timeseries.{prefix}{var}'
 
             elif var_type == 'polynomial_control_rate':
                 prefix = 'polynomial_control_rates:' if phase.timeseries_options['use_prefix'] else ''
@@ -265,7 +277,10 @@ def _configure_constraint_introspection(phase):
                 con['shape'] = control_shape
                 con['units'] = get_rate_units(control_units, time_units, deriv=1) \
                     if con['units'] is None else con['units']
-                con['constraint_path'] = f'timeseries.{prefix}{var}'
+                if birkhoff and constraint_type in ('initial', 'final'):
+                    con['constraint_path'] = f'boundary_vals.{var}'
+                else:
+                    con['constraint_path'] = f'timeseries.{prefix}{var}'
 
             elif var_type == 'polynomial_control_rate2':
                 prefix = 'polynomial_control_rates:' if phase.timeseries_options['use_prefix'] else ''
@@ -275,7 +290,10 @@ def _configure_constraint_introspection(phase):
                 con['shape'] = control_shape
                 con['units'] = get_rate_units(control_units, time_units, deriv=2) \
                     if con['units'] is None else con['units']
-                con['constraint_path'] = f'timeseries.{prefix}{var}'
+                if birkhoff and constraint_type in ('initial', 'final'):
+                    con['constraint_path'] = f'boundary_vals.{var}'
+                else:
+                    con['constraint_path'] = f'timeseries.{prefix}{var}'
 
             elif var_type == 'timeseries_exec_comp_output':
                 con['shape'] = (1,)
@@ -290,7 +308,11 @@ def _configure_constraint_introspection(phase):
 
                 con['shape'] = meta['shape']
                 con['units'] = meta['units']
-                con['constraint_path'] = f'timeseries.{con["constraint_name"]}'
+
+                if birkhoff and constraint_type in ('initial', 'final'):
+                    con['constraint_path'] = f'boundary_vals.{var}'
+                else:
+                    con['constraint_path'] = f'timeseries.{con["constraint_name"]}'
 
 
 def configure_controls_introspection(control_options, ode, time_units='s'):


Test patch:
diff --git a/dymos/examples/brachistochrone/test/test_brachistochrone_birkhoff_constraints.py b/dymos/examples/brachistochrone/test/test_brachistochrone_birkhoff_constraints.py
new file mode 100644
index 00000000..f47c5cfe
--- /dev/null
+++ b/dymos/examples/brachistochrone/test/test_brachistochrone_birkhoff_constraints.py
@@ -0,0 +1,251 @@
+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
+from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE
+
+
+@use_tempdirs
+class TestBrachistochroneBirkhoffConstraints(unittest.TestCase):
+
+    def test_brachistochrone_control_prefix(self):
+
+        p = om.Problem(model=om.Group())
+
+        p.driver = om.ScipyOptimizeDriver()
+        p.driver.options['optimizer'] = 'SLSQP'
+        p.driver.declare_coloring(tol=1.0E-12)
+
+        grid = dm.BirkhoffGrid(num_segments=1, nodes_per_seg=25, grid_type='lgl')
+        tx = dm.Birkhoff(grid=grid)
+
+        traj = dm.Trajectory()
+        phase = dm.Phase(ode_class=BrachistochroneODE, transcription=tx)
+        phase.timeseries_options['use_prefix'] = True
+        p.model.add_subsystem('traj', traj)
+        traj.add_phase('phase0', phase)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+
+        phase.add_state('x', fix_initial=True, fix_final=False)
+        phase.add_state('y', fix_initial=True, fix_final=False)
+
+        # Note that by omitting the targets here Dymos will automatically attempt to connect
+        # to a top-level input named 'v' in the ODE, and connect to nothing if it's not found.
+        phase.add_state('v', fix_initial=True, fix_final=False)
+
+        phase.add_control('theta',
+                          continuity=True, rate_continuity=True,
+                          units='deg')
+
+        phase.add_parameter('g', targets=['g'], units='m/s**2')
+
+        phase.add_path_constraint('theta', lower=0.01, upper=179.9)
+        phase.add_boundary_constraint('theta', loc='final', lower=0.01, upper=179.9)
+
+        phase.add_boundary_constraint('x', loc='final', equals=10)
+        phase.add_boundary_constraint('y', loc='final', equals=5)
+        # Minimize time at the end of the phase
+        phase.add_objective('time_phase', loc='final', scaler=10)
+
+        p.set_solver_print(0)
+
+        p.setup()
+
+        p['traj.phase0.t_initial'] = 0.0
+        p['traj.phase0.t_duration'] = 1.5
+
+        p['traj.phase0.initial_states:x'] = 0.0
+        p['traj.phase0.initial_states:y'] = 10.0
+        p['traj.phase0.initial_states:v'] = 0.0
+
+        p['traj.phase0.states:x'] = phase.interp('x', [0, 10])
+        p['traj.phase0.states:y'] = phase.interp('y', [10, 5])
+        p['traj.phase0.states:v'] = phase.interp('v', [0, 9.9])
+        p['traj.phase0.controls:theta'] = phase.interp('theta', [5, 100])
+        p['traj.phase0.parameters:g'] = 9.80665
+
+        p.run_driver()
+        assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-4)
+
+    def test_brachistochrone_control_no_prefix(self):
+
+        p = om.Problem(model=om.Group())
+
+        p.driver = om.ScipyOptimizeDriver()
+        p.driver.options['optimizer'] = 'SLSQP'
+        p.driver.declare_coloring(tol=1.0E-12)
+
+        grid = dm.BirkhoffGrid(num_segments=1, nodes_per_seg=25, grid_type='lgl')
+        tx = dm.Birkhoff(grid=grid)
+
+        traj = dm.Trajectory()
+        phase = dm.Phase(ode_class=BrachistochroneODE, transcription=tx)
+        phase.timeseries_options['use_prefix'] = False
+        p.model.add_subsystem('traj', traj)
+        traj.add_phase('phase0', phase)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+
+        phase.add_state('x', fix_initial=True, fix_final=False)
+        phase.add_state('y', fix_initial=True, fix_final=False)
+
+        # Note that by omitting the targets here Dymos will automatically attempt to connect
+        # to a top-level input named 'v' in the ODE, and connect to nothing if it's not found.
+        phase.add_state('v', fix_initial=True, fix_final=False)
+
+        phase.add_control('theta',
+                          continuity=True, rate_continuity=True,
+                          units='deg')
+
+        phase.add_parameter('g', targets=['g'], units='m/s**2')
+
+        phase.add_path_constraint('theta', lower=0.01, upper=179.9)
+        phase.add_boundary_constraint('theta', loc='final', lower=0.01, upper=179.9)
+
+        phase.add_boundary_constraint('x', loc='final', equals=10)
+        phase.add_boundary_constraint('y', loc='final', equals=5)
+        # Minimize time at the end of the phase
+        phase.add_objective('time_phase', loc='final', scaler=10)
+
+        p.set_solver_print(0)
+
+        p.setup()
+
+        p['traj.phase0.t_initial'] = 0.0
+        p['traj.phase0.t_duration'] = 1.5
+
+        p['traj.phase0.initial_states:x'] = 0.0
+        p['traj.phase0.initial_states:y'] = 10.0
+        p['traj.phase0.initial_states:v'] = 0.0
+
+        p['traj.phase0.states:x'] = phase.interp('x', [0, 10])
+        p['traj.phase0.states:y'] = phase.interp('y', [10, 5])
+        p['traj.phase0.states:v'] = phase.interp('v', [0, 9.9])
+        p['traj.phase0.controls:theta'] = phase.interp('theta', [5, 100])
+        p['traj.phase0.parameters:g'] = 9.80665
+
+        p.run_driver()
+        assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-4)
+
+    def test_brachistochrone_ode_prefix(self):
+
+        p = om.Problem(model=om.Group())
+
+        p.driver = om.ScipyOptimizeDriver()
+        p.driver.options['optimizer'] = 'SLSQP'
+        p.driver.declare_coloring(tol=1.0E-12)
+
+        grid = dm.BirkhoffGrid(num_segments=1, nodes_per_seg=25, grid_type='lgl')
+        tx = dm.Birkhoff(grid=grid)
+
+        traj = dm.Trajectory()
+        phase = dm.Phase(ode_class=BrachistochroneODE, transcription=tx)
+        phase.timeseries_options['use_prefix'] = True
+        p.model.add_subsystem('traj', traj)
+        traj.add_phase('phase0', phase)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+
+        phase.add_state('x', fix_initial=True, fix_final=False)
+        phase.add_state('y', fix_initial=True, fix_final=False)
+
+        # Note that by omitting the targets here Dymos will automatically attempt to connect
+        # to a top-level input named 'v' in the ODE, and connect to nothing if it's not found.
+        phase.add_state('v', fix_initial=True, fix_final=False)
+
+        phase.add_control('theta',
+                          continuity=True, rate_continuity=True,
+                          units='deg', lower=0.01, upper=179.9)
+
+        phase.add_parameter('g', targets=['g'], units='m/s**2')
+
+        phase.add_boundary_constraint('x', loc='final', equals=10)
+        phase.add_boundary_constraint('y', loc='final', equals=5)
+        phase.add_boundary_constraint('check', loc='final', lower=-50, upper=50)
+        phase.add_path_constraint('check', upper=100, lower=-100)
+        # Minimize time at the end of the phase
+        phase.add_objective('time_phase', loc='final', scaler=10)
+
+        p.set_solver_print(0)
+
+        p.setup()
+
+        p['traj.phase0.t_initial'] = 0.0
+        p['traj.phase0.t_duration'] = 1.5
+
+        p['traj.phase0.initial_states:x'] = 0.0
+        p['traj.phase0.initial_states:y'] = 10.0
+        p['traj.phase0.initial_states:v'] = 0.0
+
+        p['traj.phase0.states:x'] = phase.interp('x', [0, 10])
+        p['traj.phase0.states:y'] = phase.interp('y', [10, 5])
+        p['traj.phase0.states:v'] = phase.interp('v', [0, 9.9])
+        p['traj.phase0.controls:theta'] = phase.interp('theta', [5, 100])
+        p['traj.phase0.parameters:g'] = 9.80665
+
+        p.run_driver()
+        assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-4)
+
+    def test_brachistochrone_ode_no_prefix(self):
+
+        p = om.Problem(model=om.Group())
+
+        p.driver = om.ScipyOptimizeDriver()
+        p.driver.options['optimizer'] = 'SLSQP'
+        p.driver.declare_coloring(tol=1.0E-12)
+
+        grid = dm.BirkhoffGrid(num_segments=1, nodes_per_seg=25, grid_type='lgl')
+        tx = dm.Birkhoff(grid=grid)
+
+        traj = dm.Trajectory()
+        phase = dm.Phase(ode_class=BrachistochroneODE, transcription=tx)
+        phase.timeseries_options['use_prefix'] = False
+        p.model.add_subsystem('traj', traj)
+        traj.add_phase('phase0', phase)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+
+        phase.add_state('x', fix_initial=True, fix_final=False)
+        phase.add_state('y', fix_initial=True, fix_final=False)
+
+        # Note that by omitting the targets here Dymos will automatically attempt to connect
+        # to a top-level input named 'v' in the ODE, and connect to nothing if it's not found.
+        phase.add_state('v', fix_initial=True, fix_final=False)
+
+        phase.add_control('theta',
+                          continuity=True, rate_continuity=True,
+                          units='deg', lower=0.01, upper=179.9)
+
+        phase.add_parameter('g', targets=['g'], units='m/s**2')
+
+        phase.add_boundary_constraint('x', loc='final', equals=10)
+        phase.add_boundary_constraint('y', loc='final', equals=5)
+        phase.add_boundary_constraint('check', loc='final', lower=-50, upper=50)
+        phase.add_path_constraint('check', upper=100, lower=-100)
+
+        # Minimize time at the end of the phase
+        phase.add_objective('time_phase', loc='final', scaler=10)
+
+        p.set_solver_print(0)
+
+        p.setup()
+
+        p['traj.phase0.t_initial'] = 0.0
+        p['traj.phase0.t_duration'] = 1.5
+
+        p['traj.phase0.initial_states:x'] = 0.0
+        p['traj.phase0.initial_states:y'] = 10.0
+        p['traj.phase0.initial_states:v'] = 0.0
+
+        p['traj.phase0.states:x'] = phase.interp('x', [0, 10])
+        p['traj.phase0.states:y'] = phase.interp('y', [10, 5])
+        p['traj.phase0.states:v'] = phase.interp('v', [0, 9.9])
+        p['traj.phase0.controls:theta'] = phase.interp('theta', [5, 100])
+        p['traj.phase0.parameters:g'] = 9.80665
+
+        p.run_driver()
+        assert_near_equal(p.get_val('traj.phase0.timeseries.time')[-1], 1.8016, tolerance=1.0E-4)