Added ability to add expressions as constraints

dymos/examples/brachistochrone/test/test_brachistochrone_path_constraint.py

@@ -0,0 +1,75 @@
+import os
+import unittest
+
+import matplotlib
+import openmdao.api as om
+import matplotlib.pyplot as plt
+import dymos as dm
+
+from openmdao.utils.testing_utils import use_tempdirs
+from openmdao.utils.assert_utils import assert_near_equal
+from dymos.examples.brachistochrone.brachistochrone_ode import BrachistochroneODE
+
+matplotlib.use('Agg')
+plt.style.use('ggplot')
+
+
+@use_tempdirs
+class TestBrachistochroneStaticGravity(unittest.TestCase):
+
+    def _make_problem(self, tx):
+        p = om.Problem(model=om.Group())
+
+        p.driver = om.ScipyOptimizeDriver()
+        p.driver.declare_coloring()
+
+        phase = dm.Phase(ode_class=BrachistochroneODE,
+                         transcription=tx)
+
+        p.model.add_subsystem('phase0', phase)
+
+        phase.set_time_options(fix_initial=True, duration_bounds=(.5, 10))
+
+        phase.set_state_options('x', fix_initial=True)
+        phase.set_state_options('y', fix_initial=True)
+        phase.set_state_options('v', fix_initial=True)
+
+        phase.add_control('theta', continuity=True, rate_continuity=True, opt=True,
+                          units='deg', lower=0.01, upper=179.9, ref=1, ref0=0)
+
+        phase.add_boundary_constraint('x', loc='final', equals=10.0)
+        phase.add_path_constraint('pc = y-x/2-1', lower=0.0)
+
+        phase.add_parameter('g', opt=False, units='m/s**2', val=9.80665, include_timeseries=True)
+
+        phase.add_objective('time_phase', loc='final', scaler=10)
+
+        p.model.options['assembled_jac_type'] = 'csc'
+        p.model.linear_solver = om.DirectSolver()
+        p.setup(check=True)
+
+        p['phase0.t_initial'] = 0.0
+        p['phase0.t_duration'] = 2.0
+
+        p['phase0.states:x'] = phase.interp('x', [0, 10])
+        p['phase0.states:y'] = phase.interp('y', [10, 5])
+        p['phase0.states:v'] = phase.interp('v', [0, 9.9])
+        p[f'phase0.controls:theta'] = phase.interp('theta', [5, 100])
+        p['phase0.parameters:g'] = 9.80665
+        return p
+
+    def tearDown(self):
+        for filename in ['total_coloring.pkl', 'SLSQP.out', 'SNOPT_print.out', 'SNOPT_summary.out']:
+            if os.path.exists(filename):
+                os.remove(filename)
+
+    def test_brachistochrone_path_constraint(self):
+        prob = self._make_problem(tx=dm.Radau(num_segments=5, order=3, compressed=True))
+        prob.run_driver()
+        yf = prob.get_val('phase0.timeseries.states:y')[-1]
+
+        assert_near_equal(yf, 6)
+
+
+if __name__ == '__main__':  # pragma: no cover
+    unittest.main()

dymos/examples/cannonball/test/test_cannonball_expr_constraint.py

@@ -0,0 +1,124 @@
+import unittest
+import matplotlib.pyplot as plt
+
+from openmdao.utils.testing_utils import use_tempdirs
+import openmdao.api as om
+from openmdao.utils.assert_utils import assert_near_equal
+
+import dymos as dm
+from dymos.examples.cannonball.size_comp import CannonballSizeComp
+from dymos.examples.cannonball.cannonball_ode import CannonballODE
+
+
+plt.switch_backend('Agg')
+
+
+@use_tempdirs
+class TestCannonballBoundaryConstraint(unittest.TestCase):
+
+    def test_cannonball_design_boundary_constraint_expression(self):
+
+        p = om.Problem(model=om.Group())
+
+        p.driver = om.ScipyOptimizeDriver()
+        p.driver.options['optimizer'] = 'SLSQP'
+        p.driver.declare_coloring()
+
+        p.model.add_subsystem('size_comp', CannonballSizeComp(),
+                              promotes_inputs=['radius', 'dens'])
+        p.model.set_input_defaults('dens', val=7.87, units='g/cm**3')
+        p.model.add_design_var('radius', lower=0.01, upper=0.10,
+                               ref0=0.01, ref=0.10, units='m')
+
+        traj = p.model.add_subsystem('traj', dm.Trajectory())
+
+        transcription = dm.Radau(num_segments=5, order=3, compressed=True)
+        ascent = dm.Phase(ode_class=CannonballODE, transcription=transcription)
+
+        ascent = traj.add_phase('ascent', ascent)
+
+        ascent.set_time_options(fix_initial=True, duration_bounds=(1, 100),
+                                duration_ref=100, units='s')
+
+        ascent.set_time_options(fix_initial=True, duration_bounds=(1, 100), duration_ref=100, units='s')
+        ascent.add_state('r', fix_initial=True, fix_final=False, rate_source='r_dot', units='m')
+        ascent.add_state('h', fix_initial=True, fix_final=False, units='m', rate_source='h_dot')
+        ascent.add_state('gam', fix_initial=False, fix_final=True, units='rad', rate_source='gam_dot')
+        ascent.add_state('v', fix_initial=False, fix_final=False, units='m/s', rate_source='v_dot')
+
+        ascent.add_parameter('S', units='m**2', dynamic=False)
+        ascent.add_parameter('m', units='kg', dynamic=False)
+
+        # Limit the muzzle energy
+        ascent.add_boundary_constraint('ke = 0.5 * m * v**2', loc='initial',
+                                       upper=400000, lower=0, ref=100000)
+
+        # Second Phase (descent)
+        transcription = dm.GaussLobatto(num_segments=5, order=3, compressed=True)
+        descent = dm.Phase(ode_class=CannonballODE, transcription=transcription)
+
+        traj.add_phase('descent', descent)
+
+        descent.set_time_options(initial_bounds=(.5, 100), duration_bounds=(.5, 100),
+                                 duration_ref=100, units='s')
+        descent.add_state('r', units='m', rate_source='r_dot')
+        descent.add_state('h', units='m', rate_source='h_dot', fix_initial=False, fix_final=True)
+        descent.add_state('gam', units='rad', rate_source='gam_dot', fix_initial=False, fix_final=False)
+        descent.add_state('v', units='m/s', rate_source='v_dot', fix_initial=False, fix_final=False)
+
+        descent.add_parameter('S', units='m**2', dynamic=False)
+        descent.add_parameter('m', units='kg', dynamic=False)
+
+        descent.add_objective('r', loc='final', scaler=-1.0)
+
+        traj.add_parameter('CD',
+                           targets={'ascent': ['CD'], 'descent': ['CD']},
+                           val=0.5, units=None, opt=False, dynamic=False)
+
+        traj.add_parameter('m', units='kg', val=1.0,
+                           targets={'ascent': 'm', 'descent': 'm'}, dynamic=False)
+
+        traj.add_parameter('S', units='m**2', val=0.005, dynamic=False)
+
+        # Link Phases (link time and all state variables)
+        traj.link_phases(phases=['ascent', 'descent'], vars=['*'])
+
+        p.model.connect('size_comp.mass', 'traj.parameters:m')
+        p.model.connect('size_comp.S', 'traj.parameters:S')
+
+        p.model.linear_solver = om.DirectSolver()
+
+        # Finish Problem Setup
+        p.setup()
+
+        # Set constants and initial guesses
+        p.set_val('radius', 0.05, units='m')
+        p.set_val('dens', 7.87, units='g/cm**3')
+
+        p.set_val('traj.parameters:CD', 0.5)
+
+        p.set_val('traj.ascent.t_initial', 0.0)
+        p.set_val('traj.ascent.t_duration', 10.0)
+
+        p.set_val('traj.ascent.states:r', ascent.interp(ys=[0, 100], nodes='state_input'))
+        p.set_val('traj.ascent.states:h', ascent.interp(ys=[0, 100], nodes='state_input'))
+        p.set_val('traj.ascent.states:v', ascent.interp(ys=[200, 150], nodes='state_input'))
+        p.set_val('traj.ascent.states:gam', ascent.interp(ys=[25, 0], nodes='state_input'), units='deg')
+
+        p.set_val('traj.descent.t_initial', 10.0)
+        p.set_val('traj.descent.t_duration', 10.0)
+
+        p.set_val('traj.descent.states:r', descent.interp(ys=[100, 200], nodes='state_input'))
+        p.set_val('traj.descent.states:h', descent.interp(ys=[100, 0], nodes='state_input'))
+        p.set_val('traj.descent.states:v', descent.interp(ys=[150, 200], nodes='state_input'))
+        p.set_val('traj.descent.states:gam', descent.interp(ys=[0, -45], nodes='state_input'), units='deg')
+
+        # Run the optimization and final explicit simulation
+        dm.run_problem(p)
+
+        assert_near_equal(p.get_val('traj.descent.states:r')[-1],
+                          3183.25, tolerance=1.0E-2)
+
+
+if __name__ == '__main__':
+    unittest.main()

dymos/phase/options.py

@@ -688,6 +688,10 @@ def __init__(self, read_only=False):
                      desc='If True, the given indices will be treated as indices into a C-order flattened array based '
                           'on the shaped of the constrained variable at a point in time.')
 
+        self.declare(name='is_expr', types=bool, default=False,
+                     desc='If True, the given constraint is an expression that must be evaluated rather than a'
+                          ' single variable.')
+
 
 class TimeseriesOutputOptionsDictionary(om.OptionsDictionary):
     """

dymos/phase/phase.py

@@ -22,7 +22,7 @@
 from ..utils.indexing import get_constraint_flat_idxs
 from ..utils.introspection import configure_time_introspection, _configure_constraint_introspection, \
     configure_controls_introspection, configure_parameters_introspection, \
-    configure_timeseries_output_introspection, classify_var, get_promoted_vars
+    configure_timeseries_output_introspection, classify_var, configure_timeseries_expr_introspection
 from ..utils.misc import _unspecified
 from ..utils.lgl import lgl
 
@@ -1083,9 +1083,9 @@ def add_boundary_constraint(self, name, loc, constraint_name=None, units=None,
         loc : str
             The location of the boundary constraint ('initial' or 'final').
         constraint_name : str or None
-            The name of the variable as provided to the boundary constraint comp.  By
-            default this is the last element in `name` when split by dots.  The user may
-            override the constraint name if splitting the path causes name collisions.
+            The name of the boundary constraint. By default, this is the left-hand side of
+            the given expression or "var_constraint" if var is a single variable. The user may
+            override the constraint name if desired.
         units : str or None
             The units in which the boundary constraint is to be applied.  If None, use the
             units associated with the constrained output.  If provided, must be compatible with
@@ -1126,7 +1126,19 @@ def add_boundary_constraint(self, name, loc, constraint_name=None, units=None,
             raise ValueError(f'Invalid boundary constraint location "{loc}". Must be '
                              '"initial" or "final".')
 
-        if constraint_name is None:
+        expr_operators = ['(', '+', '-', '/', '*', '&', '%']
+        if '=' in name:
+            is_expr = True
+        elif '=' not in name and any(opr in name for opr in expr_operators):
+            raise ValueError(f'The expression provided {name} has invalid format. '
+                             'Expression may be a single variable or an equation'
+                             'of the form "constraint_name = func(vars)"')
+        else:
+            is_expr = False
+
+        if is_expr:
+            constraint_name = name.split('=')[0].strip()
+        elif constraint_name is None:
             constraint_name = name.rpartition('.')[-1]
 
         bc_list = self._initial_boundary_constraints if loc == 'initial' else self._final_boundary_constraints
@@ -1137,6 +1149,13 @@ def add_boundary_constraint(self, name, loc, constraint_name=None, units=None,
             raise ValueError(f'Cannot add new {loc} boundary constraint for variable `{name}` and indices {indices}. '
                              f'One already exists.')
 
+        existing_bc = [bc for bc in bc_list if bc['name'] == constraint_name and
+                       bc['indices'] is None and indices is None]
+
+        if existing_bc:
+            raise ValueError(f'Cannot add new {loc} boundary constraint named `{constraint_name}` and indices{indices}.'
+                             f' `{constraint_name}` is already in use as a {loc} boundary constraint')
+
         bc = ConstraintOptionsDictionary()
         bc_list.append(bc)
 
@@ -1154,6 +1173,7 @@ def add_boundary_constraint(self, name, loc, constraint_name=None, units=None,
         bc['linear'] = linear
         bc['units'] = units
         bc['flat_indices'] = flat_indices
+        bc['is_expr'] = is_expr
 
         # Automatically add the requested variable to the timeseries outputs if it's an ODE output.
         var_type = self.classify_var(name)
@@ -1211,7 +1231,19 @@ def add_path_constraint(self, name, constraint_name=None, units=None, shape=None
             If True, treat indices as flattened C-ordered indices of elements to constrain at each given point in time.
             Otherwise, indices should be a tuple or list giving the elements to constrain at each point in time.
         """
-        if constraint_name is None:
+        expr_operators = ['(', '+', '-', '/', '*', '&', '%']
+        if '=' in name:
+            is_expr = True
+        elif '=' not in name and any(opr in name for opr in expr_operators):
+            raise ValueError(f'The expression provided {name} has invalid format. '
+                             'Expression may be a single variable or an equation'
+                             'of the form "constraint_name = func(vars)"')
+        else:
+            is_expr = False
+
+        if is_expr:
+            constraint_name = name.split('=')[0].strip()
+        elif constraint_name is None:
             constraint_name = name.rpartition('.')[-1]
 
         existing_pc = [pc for pc in self._path_constraints
@@ -1238,6 +1270,7 @@ def add_path_constraint(self, name, constraint_name=None, units=None, shape=None
         pc['linear'] = linear
         pc['units'] = units
         pc['flat_indices'] = flat_indices
+        pc['is_expr'] = is_expr
 
         # Automatically add the requested variable to the timeseries outputs if it's an ODE output.
         var_type = self.classify_var(name)
@@ -1758,7 +1791,9 @@ def configure(self):
 
         _configure_constraint_introspection(self)
 
-        transcription._configure_boundary_constraints(self)
+        configure_timeseries_expr_introspection(self)
+
+        transcription.configure_boundary_constraints(self)
 
         transcription.configure_path_constraints(self)
 

dymos/phase/test/test_timeseries.py

@@ -483,7 +483,7 @@ def test_invalid_expr_var(self):
 
         expected = "phase0: Timeseries expression `k=units` uses the following disallowed variable " \
                    "name in its definition: 'units'."
-        self.assertEqual(str(e.exception.__cause__), expected)
+        self.assertEqual(str(e.exception), expected)
 
     def test_input_units(self):
         p = om.Problem(model=om.Group())
@@ -549,7 +549,7 @@ def test_output_units(self):
         phase.add_timeseries_output('sin_theta2=sin(theta)', sin_theta2={'units': 'unitless'})
         phase.add_timeseries_output('sin_theta3=sin(theta)', shape=(1,), sin_theta3={'shape': (1,)})
 
-        expected = 'phase0: User-provided shape for timeseries expression output `sin_theta3`using ' \
+        expected = 'phase0: User-provided shape for timeseries expression output `sin_theta3` using ' \
                    'both the\n`shape` keyword argument and the `sin_theta3` keyword argument dictionary.\n' \
                    'The latter will take precedence.'
 

dymos/trajectory/trajectory.py

@@ -1112,8 +1112,11 @@ def _print_constraints(phase, outstream):
             for loc, d in ds.items():
                 str_loc = f'[{loc}]'
                 for options in d:
-                    expr = options['name']
-                    _, shape, units, linear = tx._get_objective_src(expr, loc, phase, ode_outputs=ode_outputs)
+                    if options['is_expr']:
+                        name = options['constraint_name']
+                    else:
+                        name = options['name']
+                    _, shape, units, linear = tx._get_objective_src(name, loc, phase, ode_outputs=ode_outputs)
 
                     equals = options['equals']
                     lower = options['lower']
@@ -1146,11 +1149,11 @@ def _print_constraints(phase, outstream):
                         str_upper = ''
 
                     if equals is not None:
-                        printer(f'{2 * indent}{str_loc:<10s}{str_equals} == {expr} [{str_units}]',
+                        printer(f'{2 * indent}{str_loc:<10s}{str_equals} == {name} [{str_units}]',
                                 file=outstream)
                     else:
                         printer(
-                            f'{2 * indent}{str_loc:<10s}{str_lower} {expr} {str_upper} [{str_units}]',
+                            f'{2 * indent}{str_loc:<10s}{str_lower} {name} {str_upper} [{str_units}]',
                             file=outstream)
 
         for phase_name, phs in self._phases.items():

dymos/transcriptions/pseudospectral/pseudospectral_base.py

@@ -602,6 +602,11 @@ def _get_objective_src(self, var, loc, phase, ode_outputs=None):
             units = control_rate_units
             linear = False
             constraint_path = f'polynomial_control_rates:{var}'
+        elif var_type == 'timeseries_exec_comp_output':
+            shape = (1,)
+            units = None
+            constraint_path = f'timeseries.timeseries_exec_comp.{var}'
+            linear = False
         else:
             # Failed to find variable, assume it is in the ODE. This requires introspection.
             constraint_path = f'{self._rhs_source}.{var}'

dymos/transcriptions/solve_ivp/solve_ivp.py

@@ -460,7 +460,7 @@ def configure_path_constraints(self, phase):
         """
         pass
 
-    def _configure_boundary_constraints(self, phase):
+    def configure_boundary_constraints(self, phase):
         """
         Not used in SolveIVP.
 
@@ -686,4 +686,8 @@ def _get_num_timeseries_nodes(self):
         int
             The number of nodes in the default timeseries for this transcription.
         """
-        return self.grid_data.num_segments * self.options['output_nodes_per_seg']
+        if self.options['output_nodes_per_seg'] is None:
+            output_nodes = self.grid_data.num_segments
+        else:
+            output_nodes = self.grid_data.num_segments * self.options['output_nodes_per_seg']
+        return output_nodes