Allow boundary constraints to be defined as mathematical expressions

[robfalck]

Summary of Issue

Allow users to define boundary constraints as mathematical expressions.

Issue Type

• [ ] Bug
• [x] Enhancement
• [ ] Docs
• [ ] Miscellaneous

Blocked By

• OpenMDAO/POEMs#94
• #714

Description

Often, a user wants a boundary constraint to be defined as a relationship between two or more states, controls, parameters, or time.

The equals, lower, and upper arguments to constraints must be numeric values defined at setup time (they cannot be outputs of the ODE, state variable values, etc.

The current API

Currently, the user can implement this sort of feature in a few ways.

1. The user can define a new ODE output that captures the value to be constrained.

This works, but modifying an otherwise modular ODE is not always desirable. Some users will likely prefer that their ODE's be defined in read-only files and not changed simply to provide constraints for one particular use case. Also, defining an ODE output to capture a boundary constraint means computing that output at many nodes where it is not needed.

2. The user can add an ExecComp downstream of the ODE and use an ordinary OpenMDAO constraint.

This works, but if there is a complication that might trip up new users. Parameters are now inputs to phases and cannot be connected to downstream targets. Instead, that downstream target must have its input promoted to the same path as the parameter.

Proposed change

The new change will allow the user to define equations to be constrained as boundary constraints. For example:

phase.add_boundary_constraint(expr='remaining_prop = mass - m_dry', loc='final', equals=0)

• The name field of the current add_boundary_constraint method will be deprecated in favor of expr which will allow either plain variable names, or equations.

• If the user provides a mathematic expression without an equal sign, Dymos won't know what to name the boundary constraint output. The boundary constraint name will use the variable name of the left side of the equal sign as the constraint name. Dymos will do some checking to make sure that the given expression is either a plain variable name, or a valid expression including an equal sign.

• The expression itself will just be evaluated and compared to lower, upper, and equals. The use of logical comparison operators in the expression is not allowed (no <, >, etc)

• add_boundary_constraint will accept kwargs that are passed to the exec_comp's add_expr method. This is used to mass variable metadata.

Implementation details

• The current boundary constraint component provides a passthru that takes an input, echoes it to an output, and adds a constraint to that output.

• The new implementation will be an ExecComp.

• When a plain variable name is used, the equation will use the form {loc)_value:{name} = {loc}_value_in:{name}

• When a equation is used, each input will be renamed as {loc}_value_in:{name} and the output will be renamed as {loc)_value:{name}.

• Phase boundary constraints will not have access to trajectory parameters, those parameters must be passed to the phase as a phase parameter.

• Currently, boundary-constrained variables are automatically included in the timeseries output. This should not happen for those that are declared as equations.

• We should check the validity of the given expression in add_boundary_constraint if at all possible, rather than deferring to setup. At least check if it is a valid equaiton assuming that all of the given variable names exist.

Example

Consider the cannonball design example. The ODE in this problem includes a calculation of kinetic energy so that it may be limited at the initial time. Instead, we can use an expression in the boundary constraint to limit kinetic energy.

    def test_cannonball_design_boundary_constrint_expression(self):
        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_phase import CannonballPhase

        p = om.Problem(model=om.Group())

        p.driver = om.pyOptSparseDriver()
        p.driver.options['optimizer'] = 'SLSQP'
        p.driver.declare_coloring()

        external_params = p.model.add_subsystem('external_params', om.IndepVarComp())

        external_params.add_output('radius', val=0.10, units='m')
        external_params.add_output('dens', val=7.87, units='g/cm**3')

        external_params.add_design_var('radius', lower=0.01, upper=0.10, ref0=0.01, ref=0.10)

        p.model.add_subsystem('size_comp', CannonballSizeComp())

        traj = p.model.add_subsystem('traj', dm.Trajectory())

        transcription = dm.Radau(num_segments=5, order=3, compressed=True)
        ascent = CannonballPhase(transcription=transcription)

        ascent = traj.add_phase('ascent', ascent)

        # All initial states except flight path angle are fixed
        # Final flight path angle is fixed (we will set it to zero so that the phase ends at apogee)
        ascent.set_time_options(fix_initial=True, duration_bounds=(1, 100), duration_ref=100, units='s')
        ascent.set_state_options('r', fix_initial=True, fix_final=False)
        ascent.set_state_options('h', fix_initial=True, fix_final=False)
        ascent.set_state_options('gam', fix_initial=False, fix_final=True)
        ascent.set_state_options('v', fix_initial=False, fix_final=False)

        ascent.add_parameter('S', targets=['aero.S'], units='m**2')
        ascent.add_parameter('mass', targets=['eom.m', 'kinetic_energy.m'], units='kg')

        # Limit the muzzle energy
        ascent.add_boundary_constraint(expr='ke = 0.5 * mass * v**2', loc='initial',
                                       upper=400000, lower=0, ref=100000)

        # Second Phase (descent)
        transcription = dm.GaussLobatto(num_segments=5, order=3, compressed=True)
        descent = CannonballPhase(transcription=transcription)

        traj.add_phase('descent', descent)

        # All initial states and time are free (they will be linked to the final states of ascent.
        # Final altitude is fixed (we will set it to zero so that the phase ends at ground impact)
        descent.set_time_options(initial_bounds=(.5, 100), duration_bounds=(.5, 100),
                                 duration_ref=100, units='s')
        descent.add_state('r', )
        descent.add_state('h', fix_initial=False, fix_final=True)
        descent.add_state('gam', fix_initial=False, fix_final=False)
        descent.add_state('v', fix_initial=False, fix_final=False)

        descent.add_parameter('S', targets=['aero.S'], units='m**2')
        descent.add_parameter('mass', targets=['eom.m', 'kinetic_energy.m'], units='kg')

        descent.add_objective('r', loc='final', scaler=-1.0)

        # Add internally-managed design parameters to the trajectory.
        traj.add_parameter('CD',
                           targets={'ascent': ['aero.CD'], 'descent': ['aero.CD']},
                           val=0.5, units=None, opt=False)
        traj.add_parameter('CL',
                           targets={'ascent': ['aero.CL'], 'descent': ['aero.CL']},
                           val=0.0, units=None, opt=False)
        traj.add_parameter('T',
                           targets={'ascent': ['eom.T'], 'descent': ['eom.T']},
                           val=0.0, units='N', opt=False)
        traj.add_parameter('alpha',
                           targets={'ascent': ['eom.alpha'], 'descent': ['eom.alpha']},
                           val=0.0, units='deg', opt=False)

        # Add externally-provided design parameters to the trajectory.
        # In this case, we connect 'm' to pre-existing input parameters named 'mass' in each phase.
        traj.add_parameter('m', units='kg', val=1.0,
                           targets={'ascent': 'mass', 'descent': 'mass'})

        # In this case, by omitting targets, we're connecting these parameters to parameters
        # with the same name in each phase.
        traj.add_parameter('S', units='m**2', val=0.005)

        # Link Phases (link time and all state variables)
        traj.link_phases(phases=['ascent', 'descent'], vars=['*'])

        # Issue Connections
        p.model.connect('external_params.radius', 'size_comp.radius')
        p.model.connect('external_params.dens', 'size_comp.dens')

        p.model.connect('size_comp.mass', 'traj.parameters:m')
        p.model.connect('size_comp.S', 'traj.parameters:S')

        # A linear solver at the top level can improve performance.
        p.model.linear_solver = om.DirectSolver()

        # Finish Problem Setup
        p.setup()

        # Set Initial Guesses
        p.set_val('external_params.radius', 0.05, units='m')
        p.set_val('external_params.dens', 7.87, units='g/cm**3')

        p.set_val('traj.parameters:CD', 0.5)
        p.set_val('traj.parameters:CL', 0.0)
        p.set_val('traj.parameters:T', 0.0)

        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.interpolate(ys=[0, 100], nodes='state_input'))
        p.set_val('traj.ascent.states:h', ascent.interpolate(ys=[0, 100], nodes='state_input'))
        p.set_val('traj.ascent.states:v', ascent.interpolate(ys=[200, 150], nodes='state_input'))
        p.set_val('traj.ascent.states:gam', ascent.interpolate(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.interpolate(ys=[100, 200], nodes='state_input'))
        p.set_val('traj.descent.states:h', descent.interpolate(ys=[100, 0], nodes='state_input'))
        p.set_val('traj.descent.states:v', descent.interpolate(ys=[150, 200], nodes='state_input'))
        p.set_val('traj.descent.states:gam', descent.interpolate(ys=[0, -45], nodes='state_input'),
                  units='deg')

        dm.run_problem(p)

Environment

N/A