Repository: OpenMDAO/dymos

Issue:
Trajectory `link_phases` method should support arguments of `add_linkage_constraint`.
### Proposed feature.

The `link_phases` method is a simpler wrapper for the more complex but more capable `add_linkage_constraint` method.
In order to allow the user to control things like the scaling factors for the imposed constraints, `link_phases` should support the following options and pass them to `add_linkage_constraint`.

- `units`
-  `scaler`
- `adder`
-  `ref0`
-  `ref`
-  `linear`

If `connected` is specified, then the above options will have no effect.

Note `sign_a` and `sign_b` are explicitly left out here, since the intent of `link_phases` is to provide continuity in value.





### Example

```python)

traj.link_phases(
    phases=['ascent', 'descent'],
    vars=['time', 'r', 'h', 'gam'],
    scaler=1.0E-3)
```

Copy/paste error in `Trajectory._configure_linkages`
### Description

When connecting a parameter between two phases via linkage, the `_configure_linkages` method contains a copy/paste error when checking to see if the second variable is fixed.

https://github.com/OpenMDAO/dymos/blob/523217e13bb294694b4bd27a7960ad41cc70b1c4/dymos/trajectory/trajectory.py#L776

### Example

```python
import openmdao.api as om
import dymos as dm
from dymos.examples.finite_burn_orbit_raise.finite_burn_eom import FiniteBurnODE

p = om.Problem()

traj = dm.Trajectory()
p.model.add_subsystem('traj', subsys=traj)

# First Phase (burn)
burn1 = dm.Phase(ode_class=FiniteBurnODE, transcription=dm.GaussLobatto(num_segments=4, order=3))

traj.add_phase('burn1', burn1)

burn1.set_time_options(fix_initial=True, duration_bounds=(.5, 10), units='TU')
burn1.add_state('r', fix_initial=True, fix_final=False,
                rate_source='r_dot', units='DU')
burn1.add_state('theta', fix_initial=True, fix_final=False,
                rate_source='theta_dot', units='rad')
burn1.add_state('vr', fix_initial=True, fix_final=False,
                rate_source='vr_dot', units='DU/TU')
burn1.add_state('vt', fix_initial=True, fix_final=False,
                rate_source='vt_dot', units='DU/TU')
burn1.add_state('accel', fix_initial=True, fix_final=False,
                rate_source='at_dot', units='DU/TU**2')
burn1.add_state('deltav', fix_initial=True, fix_final=False,
                rate_source='deltav_dot', units='DU/TU')
burn1.add_control('u1',  rate_continuity=True, rate2_continuity=True,
                  units='deg', scaler=0.01, lower=-30, upper=30)
burn1.add_parameter('c', opt=False, val=1.5, units='DU/TU')

# Second Phase (Coast)

coast = dm.Phase(ode_class=FiniteBurnODE, transcription=dm.GaussLobatto(num_segments=10, order=3))

traj.add_phase('coast', coast)

coast.set_time_options(initial_bounds=(0.5, 20), duration_bounds=(.5, 10), duration_ref=10, units='TU')
coast.add_state('r', fix_initial=False, fix_final=False, defect_scaler=100.0,
                rate_source='r_dot', units='DU')
coast.add_state('theta', fix_initial=False, fix_final=False, defect_scaler=100.0,
                rate_source='theta_dot', units='rad')
coast.add_state('vr', fix_initial=False, fix_final=False, defect_scaler=100.0,
                rate_source='vr_dot', units='DU/TU')
coast.add_state('vt', fix_initial=False, fix_final=False, defect_scaler=100.0,
                rate_source='vt_dot', units='DU/TU')
coast.add_state('accel', fix_initial=True, fix_final=True, ref=1.0E-12, defect_ref=1.0E-12,
                rate_source='at_dot', units='DU/TU**2')
coast.add_state('deltav', fix_initial=False, fix_final=False,
                rate_source='deltav_dot', units='DU/TU')
coast.add_control('u1', opt=False, val=0.0, units='deg')
coast.add_parameter('c', opt=False, val=1.5, units='DU/TU')

# Third Phase (burn)

burn2 = dm.Phase(ode_class=FiniteBurnODE, transcription=dm.GaussLobatto(num_segments=3, order=3))

traj.add_phase('burn2', burn2)

burn2.set_time_options(initial_bounds=(0.5, 20), duration_bounds=(.5, 10), initial_ref=10, units='TU')
burn2.add_state('r', fix_initial=False, fix_final=True,
                rate_source='r_dot', units='DU')
burn2.add_state('theta', fix_initial=False, fix_final=False,
                rate_source='theta_dot', units='rad')
burn2.add_state('vr', fix_initial=False, fix_final=True,
                rate_source='vr_dot', units='DU/TU')
burn2.add_state('vt', fix_initial=False, fix_final=True,
                rate_source='vt_dot', units='DU/TU')
burn2.add_state('accel', fix_initial=False, fix_final=False, defect_ref=1.0E-6,
                rate_source='at_dot', units='DU/TU**2')
burn2.add_state('deltav', fix_initial=False, fix_final=False,
                rate_source='deltav_dot', units='DU/TU')
burn2.add_control('u1', rate_continuity=True, rate2_continuity=True,
                  units='deg', scaler=0.01, lower=-30, upper=30)
burn2.add_parameter('c', opt=False, val=1.5, units='DU/TU')

burn2.add_objective('deltav', loc='final')

# Link Phases
traj.link_phases(phases=['burn1', 'coast', 'burn2'],
                 vars=['time', 'r', 'theta', 'vr', 'vt', 'deltav', 'c'])

with self.assertRaises(ValueError) as e:
    p.setup()
```


### Dymos Version

1.6.1-dev

### Relevant environment information

_No response_

Discussion:

Patch:
diff --git a/dymos/phase/phase.py b/dymos/phase/phase.py
index 023549b3..69dac6ce 100644
--- a/dymos/phase/phase.py
+++ b/dymos/phase/phase.py
@@ -2316,15 +2316,39 @@ class Phase(om.Group):
         bool
             True if the state of the given name is guaranteed to be fixed at the given location.
         """
-        if loc == 'initial':
-            res = self.control_options[name]['fix_initial']
-        elif loc == 'final':
-            res = self.control_options[name]['fix_final']
+        control_opts = self.control_options[name]
+        if loc == 'initial' and control_opts['fix_initial']:
+            res = True
+        elif loc == 'final' and control_opts['fix_final']:
+            res = True
         else:
-            raise ValueError(f'Unknown value for argument "loc": must be either "initial" or '
-                             f'"final" but got {loc}')
+            res = not control_opts['opt']
         return res
 
+    def is_control_rate_fixed(self, name, loc):
+        """
+        Test if the control rate of the given name is guaranteed to be fixed at the initial or final time.
+
+        Parameters
+        ----------
+        name : str
+            The name of the control to be tested.
+        loc : str
+            The location of time to be tested: either 'initial' or 'final'.
+
+        Returns
+        -------
+        bool
+            True if the state of the given name is guaranteed to be fixed at the given location.
+        """
+        if name.endswith('_rate') and self.control_options is not None and \
+                name[:-5] in self.control_options:
+            control_name = name[:-5]
+        elif name.endswith('_rate2') and self.control_options is not None and \
+                name[:-6] in self.control_options:
+            control_name = name[:-6]
+        return self.is_control_fixed(control_name, loc)
+
     def is_polynomial_control_fixed(self, name, loc):
         """
         Test if the polynomial control of the given name is guaranteed to be fixed at the initial or final time.
@@ -2350,6 +2374,30 @@ class Phase(om.Group):
                              f'"final" but got {loc}')
         return res
 
+    def is_polynomial_control_rate_fixed(self, name, loc):
+        """
+        Test if the polynomial control rate of the given name is guaranteed to be fixed at the initial or final time.
+
+        Parameters
+        ----------
+        name : str
+            The name of the control to be tested.
+        loc : str
+            The location of time to be tested: either 'initial' or 'final'.
+
+        Returns
+        -------
+        bool
+            True if the state of the given name is guaranteed to be fixed at the given location.
+        """
+        if name.endswith('_rate') and self.polynomial_control_options is not None and \
+                name[:-5] in self.polynomial_control_options:
+            control_name = name[:-5]
+        elif name.endswith('_rate2') and self.polynomial_control_options is not None and \
+                name[:-6] in self.options['polynomial_control_options']:
+            control_name = name[:-6]
+        return self.is_polynomial_control_fixed(control_name, loc)
+
     def _indices_in_constraints(self, name, loc):
         """
         Returns a set of the C-order flattened indices involving constraint of the given name at the given loc.
diff --git a/dymos/trajectory/trajectory.py b/dymos/trajectory/trajectory.py
index 73ce8335..2ce571ff 100644
--- a/dymos/trajectory/trajectory.py
+++ b/dymos/trajectory/trajectory.py
@@ -632,7 +632,7 @@ class Trajectory(om.Group):
                                                     sign_b=options['sign_b'])
                     self._linkages[phase_pair].pop(var_pair)
 
-    def _is_valid_linkage(self, phase_name_a, phase_name_b, loc_a, loc_b, var_a, var_b):
+    def _is_valid_linkage(self, phase_name_a, phase_name_b, loc_a, loc_b, var_a, var_b, fixed_a, fixed_b):
         """
         Validates linkage constraints.
 
@@ -655,6 +655,10 @@ class Trajectory(om.Group):
             The variable name of the first side of the linkage.
         var_b : str
             The variable name of the second side of the linkage.
+        fixed_a : bool
+            True if variable a is fixed at the linkage location.
+        fixed_b : bool
+            True if variable b is fixed at the linkage location
 
         Returns
         -------
@@ -759,10 +763,12 @@ class Trajectory(om.Group):
                     fixed_a = phase_a.is_control_fixed(var_a, loc_a)
                 elif class_a in {'input_polynomial_control', 'indep_polynomial_control'}:
                     fixed_a = phase_a.is_polynomial_control_fixed(var_a, loc_a)
+                elif class_a in {'control_rate', 'control_rate2'}:
+                    fixed_a = phase_a.is_control_rate_fixed(var_a, loc_a)
                 elif class_a == 'parameter':
                     fixed_a = not phase_a.parameter_options[var_a]['opt']
                 else:
-                    fixed_a = True
+                    fixed_a = False  # No way to know so we allow these to go through
 
                 if class_b == 'time':
                     fixed_b = phase_b.is_time_fixed(loc_b)
@@ -772,10 +778,12 @@ class Trajectory(om.Group):
                     fixed_b = phase_b.is_control_fixed(var_b, loc_b)
                 elif class_b in {'input_polynomial_control', 'indep_polynomial_control'}:
                     fixed_b = phase_b.is_polynomial_control_fixed(var_b, loc_b)
-                elif class_a == 'parameter':
-                    fixed_a = not phase_b.parameter_options[var_a]['opt']
+                elif class_b in {'control_rate', 'control_rate2'}:
+                    fixed_b = phase_b.is_control_rate_fixed(var_b, loc_b)
+                elif class_b == 'parameter':
+                    fixed_b = not phase_b.parameter_options[var_b]['opt']
                 else:
-                    fixed_b = True
+                    fixed_b = False  # No way to know so we allow these to go through
 
                 prefixed_a = _get_prefixed_var(var_a, phase_a)
                 prefixed_b = _get_prefixed_var(var_b, phase_b)
@@ -808,10 +816,13 @@ class Trajectory(om.Group):
                     _print_on_rank(f'{indent * 2}{prefixed_a:<{padding_a}s} [{loc_a}{str_fixed_a}] ->  '
                                    f'{prefixed_b:<{padding_b}s} [{loc_b}{str_fixed_b}]')
                 else:
-                    is_valid, msg = self._is_valid_linkage(phase_name_a, phase_name_b,
-                                                           loc_a, loc_b, var_a, var_b)
 
-                    if not is_valid:
+                    if fixed_a and fixed_b:
+                        msg = f'Cannot link {loc_a} value of "{var_a}" in {phase_name_a} to {loc_b} ' \
+                              f'value of "{var_b}" in {phase_name_b}. Values on both sides of the linkage are fixed ' \
+                              'and the linkage is enforced via constraint. Either link the variables via connection ' \
+                              'or make the variables design variables on at least one side of the connection.'
+
                         raise ValueError(f'Invalid linkage in Trajectory {self.pathname}: {msg}')
 
                     linkage_comp.add_linkage_configure(options)
@@ -924,8 +935,8 @@ class Trajectory(om.Group):
         sign_b : float
             The sign applied to the variable from the second phase in the linkage constraint.
         units : str or None or _unspecified
-            Units of the linkage.  If _unspecified, dymos will use the units from the variable
-            in the first phase of the linkage.  Units of the two specified variables must be
+            Units of the linkage. If _unspecified, dymos will use the units from the variable
+            in the first phase of the linkage. Units of the two specified variables must be
             compatible.
         lower : float or array or None
             The lower bound applied as a constraint on the linkage equation.
@@ -942,16 +953,16 @@ class Trajectory(om.Group):
         ref : float or array or None
             The unit-reference value of the linkage constraint.
         linear : bool
-            If True, treat this variable as a linear constraint, otherwise False.  Linear
+            If True, treat this variable as a linear constraint, otherwise False. Linear
             constraints should only be applied if the variable on each end of the linkage is a
             design variable or a linear function of one.
         connected : bool
             If True, this constraint is enforced by direct connection rather than a constraint
-            for the optimizer.  This is only valid for states and time.
+            for the optimizer. This is only valid for states and time.
         """
         if connected:
             invalid_options = []
-            for arg in ['lower', 'upper', 'equals', 'scaler', 'adder', 'ref0', 'ref']:
+            for arg in ['lower', 'upper', 'equals', 'scaler', 'adder', 'ref0', 'ref', 'units']:
                 if locals()[arg] is not None:
                     invalid_options.append(arg)
             if locals()['linear']:
@@ -961,7 +972,7 @@ class Trajectory(om.Group):
                       f'in trajectory {self.pathname}. The following options for ' \
                       f'add_linkage_constraint were specified but not valid when ' \
                       f'option \'connected\' is True: ' + ' '.join(invalid_options)
-                warnings.warn(msg)
+                om.issue_warning(msg, category=om.UnusedOptionWarning)
 
         if (phase_a, phase_b) not in self._linkages:
             self._linkages[phase_a, phase_b] = OrderedDict()
@@ -986,7 +997,8 @@ class Trajectory(om.Group):
         d['linear'] = linear
         d['connected'] = connected
 
-    def link_phases(self, phases, vars=None, locs=('final', 'initial'), connected=False):
+    def link_phases(self, phases, vars=None, locs=('final', 'initial'), connected=False,
+                    units=_unspecified, scaler=None, adder=None, ref0=None, ref=None, linear=False):
         """
         Specify that phases in the given sequence are to be assume continuity of the given variables.
 
@@ -1015,6 +1027,22 @@ class Trajectory(om.Group):
         connected : bool
             Set to True to directly connect the phases being linked. Otherwise, create constraints
             for the optimizer to solve.
+        units : str or None or _unspecified
+            Units of the linkage.  If _unspecified, dymos will use the units from the variable
+            in the first phase of the linkage.  Units of the two specified variables must be
+            compatible.
+        scaler : float or array or None
+            The scaler of the linkage constraint.
+        adder : float or array or None
+            The adder of the linkage constraint.
+        ref0 : float or array or None
+            The zero-reference value of the linkage constraint.
+        ref : float or array or None
+            The unit-reference value of the linkage constraint.
+        linear : bool
+            If True, treat this variable as a linear constraint, otherwise False.  Linear
+            constraints should only be applied if the variable on each end of the linkage is a
+            design variable or a linear function of one.
 
         See Also
         --------
@@ -1050,7 +1078,8 @@ class Trajectory(om.Group):
             for var in _vars:
                 self.add_linkage_constraint(phase_a=phase_name_a, phase_b=phase_name_b,
                                             var_a=var, var_b=var, loc_a=loc_a, loc_b=loc_b,
-                                            connected=connected)
+                                            connected=connected, units=units,
+                                            scaler=scaler, adder=adder, ref0=ref0, ref=ref, linear=linear)
 
     def _constraint_report(self, outstream=sys.stdout):
         if self.options['sim_mode']:
diff --git a/dymos/visualization/linkage/report.py b/dymos/visualization/linkage/report.py
index 7c298efa..0e943798 100644
--- a/dymos/visualization/linkage/report.py
+++ b/dymos/visualization/linkage/report.py
@@ -197,8 +197,15 @@ def _linkages_to_list(traj, model_data):
             loc_a = options['loc_a']
             loc_b = options['loc_b']
 
-            tree_var_a = tree[CBN][phase_name_a][CBN][loc_a][CBN][var_a]
-            tree_var_b = tree[CBN][phase_name_b][CBN][loc_b][CBN][var_b]
+            try:
+                tree_var_a = tree[CBN][phase_name_a][CBN][loc_a][CBN][var_a]
+                tree_var_b = tree[CBN][phase_name_b][CBN][loc_b][CBN][var_b]
+            except KeyError:
+                # When linking two ODE outputs, it's possible for the variable to not show up in the tree.
+                # For example, see the 'ke' linkage in
+                # TestTwoPhaseCannonballODEOutputLinkage.test_traj_param_target_unspecified_units.
+                # For now, let these linkages go undisplayed in the report.
+                continue
 
             tree_var_a['linked'] = tree_var_b['linked'] = True
             tree_var_a['connected'] = tree_var_b['connected'] = options['connected']
diff --git a/readme.md b/readme.md
index 299212a4..955fab57 100644
--- a/readme.md
+++ b/readme.md
@@ -106,7 +106,7 @@ import openmdao.api as om
 import dymos as dm
 import matplotlib.pyplot as plt
 
-
+# First define a system which computes the equations of motion
 class BrachistochroneEOM(om.ExplicitComponent):
     def initialize(self):
         self.options.declare('num_nodes', types=int)
@@ -115,68 +115,37 @@ class BrachistochroneEOM(om.ExplicitComponent):
         nn = self.options['num_nodes']
 
         # Inputs
-        self.add_input('v',
-                       shape=(nn,),
-                       desc='velocity',
-                       units='m/s')
-
-        self.add_input('theta',
-                       shape=(nn,),
-                       desc='angle of wire',
-                       units='rad')
-
-        self.add_output('xdot',
-                        shape=(nn,),
-                        desc='velocity component in x',
-                        units='m/s')
-
-        self.add_output('ydot',
-                        shape=(nn,),
-                        desc='velocity component in y',
-                        units='m/s')
-
-        self.add_output('vdot',
-                        val=np.zeros(nn),
-                        desc='acceleration magnitude',
-                        units='m/s**2')
-
-        self.add_output('check',
-                        val=np.zeros(nn),
-                        desc='A check on the solution: v/sin(theta) = constant',
-                        units='m/s')
-
-        # Use OpenMDAO's ability to automatically determine a sparse "coloring" of the jacobian
-        # for this ODE component.
-        self.declare_coloring(wrt='*', method='cs')
-
-    def compute(self, inputs, outputs):
-        theta = inputs['theta']
-        cos_theta = np.cos(theta)
-        sin_theta = np.sin(theta)
-        v = inputs['v']
+        self.add_input('v', val=np.zeros(nn), units='m/s', desc='velocity')
+        self.add_input('theta', val=np.zeros(nn), units='rad', desc='angle of wire')
+        self.add_output('xdot', val=np.zeros(nn), units='m/s', desc='x rate of change')
+        self.add_output('ydot', val=np.zeros(nn), units='m/s', desc='y rate of change')
+        self.add_output('vdot', val=np.zeros(nn), units='m/s**2', desc='v rate of change')
 
-        outputs['vdot'] = 9.80665 * cos_theta
-        outputs['xdot'] = v * sin_theta
-        outputs['ydot'] = -v * cos_theta
-        outputs['check'] = v / sin_theta
+        # Ask OpenMDAO to compute the partial derivatives using complex-step
+        # with a partial coloring algorithm for improved performance
+        self.declare_partials(of='*', wrt='*', method='cs')
+        self.declare_coloring(wrt='*', method='cs', show_summary=True)
 
+    def compute(self, inputs, outputs):
+        v, theta = inputs.values()
+        outputs['vdot'] = 9.80665 * np.cos(theta)
+        outputs['xdot'] = v * np.sin(theta)
+        outputs['ydot'] = -v * np.cos(theta)
 
-# Define the OpenMDAO problem
-p = om.Problem(model=om.Group())
+p = om.Problem()
 
 # Define a Trajectory object
-traj = dm.Trajectory()
-p.model.add_subsystem('traj', subsys=traj)
+traj = p.model.add_subsystem('traj', dm.Trajectory())
 
 # Define a Dymos Phase object with GaussLobatto Transcription
-phase = dm.Phase(ode_class=BrachistochroneEOM,
-                 transcription=dm.Radau(num_segments=10, order=3))
+tx = dm.GaussLobatto(num_segments=10, order=3)
+phase = dm.Phase(ode_class=BrachistochroneEOM, transcription=tx)
+
 traj.add_phase(name='phase0', phase=phase)
 
 # Set the time options
 phase.set_time_options(fix_initial=True,
                        duration_bounds=(0.5, 10.0))
-
 # Set the state options
 phase.add_state('x', rate_source='xdot',
                 fix_initial=True, fix_final=True)
@@ -184,20 +153,19 @@ phase.add_state('y', rate_source='ydot',
                 fix_initial=True, fix_final=True)
 phase.add_state('v', rate_source='vdot',
                 fix_initial=True, fix_final=False)
-
 # Define theta as a control.
 phase.add_control(name='theta', units='rad',
                   lower=0, upper=np.pi)
-
 # Minimize final time.
 phase.add_objective('time', loc='final')
 
 # Set the driver.
 p.driver = om.ScipyOptimizeDriver()
 
-# Allow OpenMDAO to automatically determine our sparsity pattern
-# for TOTAL derivatives. Doing so can significantly speed up the
-# execution of dymos.
+# Allow OpenMDAO to automatically determine total
+# derivative sparsity pattern.
+# This works in conjunction with partial derivative
+# coloring to give a large speedup
 p.driver.declare_coloring()
 
 # Setup the problem
@@ -210,24 +178,20 @@ p.set_val('traj.phase0.t_duration', 2.0)
 # States and controls here use a linearly interpolated
 # initial guess along the trajectory.
 p.set_val('traj.phase0.states:x',
-          phase.interpolate(ys=[0, 10], nodes='state_input'),
+          phase.interp('x', ys=[0, 10]),
           units='m')
-
 p.set_val('traj.phase0.states:y',
-          phase.interpolate(ys=[10, 5], nodes='state_input'),
+          phase.interp('y', ys=[10, 5]),
           units='m')
-
 p.set_val('traj.phase0.states:v',
-          phase.interpolate(ys=[0, 5], nodes='state_input'),
+          phase.interp('v', ys=[0, 5]),
           units='m/s')
+# constant initial guess for control
+p.set_val('traj.phase0.controls:theta', 90, units='deg')
 
-p.set_val('traj.phase0.controls:theta',
-          phase.interpolate(ys=[5, 45], nodes='control_input'),
-          units='deg')
-
-# Use Dymos' run_problem method to run the driver, simulate the results,
-# and record the results to 'dymos_solution.db' and 'dymos_simulation.db'.
-dm.run_problem(p, simulate=True)
+# Run the driver to solve the problem and generate default plots of
+# state and control values vs time
+dm.run_problem(p, make_plots=True, simulate=True)
 
 # Load the solution and simulation files.
 sol_case = om.CaseReader('dymos_solution.db').get_case('final')


Test patch:
diff --git a/dymos/examples/cannonball/test/test_phase_linkage_compliance.py b/dymos/examples/cannonball/test/test_phase_linkage_compliance.py
index d9656f5a..008a25dd 100644
--- a/dymos/examples/cannonball/test/test_phase_linkage_compliance.py
+++ b/dymos/examples/cannonball/test/test_phase_linkage_compliance.py
@@ -104,7 +104,10 @@ class TestTwoPhaseCannonballODEOutputLinkage(unittest.TestCase):
 
         self.assertEqual(str(e.exception), 'Invalid linkage in Trajectory traj: Cannot link final '
                                            'value of "gam" in ascent to initial value of "gam" in '
-                                           'descent.  Values on both sides of the linkage are fixed.')
+                                           'descent. Values on both sides of the '
+                                           'linkage are fixed and the linkage is enforced via constraint. Either link '
+                                           'the variables via connection or make the variables design variables on at '
+                                           'least one side of the connection.')
 
     @require_pyoptsparse(optimizer='SLSQP')
     def test_link_fixed_states_final_to_final(self):
@@ -205,7 +208,10 @@ class TestTwoPhaseCannonballODEOutputLinkage(unittest.TestCase):
 
         self.assertEqual(str(e.exception), 'Invalid linkage in Trajectory traj: Cannot link final '
                                            'value of "h" in ascent to final value of "h" in '
-                                           'descent.  Values on both sides of the linkage are fixed.')
+                                           'descent. Values on both sides of the '
+                                           'linkage are fixed and the linkage is enforced via constraint. Either link '
+                                           'the variables via connection or make the variables design variables on at '
+                                           'least one side of the connection.')
 
     @require_pyoptsparse(optimizer='SLSQP')
     def test_link_fixed_states_initial_to_initial(self):
@@ -306,7 +312,10 @@ class TestTwoPhaseCannonballODEOutputLinkage(unittest.TestCase):
 
         self.assertEqual(str(e.exception), 'Invalid linkage in Trajectory traj: Cannot link initial '
                                            'value of "h" in ascent to initial value of "h" in '
-                                           'descent.  Values on both sides of the linkage are fixed.')
+                                           'descent. Values on both sides of the '
+                                           'linkage are fixed and the linkage is enforced via constraint. Either link '
+                                           'the variables via connection or make the variables design variables on at '
+                                           'least one side of the connection.')
 
     @require_pyoptsparse(optimizer='SLSQP')
     def test_link_fixed_times_final_to_initial(self):
@@ -403,12 +412,14 @@ class TestTwoPhaseCannonballODEOutputLinkage(unittest.TestCase):
 
         self.assertEqual(str(e.exception), 'Invalid linkage in Trajectory traj: Cannot link final '
                                            'value of "time" in ascent to initial value of "time" in '
-                                           'descent.  Values on both sides of the linkage are fixed.')
+                                           'descent. Values on both sides of the '
+                                           'linkage are fixed and the linkage is enforced via constraint. Either link '
+                                           'the variables via connection or make the variables design variables on at '
+                                           'least one side of the connection.')
 
     @require_pyoptsparse(optimizer='SLSQP')
     def test_link_bounded_times_final_to_initial(self):
         """ Test that linking phases with times that are fixed at the linkage point raises an exception. """
-
         import openmdao.api as om
 
         import dymos as dm
@@ -501,7 +512,105 @@ class TestTwoPhaseCannonballODEOutputLinkage(unittest.TestCase):
 
         self.assertEqual(str(e.exception), 'Invalid linkage in Trajectory traj: Cannot link final '
                                            'value of "time" in ascent to initial value of "time" in '
-                                           'descent.  Values on both sides of the linkage are fixed.')
+                                           'descent. Values on both sides of the '
+                                           'linkage are fixed and the linkage is enforced via constraint. Either link '
+                                           'the variables via connection or make the variables design variables on at '
+                                           'least one side of the connection.')
+
+    def test_link_parameters_fixed_to_fixed(self):
+        """ Linking non-opt parameters across phases by connection should work.
+         But linking them via constraint should not (no freedom to move them to satisfy the constraint).
+         """
+        import openmdao.api as om
+        import dymos as dm
+        from dymos.examples.finite_burn_orbit_raise.finite_burn_eom import FiniteBurnODE
+
+        p = om.Problem()
+
+        traj = dm.Trajectory()
+        p.model.add_subsystem('traj', subsys=traj)
+
+        # First Phase (burn)
+        burn1 = dm.Phase(ode_class=FiniteBurnODE, transcription=dm.GaussLobatto(num_segments=4, order=3))
+
+        traj.add_phase('burn1', burn1)
+
+        burn1.set_time_options(fix_initial=True, duration_bounds=(.5, 10), units='TU')
+        burn1.add_state('r', fix_initial=True, fix_final=False,
+                        rate_source='r_dot', units='DU')
+        burn1.add_state('theta', fix_initial=True, fix_final=False,
+                        rate_source='theta_dot', units='rad')
+        burn1.add_state('vr', fix_initial=True, fix_final=False,
+                        rate_source='vr_dot', units='DU/TU')
+        burn1.add_state('vt', fix_initial=True, fix_final=False,
+                        rate_source='vt_dot', units='DU/TU')
+        burn1.add_state('accel', fix_initial=True, fix_final=False,
+                        rate_source='at_dot', units='DU/TU**2')
+        burn1.add_state('deltav', fix_initial=True, fix_final=False,
+                        rate_source='deltav_dot', units='DU/TU')
+        burn1.add_control('u1',  rate_continuity=True, rate2_continuity=True,
+                          units='deg', scaler=0.01, lower=-30, upper=30)
+        burn1.add_parameter('c', opt=False, val=1.5, units='DU/TU')
+
+        # Second Phase (Coast)
+
+        coast = dm.Phase(ode_class=FiniteBurnODE, transcription=dm.GaussLobatto(num_segments=10, order=3))
+
+        traj.add_phase('coast', coast)
+
+        coast.set_time_options(initial_bounds=(0.5, 20), duration_bounds=(.5, 10), duration_ref=10, units='TU')
+        coast.add_state('r', fix_initial=False, fix_final=False, defect_scaler=100.0,
+                        rate_source='r_dot', units='DU')
+        coast.add_state('theta', fix_initial=False, fix_final=False, defect_scaler=100.0,
+                        rate_source='theta_dot', units='rad')
+        coast.add_state('vr', fix_initial=False, fix_final=False, defect_scaler=100.0,
+                        rate_source='vr_dot', units='DU/TU')
+        coast.add_state('vt', fix_initial=False, fix_final=False, defect_scaler=100.0,
+                        rate_source='vt_dot', units='DU/TU')
+        coast.add_state('accel', fix_initial=True, fix_final=True, ref=1.0E-12, defect_ref=1.0E-12,
+                        rate_source='at_dot', units='DU/TU**2')
+        coast.add_state('deltav', fix_initial=False, fix_final=False,
+                        rate_source='deltav_dot', units='DU/TU')
+        coast.add_control('u1', opt=False, val=0.0, units='deg')
+        coast.add_parameter('c', opt=False, val=1.5, units='DU/TU')
+
+        # Third Phase (burn)
+
+        burn2 = dm.Phase(ode_class=FiniteBurnODE, transcription=dm.GaussLobatto(num_segments=3, order=3))
+
+        traj.add_phase('burn2', burn2)
+
+        burn2.set_time_options(initial_bounds=(0.5, 20), duration_bounds=(.5, 10), initial_ref=10, units='TU')
+        burn2.add_state('r', fix_initial=False, fix_final=True,
+                        rate_source='r_dot', units='DU')
+        burn2.add_state('theta', fix_initial=False, fix_final=False,
+                        rate_source='theta_dot', units='rad')
+        burn2.add_state('vr', fix_initial=False, fix_final=True,
+                        rate_source='vr_dot', units='DU/TU')
+        burn2.add_state('vt', fix_initial=False, fix_final=True,
+                        rate_source='vt_dot', units='DU/TU')
+        burn2.add_state('accel', fix_initial=False, fix_final=False, defect_ref=1.0E-6,
+                        rate_source='at_dot', units='DU/TU**2')
+        burn2.add_state('deltav', fix_initial=False, fix_final=False,
+                        rate_source='deltav_dot', units='DU/TU')
+        burn2.add_control('u1', rate_continuity=True, rate2_continuity=True,
+                          units='deg', scaler=0.01, lower=-30, upper=30)
+        burn2.add_parameter('c', opt=False, val=1.5, units='DU/TU')
+
+        burn2.add_objective('deltav', loc='final')
+
+        # Link Phases
+        traj.link_phases(phases=['burn1', 'coast', 'burn2'],
+                         vars=['time', 'r', 'theta', 'vr', 'vt', 'deltav', 'c'])
+
+        with self.assertRaises(ValueError) as e:
+            p.setup()
+
+        self.assertEqual(str(e.exception), 'Invalid linkage in Trajectory traj: Cannot link final value of "c" in '
+                                           'burn1 to initial value of "c" in coast. Values on both sides of the '
+                                           'linkage are fixed and the linkage is enforced via constraint. Either link '
+                                           'the variables via connection or make the variables design variables on at '
+                                           'least one side of the connection.')
 
 
 if __name__ == '__main__':  # pragma: no cover
diff --git a/dymos/examples/finite_burn_orbit_raise/test/test_two_burn_orbit_raise_linkages.py b/dymos/examples/finite_burn_orbit_raise/test/test_two_burn_orbit_raise_linkages.py
index 1a51ea97..754ce937 100644
--- a/dymos/examples/finite_burn_orbit_raise/test/test_two_burn_orbit_raise_linkages.py
+++ b/dymos/examples/finite_burn_orbit_raise/test/test_two_burn_orbit_raise_linkages.py
@@ -171,7 +171,7 @@ class TestTwoBurnOrbitRaiseLinkages(unittest.TestCase):
                         rate_source='at_dot', targets=['accel'], units='DU/TU**2')
         burn1.add_state('deltav', fix_initial=True, fix_final=False,
                         rate_source='deltav_dot', units='DU/TU')
-        burn1.add_control('u1', rate_continuity=True, rate2_continuity=True, units='deg',
+        burn1.add_control('u1', rate_continuity=False, rate2_continuity=False, units='deg',
                           scaler=0.01,
                           rate_continuity_scaler=0.001, rate2_continuity_scaler=0.001,
                           lower=-30, upper=30, targets=['u1'])
@@ -218,7 +218,7 @@ class TestTwoBurnOrbitRaiseLinkages(unittest.TestCase):
                         rate_source='at_dot', targets=['accel'], units='DU/TU**2')
         burn2.add_state('deltav', fix_initial=False, fix_final=False,
                         rate_source='deltav_dot', units='DU/TU')
-        burn2.add_control('u1', targets=['u1'], rate_continuity=True, rate2_continuity=True,
+        burn2.add_control('u1', targets=['u1'], rate_continuity=False, rate2_continuity=False,
                           units='deg', scaler=0.01, lower=-30, upper=30)
 
         burn2.add_objective('deltav', loc='final', scaler=1.0)
@@ -233,9 +233,9 @@ class TestTwoBurnOrbitRaiseLinkages(unittest.TestCase):
 
         # Link Phases
         traj.link_phases(phases=['burn1', 'coast', 'burn2'],
-                         vars=['time', 'r', 'vr', 'vt', 'deltav'])
+                         vars=['time', 'r', 'vr', 'vt', 'deltav'], scaler=1.0E-6, linear=True)
         traj.link_phases(phases=['burn1', 'burn2'], vars=['accel'],
-                         locs=('final', 'initial'))
+                         locs=('final', 'initial'), scaler=1.0E-6, linear=True)
 
         # Finish Problem Setup
         p.model.linear_solver = om.DirectSolver()
diff --git a/dymos/trajectory/test/test_trajectory.py b/dymos/trajectory/test/test_trajectory.py
index 528d66b6..254dbad6 100644
--- a/dymos/trajectory/test/test_trajectory.py
+++ b/dymos/trajectory/test/test_trajectory.py
@@ -358,7 +358,7 @@ class TestLinkages(unittest.TestCase):
         burn2.add_state('deltav', fix_initial=False, fix_final=False,
                         rate_source='deltav_dot', units='DU/TU')
         burn2.add_control('u1', units='deg', scaler=0.01, lower=-30, upper=30)
-        burn2.add_parameter('c', opt=False, val=1.5, units='DU/TU')
+        burn2.add_parameter('c', opt=True, val=1.5, units='DU/TU')
 
         burn2.add_objective('deltav', loc='final')
 
@@ -1330,12 +1330,12 @@ class TestInvalidLinkages(unittest.TestCase):
                          vars=['time', 'r', 'theta', 'vr', 'vt', 'deltav'])
 
         with self.assertWarns(UserWarning) as w:
-            traj.add_linkage_constraint(phase_a='burn1', phase_b='burn2', var_a='accel', var_b='accel',
-                                        lower=-5, upper=5, ref0=-5, ref=5, linear=True, connected=True)
+            traj.add_linkage_constraint(phase_a='burn1', phase_b='burn2', var_a='accel', var_b='accel', connected=True,
+                                        units='DU/TU**2', lower=-5, upper=5, ref0=-5, ref=5, linear=True)
 
         expected_warning = 'Invalid option in linkage between burn1:accel and burn2:accel in trajectory traj. ' \
                            'The following options for add_linkage_constraint were specified but ' \
-                           'not valid when option \'connected\' is True: lower upper ref0 ref linear'
+                           'not valid when option \'connected\' is True: lower upper ref0 ref units linear'
 
         self.assertEqual(expected_warning, str(w.warning))
 
diff --git a/joss/test/test_brachistochrone_for_joss.py b/joss/test/test_brachistochrone_for_joss.py
index db4ef378..42351267 100644
--- a/joss/test/test_brachistochrone_for_joss.py
+++ b/joss/test/test_brachistochrone_for_joss.py
@@ -86,13 +86,13 @@ class TestBrachistochroneForJOSS(unittest.TestCase):
         # States and controls here use a linearly interpolated
         # initial guess along the trajectory.
         p.set_val('traj.phase0.states:x',
-                  phase.interpolate(ys=[0, 10], nodes='state_input'),
+                  phase.interp('x', ys=[0, 10]),
                   units='m')
         p.set_val('traj.phase0.states:y',
-                  phase.interpolate(ys=[10, 5], nodes='state_input'),
+                  phase.interp('y', ys=[10, 5]),
                   units='m')
         p.set_val('traj.phase0.states:v',
-                  phase.interpolate(ys=[0, 5], nodes='state_input'),
+                  phase.interp('v', ys=[0, 5]),
                   units='m/s')
         # constant initial guess for control
         p.set_val('traj.phase0.controls:theta', 90, units='deg')