Repository: OpenMDAO/dymos

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_

Hints from issue comments:

Solution patch:
diff --git a/dymos/transcriptions/common/control_comp.py b/dymos/transcriptions/common/control_comp.py
index 74a1449b4..ca4a18c08 100644
--- a/dymos/transcriptions/common/control_comp.py
+++ b/dymos/transcriptions/common/control_comp.py
@@ -5,7 +5,7 @@
 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 @@ def _configure_controls(self):
             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 @@ def _configure_controls(self):
                 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 @@ def _configure_controls(self):
                 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 @@ def _configure_desvars(self):
             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 @@ def _configure_desvars(self):
                                         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 @@ def _configure_desvars(self):
                                             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 d8c859a65..d9e02abfb 100644
--- a/dymos/transcriptions/explicit_shooting/explicit_shooting.py
+++ b/dymos/transcriptions/explicit_shooting/explicit_shooting.py
@@ -9,7 +9,7 @@
 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 @@ def configure_controls(self, phase):
                 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 @@ def _requires_continuity_constraints(self, phase):
         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 acf93d304..f37d1b93a 100644
--- a/dymos/transcriptions/explicit_shooting/ode_evaluation_group.py
+++ b/dymos/transcriptions/explicit_shooting/ode_evaluation_group.py
@@ -225,7 +225,9 @@ def _configure_states(self):
 
             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 @@ def _configure_controls(self):
 
                 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 c0e115063..338686076 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 16f68e1a3..79b0e2151 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 @@ def test_restart_from_sim(self):
 @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 fe90123e5..f4e1ccb9c 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 @@ def test_control_interp_scalar(self):
 
                 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 @@ def test_control_interp_vector(self, transcription='gauss-lobatto', compressed=T
 
                 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 @@ def test_control_interp_matrix_3x1(self, transcription='gauss-lobatto', compress
 
                 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 @@ def test_control_interp_matrix_2x2(self, transcription='gauss-lobatto', compress
 
                 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=['*'])