Repository: OpenMDAO/dymos (python)

Issue:
Default val set in `add_control` or `set_control_options` is not applied.
### Description

Specifying the option `val=...` in `add_control` or `set_control_options` has no impact in dymos 1.13.0.


### Example

The test case to see if default values are "sticking" is no longer being applied, because the function upon which it relies overrides default values for other tests.  Removing those set_xxx_val calls will expose the error.

### Dymos Version

3.13.1-dev

### Relevant environment information

_No response_

Pull request description:
Fixed an issue where default values in add_control or set_control_options were no longer being respected.
### Summary

The default control values were set using a call to set_input_defaults, which went away because the ControlInterpComp is no longer housed in a ControlGroup. Values are now applied at add_input, and the test to make sure these default values stick has been fixed.

### Related Issues

- Resolves #1166

### Backwards incompatibilities

None

### New Dependencies

None


Patch:
diff --git a/dymos/transcriptions/common/control_comp.py b/dymos/transcriptions/common/control_comp.py
index 74a1449b..ca4a18c0 100644
--- a/dymos/transcriptions/common/control_comp.py
+++ b/dymos/transcriptions/common/control_comp.py
@@ -5,7 +5,7 @@ from scipy.linalg import block_diag
 import openmdao.api as om
 
 from ..grid_data import GridData
-from ...utils.misc import get_rate_units, CoerceDesvar
+from ...utils.misc import get_rate_units, CoerceDesvar, reshape_val
 from ...utils.lgl import lgl
 from ...utils.lagrange import lagrange_matrices
 from ...utils.indexing import get_desvar_indices
@@ -95,15 +95,16 @@ class ControlInterpComp(om.ExplicitComponent):
             if 'control_type' not in options:
                 options['control_type'] = 'full'
             if options['control_type'] == 'polynomial':
-                disc_nodes, _ = lgl(options['order'] + 1)
+                num_input_nodes = options['order'] + 1
+                disc_nodes, _ = lgl(num_input_nodes)
                 num_control_input_nodes = len(disc_nodes)
                 shape = options['shape']
                 size = np.prod(shape)
                 units = options['units']
                 rate_units = get_rate_units(units, self.options['time_units'], deriv=1)
                 rate2_units = get_rate_units(units, self.options['time_units'], deriv=2)
+                default_val = reshape_val(options['val'], shape, num_input_nodes)
 
-                input_shape = (num_control_input_nodes,) + shape
                 output_shape = (num_output_nodes,) + shape
 
                 L_de, D_de = lagrange_matrices(disc_nodes, eval_nodes)
@@ -117,7 +118,7 @@ class ControlInterpComp(om.ExplicitComponent):
                 self._output_rate_names[name] = f'control_rates:{name}_rate'
                 self._output_rate2_names[name] = f'control_rates:{name}_rate2'
 
-                self.add_input(self._input_names[name], val=np.ones(input_shape), units=units)
+                self.add_input(self._input_names[name], val=default_val, units=units)
                 self.add_output(self._output_val_names[name], shape=output_shape, units=units)
                 self.add_output(self._output_rate_names[name], shape=output_shape, units=rate_units)
                 self.add_output(self._output_rate2_names[name], shape=output_shape, units=rate2_units)
@@ -178,14 +179,14 @@ class ControlInterpComp(om.ExplicitComponent):
                 self._output_rate2_names[name] = f'control_rates:{name}_rate2'
                 shape = options['shape']
                 num_control_input_nodes = gd.subset_num_nodes['control_input']
-                input_shape = (num_control_input_nodes,) + shape
                 output_shape = (num_output_nodes,) + shape
+                default_val = reshape_val(options['val'], shape, num_control_input_nodes)
 
                 units = options['units']
                 rate_units = get_rate_units(units, time_units)
                 rate2_units = get_rate_units(units, time_units, deriv=2)
 
-                self.add_input(self._input_names[name], val=np.ones(input_shape), units=units)
+                self.add_input(self._input_names[name], val=default_val, units=units)
 
                 self.add_output(self._output_val_names[name], shape=output_shape, units=units)
 
@@ -244,74 +245,6 @@ class ControlInterpComp(om.ExplicitComponent):
             if options['control_type'] == 'polynomial':
                 num_input_nodes = options['order'] + 1
                 shape = options['shape']
-                # default_val = reshape_val(options['val'], shape, num_input_nodes)
-                if options['opt']:
-
-                    desvar_indices = np.arange(num_input_nodes, dtype=int)
-                    if options['fix_initial']:
-                        desvar_indices = desvar_indices[1:]
-                    if options['fix_final']:
-                        desvar_indices = desvar_indices[:-1]
-
-                    lb = -INF_BOUND if options['lower'] is None else options['lower']
-                    ub = INF_BOUND if options['upper'] is None else options['upper']
-
-                    self.add_design_var(f'controls:{name}',
-                                        lower=lb,
-                                        upper=ub,
-                                        ref=options['ref'],
-                                        ref0=options['ref0'],
-                                        adder=options['adder'],
-                                        scaler=options['scaler'],
-                                        indices=desvar_indices,
-                                        flat_indices=True)
-
-                # self.set_input_defaults(name=f'controls:{name}', val=default_val, units=options['units'])
-            else:
-                num_input_nodes = gd.subset_num_nodes['control_input']
-
-                dvname = f'controls:{name}'
-                shape = options['shape']
-                size = np.prod(shape)
-                if options['opt']:
-                    desvar_indices = get_desvar_indices(size, num_input_nodes,
-                                                        options['fix_initial'], options['fix_final'])
-
-                    if len(desvar_indices) > 0:
-                        coerce_desvar_option = CoerceDesvar(num_input_nodes, desvar_indices,
-                                                            options=options)
-
-                        lb = np.zeros_like(desvar_indices, dtype=float)
-                        lb[:] = -INF_BOUND if coerce_desvar_option('lower') is None else \
-                            coerce_desvar_option('lower')
-
-                        ub = np.zeros_like(desvar_indices, dtype=float)
-                        ub[:] = INF_BOUND if coerce_desvar_option('upper') is None else \
-                            coerce_desvar_option('upper')
-
-                        self.add_design_var(name=dvname,
-                                            lower=lb,
-                                            upper=ub,
-                                            scaler=coerce_desvar_option('scaler'),
-                                            adder=coerce_desvar_option('adder'),
-                                            ref0=coerce_desvar_option('ref0'),
-                                            ref=coerce_desvar_option('ref'),
-                                            indices=desvar_indices,
-                                            flat_indices=True)
-
-                # default_val = reshape_val(options['val'], shape, num_input_nodes)
-
-                # self.set_input_defaults(name=dvname, val=default_val, units=options['units'])
-
-    def _configure_desvars(self):
-        control_options = self.options['control_options']
-        gd = self.options['grid_data']
-
-        for name, options in control_options.items():
-            if options['control_type'] == 'polynomial':
-                num_input_nodes = options['order'] + 1
-                shape = options['shape']
-                # default_val = reshape_val(options['val'], shape, num_input_nodes)
                 if options['opt']:
 
                     desvar_indices = np.arange(num_input_nodes, dtype=int)
@@ -333,7 +266,6 @@ class ControlInterpComp(om.ExplicitComponent):
                                         indices=desvar_indices,
                                         flat_indices=True)
 
-                # self.set_input_defaults(name=f'controls:{name}', val=default_val, units=options['units'])
             else:
                 num_input_nodes = gd.subset_num_nodes['control_input']
 
@@ -366,10 +298,6 @@ class ControlInterpComp(om.ExplicitComponent):
                                             indices=desvar_indices,
                                             flat_indices=True)
 
-                # default_val = reshape_val(options['val'], shape, num_input_nodes)
-
-                # self.set_input_defaults(name=dvname, val=default_val, units=options['units'])
-
     def configure_io(self):
         """
         I/O creation is delayed until configure so we can determine shape and units for the states.
diff --git a/dymos/transcriptions/explicit_shooting/explicit_shooting.py b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
index d8c859a6..d9e02abf 100644
--- a/dymos/transcriptions/explicit_shooting/explicit_shooting.py
+++ b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
@@ -9,7 +9,7 @@ from .explicit_shooting_continuity_comp import ExplicitShootingContinuityComp
 from ..transcription_base import TranscriptionBase
 from ..grid_data import BirkhoffGrid, GaussLobattoGrid, RadauGrid, UniformGrid, ChebyshevGaussLobattoGrid
 from .ode_integration_comp import ODEIntegrationComp
-from ...utils.misc import get_rate_units, CoerceDesvar
+from ...utils.misc import get_rate_units, CoerceDesvar, reshape_val
 from ...utils.indexing import get_src_indices_by_row
 from ...utils.introspection import get_promoted_vars, get_source_metadata, get_targets, _get_targets_metadata
 from ...utils.constants import INF_BOUND
@@ -414,6 +414,8 @@ class ExplicitShooting(TranscriptionBase):
                 ncin = self.options['grid'].subset_num_nodes['control_input']
 
             phase.promotes('integrator', inputs=[f'controls:{control_name}'])
+            default_val = reshape_val(options['val'], options['shape'], ncin)
+            phase.set_input_defaults(f'controls:{control_name}', val=default_val)
 
             if options['opt']:
                 coerce_desvar_option = CoerceDesvar(num_input_nodes=ncin, options=options)
@@ -708,6 +710,11 @@ class ExplicitShooting(TranscriptionBase):
         control_rate_continuity : bool
             True if any control rate continuity is required to be enforced.
         """
+        if not self.options['propagate_derivs']:
+            # We're not propagating derivatives because we're just doing a simulation run_model.
+            # No continuity is needed.
+            return False, False, False
+
         num_seg = self.options['grid'].num_segments
         compressed = self.options['grid'].compressed
         transcription = self.options['grid'].transcription
diff --git a/dymos/transcriptions/explicit_shooting/ode_evaluation_group.py b/dymos/transcriptions/explicit_shooting/ode_evaluation_group.py
index acf93d30..f37d1b93 100644
--- a/dymos/transcriptions/explicit_shooting/ode_evaluation_group.py
+++ b/dymos/transcriptions/explicit_shooting/ode_evaluation_group.py
@@ -225,7 +225,9 @@ class ODEEvaluationGroup(om.Group):
 
             self.connect(rate_path, f'state_rate_collector.state_rates_in:{name}_rate')
 
-            self.add_constraint(f'state_rate_collector.state_rates:{name}_rate')
+            if self._compute_derivs:
+                # Adding the constraint/responds lets use compute the derivatives for this.
+                self.add_constraint(f'state_rate_collector.state_rates:{name}_rate')
 
     def _configure_params(self):
         ode_inputs = get_promoted_vars(self.ode, iotypes='input', metadata_keys=['shape', 'units', 'val', 'tags'])
@@ -290,9 +292,12 @@ class ODEEvaluationGroup(om.Group):
 
                 self._ivc.add_output(uhat_name, shape=(num_control_input_nodes,) + shape, units=units)
                 self.add_design_var(uhat_name)
-                self.add_constraint(u_name)
-                self.add_constraint(u_rate_name)
-                self.add_constraint(u_rate2_name)
+
+                if self._compute_derivs:
+                    # Adding the constraint/responds lets use compute the derivatives for this.
+                    self.add_constraint(u_name)
+                    self.add_constraint(u_rate_name)
+                    self.add_constraint(u_rate2_name)
 
                 self.promotes('control_interp', inputs=[uhat_name],
                               outputs=[u_name, u_rate_name, u_rate2_name])
diff --git a/dymos/utils/misc.py b/dymos/utils/misc.py
index c0e11506..33868607 100644
--- a/dymos/utils/misc.py
+++ b/dymos/utils/misc.py
@@ -70,7 +70,7 @@ def reshape_val(val, shape, num_input_nodes):
         The given value of the correct shape.
     """
     if np.isscalar(val) or np.prod(np.asarray(val).shape) == 1:
-        shaped_val = float(val) * np.ones((num_input_nodes,) + shape)
+        shaped_val = float(np.asarray(val).ravel()[0]) * np.ones((num_input_nodes,) + shape)
     elif np.asarray(val).shape == shape:
         shaped_val = np.repeat(val[np.newaxis, ...], num_input_nodes, axis=0)
     else:


Test patch:
diff --git a/dymos/examples/balanced_field/test/test_balanced_field_length.py b/dymos/examples/balanced_field/test/test_balanced_field_length.py
index 16f68e1a..79b0e215 100644
--- a/dymos/examples/balanced_field/test/test_balanced_field_length.py
+++ b/dymos/examples/balanced_field/test/test_balanced_field_length.py
@@ -88,18 +88,148 @@ class TestBalancedFieldLengthRestart(unittest.TestCase):
 @use_tempdirs
 class TestBalancedFieldLengthDefaultValues(unittest.TestCase):
 
-    @require_pyoptsparse(optimizer='IPOPT')
     def test_default_vals_stick(self):
         """
         Make the balanced field problem without any set_val calls after setup.
         """
-        p = make_balanced_field_length_problem(ode_class=BalancedFieldODEComp, tx=dm.GaussLobatto(num_segments=3))
+        ode_class = BalancedFieldODEComp
+        tx = dm.GaussLobatto(num_segments=5, order=3, compressed=True)
+
+        p = om.Problem()
+
+        # First Phase: Brake release to V1 - both engines operable
+        br_to_v1 = dm.Phase(ode_class=ode_class, transcription=tx,
+                            ode_init_kwargs={'mode': 'runway'})
+        br_to_v1.set_time_options(fix_initial=True, duration_bounds=(1, 1000), duration_ref=10.0)
+        br_to_v1.add_state('r', fix_initial=True, lower=0, ref=1000.0, defect_ref=1000.0)
+        br_to_v1.add_state('v', fix_initial=True, lower=0, ref=100.0, defect_ref=100.0)
+        br_to_v1.add_parameter('alpha', val=0.0, opt=False, units='deg')
+        br_to_v1.add_timeseries_output('*')
+
+        # Second Phase: Rejected takeoff at V1 - no engines operable
+        rto = dm.Phase(ode_class=ode_class, transcription=tx,
+                       ode_init_kwargs={'mode': 'runway'})
+        rto.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0)
+        rto.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0)
+        rto.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0)
+        rto.add_parameter('alpha', val=0.0, opt=False, units='deg')
+        rto.add_timeseries_output('*')
+
+        # Third Phase: V1 to Vr - single engine operable
+        v1_to_vr = dm.Phase(ode_class=ode_class, transcription=tx,
+                            ode_init_kwargs={'mode': 'runway'})
+        v1_to_vr.set_time_options(fix_initial=False, duration_bounds=(1, 1000), duration_ref=1.0)
+        v1_to_vr.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0)
+        v1_to_vr.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0)
+        v1_to_vr.add_parameter('alpha', val=0.0, opt=False, units='deg')
+        v1_to_vr.add_timeseries_output('*')
+
+        # Fourth Phase: Rotate - single engine operable
+        rotate = dm.Phase(ode_class=ode_class, transcription=tx,
+                          ode_init_kwargs={'mode': 'runway'})
+        rotate.set_time_options(fix_initial=False, initial_val=35.0, duration_val=5.0,
+                                duration_bounds=(1.0, 5), duration_ref=1.0)
+        rotate.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0)
+        rotate.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0)
+        rotate.add_control('alpha', order=1, opt=True, units='deg', lower=0, upper=10, ref=10,
+                           val=[0, 10], control_type='polynomial')
+        rotate.add_timeseries_output('*')
+
+        # Fifth Phase: Climb to target speed and altitude at end of runway.
+        climb = dm.Phase(ode_class=ode_class, transcription=tx,
+                         ode_init_kwargs={'mode': 'climb'})
+        climb.set_time_options(fix_initial=False, duration_bounds=(1, 100), duration_ref=1.0)
+        climb.add_state('r', fix_initial=False, lower=0, ref=1000.0, defect_ref=1000.0)
+        climb.add_state('h', fix_initial=True, lower=0, ref=1.0, defect_ref=1.0,
+                        val=np.linspace(0.0, 35.0, 6) * 0.3048)
+        climb.add_state('v', fix_initial=False, lower=0, ref=100.0, defect_ref=100.0)
+        climb.add_state('gam', fix_initial=True, lower=0, ref=0.05, defect_ref=0.05,
+                        val=np.radians(np.linspace(0.0, 5.0, 6)))
+        climb.add_control('alpha', opt=True, val=5.0, units='deg', lower=-10, upper=15, ref=10)
+        climb.add_timeseries_output('*')
+
+        # Instantiate the trajectory and add phases
+        traj = dm.Trajectory()
+        p.model.add_subsystem('traj', traj)
+        traj.add_phase('br_to_v1', br_to_v1)
+        traj.add_phase('rto', rto)
+        traj.add_phase('v1_to_vr', v1_to_vr)
+        traj.add_phase('rotate', rotate)
+        traj.add_phase('climb', climb)
+
+        all_phases = ['br_to_v1', 'v1_to_vr', 'rto', 'rotate', 'climb']
+        groundroll_phases = ['br_to_v1', 'v1_to_vr', 'rto', 'rotate']
+
+        # Add parameters common to multiple phases to the trajectory
+        traj.add_parameter('m', val=174200., opt=False, units='lbm',
+                           desc='aircraft mass',
+                           targets={phase: ['m'] for phase in all_phases})
+
+        # Handle parameters which change from phase to phase.
+        traj.add_parameter('T_nominal', val=27000 * 2, opt=False, units='lbf', static_target=True,
+                           desc='nominal aircraft thrust',
+                           targets={'br_to_v1': ['T']})
+
+        traj.add_parameter('T_engine_out', val=27000, opt=False, units='lbf', static_target=True,
+                           desc='thrust under a single engine',
+                           targets={'v1_to_vr': ['T'], 'rotate': ['T'], 'climb': ['T']})
+
+        traj.add_parameter('T_shutdown', val=0.0, opt=False, units='lbf', static_target=True,
+                           desc='thrust when engines are shut down for rejected takeoff',
+                           targets={'rto': ['T']})
+
+        traj.add_parameter('mu_r_nominal', val=0.03, opt=False, units=None, static_target=True,
+                           desc='nominal runway friction coefficient',
+                           targets={'br_to_v1': ['mu_r'], 'v1_to_vr': ['mu_r'],  'rotate': ['mu_r']})
+
+        traj.add_parameter('mu_r_braking', val=0.3, opt=False, units=None, static_target=True,
+                           desc='runway friction coefficient under braking',
+                           targets={'rto': ['mu_r']})
+
+        traj.add_parameter('h_runway', val=0., opt=False, units='ft',
+                           desc='runway altitude',
+                           targets={phase: ['h'] for phase in groundroll_phases})
+
+        # Here we're omitting some constants that are common throughout all phases for the sake of brevity.
+        # Their correct defaults are specified in add_input calls to `wrap_ode_func`.
+
+        # Standard "end of first phase to beginning of second phase" linkages
+        # Alpha changes from being a parameter in v1_to_vr to a polynomial control
+        # in rotate, to a dynamic control in `climb`.
+        traj.link_phases(['br_to_v1', 'v1_to_vr'], vars=['time', 'r', 'v'])
+        traj.link_phases(['v1_to_vr', 'rotate'], vars=['time', 'r', 'v', 'alpha'])
+        traj.link_phases(['rotate', 'climb'], vars=['time', 'r', 'v', 'alpha'])
+        traj.link_phases(['br_to_v1', 'rto'], vars=['time', 'r', 'v'])
+
+        # Less common "final value of r must match at ends of two phases".
+        traj.add_linkage_constraint(phase_a='rto', var_a='r', loc_a='final',
+                                    phase_b='climb', var_b='r', loc_b='final',
+                                    ref=1000)
+
+        # Define the constraints and objective for the optimal control problem
+        v1_to_vr.add_boundary_constraint('v_over_v_stall', loc='final', lower=1.2, ref=100)
+
+        rto.add_boundary_constraint('v', loc='final', equals=0., ref=100, linear=True)
+
+        rotate.add_boundary_constraint('F_r', loc='final', equals=0, ref=100000)
+
+        climb.add_boundary_constraint('h', loc='final', equals=35, ref=35, units='ft', linear=True)
+        climb.add_boundary_constraint('gam', loc='final', equals=5, ref=5, units='deg', linear=True)
+        climb.add_path_constraint('gam', lower=0, upper=5, ref=5, units='deg')
+        climb.add_boundary_constraint('v_over_v_stall', loc='final', lower=1.25, ref=1.25)
+
+        rto.add_objective('r', loc='final', ref=1000.0)
+
+        #
+        # Setup the problem and set the initial guess
+        #
+        p.setup(check=False)
 
         p.run_model()
 
         assert_near_equal(p.get_val('traj.rotate.t_initial'), 35)
         assert_near_equal(p.get_val('traj.rotate.t_duration'), 5)
-        assert_near_equal(p.get_val('traj.rotate.controls:alpha'), np.array([[0, 0]]).T)
+        assert_near_equal(p.get_val('traj.rotate.controls:alpha'), np.array([[0, 10]]).T)
         assert_near_equal(p.get_val('traj.climb.controls:alpha', units='deg'),
                           p.model.traj.phases.climb.interp('', [5, 5], nodes='control_input'))
         assert_near_equal(p.get_val('traj.climb.states:gam', units='deg'),
diff --git a/dymos/transcriptions/common/test/test_control_interp_comp.py b/dymos/transcriptions/common/test/test_control_interp_comp.py
index fe90123e..f4e1ccb9 100644
--- a/dymos/transcriptions/common/test/test_control_interp_comp.py
+++ b/dymos/transcriptions/common/test/test_control_interp_comp.py
@@ -92,8 +92,10 @@ class TestControlRateComp(unittest.TestCase):
 
                 p = om.Problem(model=om.Group())
 
-                controls = {'a': {'units': 'm', 'shape': (1,), 'dynamic': True, 'opt': False},
-                            'b': {'units': 'm', 'shape': (1,), 'dynamic': True, 'opt': False}}
+                controls = {'a': {'units': 'm', 'shape': (1,), 'val': 1.0,
+                                  'dynamic': True, 'opt': False},
+                            'b': {'units': 'm', 'shape': (1,), 'val': 1.0,
+                                  'dynamic': True, 'opt': False}}
 
                 ivc = om.IndepVarComp()
                 p.model.add_subsystem('ivc', ivc, promotes_outputs=['*'])
@@ -183,7 +185,8 @@ class TestControlRateComp(unittest.TestCase):
 
                 p = om.Problem(model=om.Group())
 
-                controls = {'a': {'units': 'm', 'shape': (3,), 'dynamic': True, 'opt': False}}
+                controls = {'a': {'units': 'm', 'shape': (3,), 'val': np.array([1, 2, 3]),
+                                  'dynamic': True, 'opt': False}}
 
                 ivc = om.IndepVarComp()
                 p.model.add_subsystem('ivc', ivc, promotes_outputs=['*'])
@@ -283,7 +286,8 @@ class TestControlRateComp(unittest.TestCase):
 
                 p = om.Problem(model=om.Group())
 
-                controls = {'a': {'units': 'm', 'shape': (3, 1), 'dynamic': True, 'opt': False}}
+                controls = {'a': {'units': 'm', 'shape': (3, 1), 'val': 1.0,
+                                  'dynamic': True, 'opt': False}}
 
                 ivc = om.IndepVarComp()
                 p.model.add_subsystem('ivc', ivc, promotes_outputs=['*'])
@@ -382,7 +386,8 @@ class TestControlRateComp(unittest.TestCase):
 
                 p = om.Problem(model=om.Group())
 
-                controls = {'a': {'units': 'm', 'shape': (2, 2), 'dynamic': True, 'opt': False}}
+                controls = {'a': {'units': 'm', 'shape': (2, 2), 'val': np.array([[1, 2], [3, 4]]),
+                                  'dynamic': True, 'opt': False}}
 
                 ivc = om.IndepVarComp()
                 p.model.add_subsystem('ivc', ivc, promotes_outputs=['*'])