Skip to content

Defining ODEs

Optimal control problems contain ordinary differential equations (ODE) or differential algebraic equations (DAE) that dictate the evolution of the state of the system. Typically this evolution occurs in time, and the ODE represents equations of motion (EOM). The equations of motion can define a variety of systems, not just mechanical ones. In other fields they are sometimes referred to as process equations.

\begin{align} \dot{\bar{x}} = f_{ode}(\bar{x},t,\bar{u},\bar{d}) \end{align}

To represent EOM, Dymos uses a standard OpenMDAO System (a Group or Component). This System takes some set of variables as input and computes outputs that include the time-derivatives of the state variables \bar{x}. The ODE may also be a function of the current time t.

Finally, the dynamics may be subject to some set of controllable parameters. In Dymos these are broken into the dynamic controls \bar{u} and the static parameters \bar{d}.

System Options

ODE Systems in Dymos need provide values at numerous points in time we call nodes. For performance reasons, it's best if it can do so using vectorized mathematical operations to compute the values simultaneously rather than using the loop to perform the calculation at each node. Different optimal control transcriptions will need to have computed ODE values at different numbers of nodes, so each ODE system in Dymos is required to support the option num_nodes, which is defined in the initialize method.

ODE system may define initial options as well. Since these options in OpenMDAO are typically provided as arguments to the instantiation of the ODE system, the user has the ability to provide additional input keyword arguments using the ode_init_kwargs option on Phase.

Variables of the Optimal Control Problem

Dymos needs to know how state, time, and control variables are to be connected to the System, and needs to know which outputs to use as state time-derivatives.

Time

The following time options can be set via the set_time_options method of Phase:

Options for Time

Option Default Acceptable Values Acceptable Types Description
duration_adder None N/A ['Number'] Adder for the duration of the integration variable across the phase.
duration_bounds (None, None) N/A ['Iterable'] Tuple of (lower, upper) bounds for the duration of the integration variable across the phase.
duration_ref None N/A ['Number'] Unit-reference value for the duration of the integration variable across the phase.
duration_ref0 None N/A ['Number'] Zero-reference value for the duration of the integration variable across the phase.
duration_scaler None N/A ['Number'] Scalar for the duration of the integration variable across the phase.
duration_val 1.0 N/A ['Number'] Value of the duration of the integration variable across the phase.
fix_duration False [True, False] ['bool'] If True, the phase duration is not a design variable.
fix_initial False [True, False] ['bool'] If True, the initial value of time is not a design variable.
initial_adder None N/A ['Number'] Adder for the initial value of the integration variable.
initial_bounds (None, None) N/A ['Iterable'] Tuple of (lower, upper) bounds for the integration variable at the start of the phase.
initial_ref None N/A ['Number'] Unit-reference value for the initial value of the integration variable.
initial_ref0 None N/A ['Number'] Zero-reference value for the initial value of the integration variable.
initial_scaler None N/A ['Number'] Scalar for the initial value of the integration variable.
initial_val 0.0 N/A ['Number'] Value of the integration variable at the start of the phase.
input_duration False [True, False] ['bool'] If True, the phase duration (t_duration) is expected to be connected to an external output source.
input_initial False [True, False] ['bool'] If True, the initial value of time (t_initial) is expected to be connected to an external output source.
t_duration_targets [] N/A ['Iterable'] targets in the ODE to which the total duration of the phase is connected
t_initial_targets [] N/A ['Iterable'] targets in the ODE to which the initial time of the phase is connected
targets [] N/A ['Iterable'] targets in the ODE to which the integration variable is connected
time_phase_targets [] N/A ['Iterable'] targets in the ODE to which the elapsed duration of the phase is connected
units s N/A ['str'] Units for the integration variable

States

States have the following options set via the add_state and set_state_options methods of Phase:

Options for States

Option Default Acceptable Values Acceptable Types Description
adder None N/A ['Iterable', 'Number'] Adder of the state variable at the discretization nodes. This option is invalid if opt=False.
connected_initial False [True, False] ['bool'] Whether an input is created to pass in the initial state. This may be set by a trajectory that links phases.
continuity True N/A ['bool', 'dict'] Enforce continuity of state values at segment boundaries. This option is invalid if opt=False.
defect_ref None N/A ['Iterable', 'Number'] Unit-reference value of the state defects at the collocation nodes. This option is invalid if opt=False. If provided, this option overrides defect_scaler. If defect_scaler and defect_ref are both None but the state scaler or ref are provided, use those values as the defect scaler or ref.
defect_scaler None N/A ['Iterable', 'Number'] Scaler of the state variable defects at the collocation nodes. If defect_scaler and defect_ref are both None but the state scaler or ref are provided, use those values as the defect scaler or ref. This option is invalid if opt=False.
desc N/A ['str'] description of the state variable
final_bounds None N/A ['Iterable'] Bounds on the value of the state at the end of the phase. This option is invalid if opt=False.
fix_final False N/A ['bool', 'Iterable'] If True, the final value of this state is fixed and not a design variable. If the state variable has a non-scalar shape, this may be an iterable of bool for each index. This option is invalid if opt=False.
fix_initial False N/A ['bool', 'Iterable'] If True, the initial value of this state is fixed and not a design variable. If the state variable has a non-scalar shape, this may be an iterable of bool for each index. This option is invalid if opt=False.
initial_bounds None N/A ['Iterable'] Bounds on the value of the state at the start of the phase. This option is invalid if opt=False.
lower None N/A ['Iterable', 'Number'] Lower bound of the state variable at the discretization nodes. This option is invalid if opt=False.
name Required N/A ['str'] name of ODE state variable
opt True [True, False] ['bool'] If true, the values of this state are a design variable for the optimizer. Otherwise it exists as an unconnected input.
rate_source None N/A ['str'] ODE-path to the derivative of the state variable
ref None N/A ['Iterable', 'Number'] Unit-reference value of the state variable at the discretization nodes. This option is invalid if opt=False.
ref0 None N/A ['Iterable', 'Number'] Zero-reference value of the state variable at the discretization nodes. This option is invalid if opt=False.
scaler None N/A ['Iterable', 'Number'] Scaler of the state variable at the discretization nodes. This option is invalid if opt=False.
shape None N/A ['Iterable'] shape of the state variable, as determined by introspection
solve_segments None [True, False, 'forward', 'backward'] N/A If True (deprecated) or 'forward', collocation defects within eachsegment are solved with a Newton solver by fixing the initial value in thephase (if using compressed transcription) or segment (if not using compressed transcription). This provides a forward shooting (or multiple shooting)method. If 'backward', the final value in the phase or segment is fixedand a solver finds the other ones to mimic reverse propagation. If None, (the default) use the value of solve_segments in the transcription. Set to False to explicitly disable the use of a solver to converge the statetime history.
targets unspecified N/A N/A Targets in the ODE to which the state is connected
units unspecified N/A N/A units in which the state variable is defined
upper None N/A ['Iterable', 'Number'] Upper bound of the state variable at the discretization nodes. This option is invalid if opt=False.
val 1.0 N/A ['Iterable', 'Number'] Default value of the state variable at the discretization nodes

State Discovery

Dymos will also automatically find and add any states that have been declared in components in the ODE. The syntax for declaring them is as follows.

        self.add_output('vdot', val=np.zeros(nn), desc='acceleration magnitude', units='m/s**2',
                        tags=['state_rate_source:v', 'state_units:m/s'])

The state is defined by adding tags to the state rate's output. The tag 'state_rate_source:v' declares that 'v' is the state for which this output ('vdot') is the rate source. You can also optionally use a tag in the format 'state_units:m/s' to define units for that state. If you need to set any other options, then use set_state_options at the phase level.

Controls

Inputs to the ODE which are to be dynamic controls are added via the add_control and set_control_options methods of Phase. The available options are as follows:

Options for Control

Option Default Acceptable Values Acceptable Types Description
adder None N/A ['Iterable', 'Number'] The adder of the control variable at the nodes. This option is invalid if opt=False.
continuity True N/A ['bool', 'dict'] Enforce continuity of control values at segment boundaries. This option is invalid if opt=False.
desc N/A ['str'] The description of the control variable.
dynamic True [True, False] ['bool'] If True, the value of the shape of the parameter will be (num_nodes, ...), allowing the variable to be used as either a static or dynamic control. This impacts the shape of the partial derivatives matrix. Unless a parameter is large and broadcasting a value to each individual node would be inefficient, users should stick to the default value of True.
fix_final False [True, False] ['bool'] If True, the final value of this control is fixed and not a design variable. This option is invalid if opt=False.
fix_initial False [True, False] ['bool'] If True, the initial value of this control is fixed and not a design variable. This option is invalid if opt=False.
lower None N/A ['Iterable', 'Number'] The lower bound of the control variable at the nodes. This option is invalid if opt=False.
name Required N/A ['str'] The name of ODE system parameter to be controlled.
opt True [True, False] ['bool'] If True, the control value will be a design variable for the optimization problem. If False, allow the control to be connected externally.
rate2_continuity False N/A ['bool', 'dict'] Enforce continuity of control second derivatives at segment boundaries. This option is invalid if opt=False.
rate2_continuity_scaler 1.0 N/A ['Number'] Scaler of the dimensionless rate continuity constraint at segment boundaries. This option is invalid if opt=False.
rate2_targets unspecified N/A N/A The targets in the ODE to which the control 2nd derivative is connected
rate_continuity True N/A ['bool', 'dict'] Enforce continuity of control first derivatives in dimensionless time at segment boundaries. This option is invalid if opt=False.
rate_continuity_scaler 1.0 N/A ['Number'] Scaler of the dimensionless rate continuity constraint at segment boundaries. This option is invalid if opt=False.
rate_targets unspecified N/A N/A The targets in the ODE to which the control rate is connected
ref None N/A ['Iterable', 'Number'] The unit-reference value of the control variable at the nodes. This option is invalid if opt=False.
ref0 None N/A ['Iterable', 'Number'] The zero-reference value of the control variable at the nodes. This option is invalid if opt=False.
scaler None N/A ['Iterable', 'Number'] The scaler of the control variable at the nodes. This option is invalid if opt=False.
shape None N/A ['Iterable'] The shape of the control variable at each point in time.
targets unspecified N/A N/A Targets in the ODE to which the state is connected
units unspecified N/A N/A The units in which the control variable is defined.
upper None N/A ['Iterable', 'Number'] The upper bound of the control variable at the nodes. This option is invalid if opt=False.
val [0.] N/A ['Iterable', 'ndarray', 'Number'] The default value of the control variable at the control discretization nodes.

Parameters

Inputs to the ODE which are non-time-varying can be specified using the add_parameter method of Phase. Parameters may be used as design variables (by setting opt = True), or they may be connected to an output external to the Phase or Trajectory. The available options are as follows:

Options for Parameters

Option Default Acceptable Values Acceptable Types Description
adder None N/A ['Iterable', 'Number'] The adder of the parameter. This option is invalid if opt=False.
desc N/A ['str'] The description of the parameter.
dynamic True [True, False] ['bool'] True if this parameter can be used as a dynamic control, else False
include_timeseries True [True, False] ['bool'] True if the static parameters should be included in output timeseries, else False
lower None N/A ['Iterable', 'Number'] The lower bound of the parameter. This option is invalid if opt=False.
name Required N/A ['str'] The name of ODE system parameter to be set via parameter.
opt True [True, False] ['bool'] If True, the control value will be a design variable for the optimization problem. If False, allow the control to be connected externally.
ref None N/A ['Iterable', 'Number'] The unit-reference value of the parameter. This option is invalid if opt=False.
ref0 None N/A ['Iterable', 'Number'] The zero-reference value of the parameter. This option is invalid if opt=False.
scaler None N/A ['Iterable', 'Number'] The scaler of the parameter. This option is invalid if opt=False.
shape unspecified N/A N/A The shape of the parameter.
targets unspecified N/A N/A Targets in the ODE to which the state is connected
units unspecified N/A N/A The units in which the parameter is defined.
upper None N/A ['Iterable', 'Number'] The upper bound of the parameter. This option is invalid if opt=False.
val [0.] N/A ['Iterable', 'ndarray', 'Number'] The default value of the parameter in the phase.

By default, Dymos assumes that the ODE is defined such that a value of the parameter is provided at each node. This makes it easier to define the partials for the ODE in a way such that some inputs may be used interchangeably as either controls or parameters. If an ODE input is only to be used as a static variable (and sized as (1,) instead of by the number of nodes, then the user may specify option dynamic = False to override this behavior.

Vectorizing the ODE

In addition to specifying the ODE Options, a system used as an ODE is required to have a metadata entry called num_nodes. When evaluating the dynamics, these systems will receive time, states, controls, and other inputs as vectorized values, where item in the vector represents the variable value at a discrete time in the trajectory.

The nodes are discretization or collocation locations in the polynomials which represent each segment. The number of nodes in a given phase (to be evaluated by the ODE system) is determined by the number of segments in the phase and the polynomial order in each segment. When Dymos instantiates the ODE system it provides the total number of nodes at which evaluation is required to the ODE system. Thus, at a minimum, the initialize method of components for an ODE system typically look something like this:

The inputs and outputs of the system are expected to provide a scalar or dimensioned value at each node. Vectorization of the component via numpy adds a significant performance increase compared to using a for loop to cycle through calculations at each node. It's important to remember that vectorized data is going to be coming in, this is especially important for defining partials. From the perspective of the ODE system, the outputs at some time t only depend on the values of the input variables at time t. When the output variables are scalar at any given time, this results in components whose Jacobian matrices are diagonal. This large degree of sparsity leads to computational advantages when using sparse-aware optimizers like SNOPT. Users should declare the partial derivatives of their components to be sparse (by specifying nonzero rows and columns) whenever possible.

    class MyODESystem(ExplicitComponent):

        def initialize(self):
            self.metadata.declare('num_nodes', types=int)

For example, if MyODEComponent is to compute the linear function y = a * x + b then the setup, compute, and compute partials methods might look like this:

    def setup(self):
        nn = self.metadata['num_nodes']

        self.add_input('a', shape=(nn,), units='m')
        self.add_input('x', shape=(nn,), units='1/s')
        self.add_input('b', shape=(nn,), units='m/s')

        self.add_output('y', shape=(nn,), units='m/s')

        r = c = np.arange(nn)
        self.declare_partials(of='y', wrt='a', rows=r, cols=c)
        self.declare_partials(of='y', wrt='x', rows=r, cols=c)
        self.declare_partials(of='y', wrt='b', rows=r, cols=c, val=1.0)

    def compute(self, inputs, outputs):
        a = inputs['a']
        x = inputs['x']
        b = inputs['b']

        outputs['y'] = a * x + b

    def compute_partials(self, inputs, outputs, partials):
        a = inputs['a']
        x = inputs['x']
        b = inputs['b']

        partials['y', 'a'] = x
        partials['y', 'x'] = a

A few things to note here. We can use the shape or val argument of add_input and add_output to dimension each variable. In this case each variable is assumed to be a scalar at each point in time (each node). We use the rows and cols arguments of declare_partials to provide the sparsity. Here using arange(nn) for both gives us a diagonal jacobian with nn rows and nn columns. Since the number of nonzero values in the jacobian is nn, we only need to provide nn values in the compute_partials method. It will automatically fill them into the sparse jacobian matrix, in row-major order.

In this example, the partial of y with respect to b is linear, so we can simply provide it in the declare_partials call rather than reassigning it every time compute_partials is called. The provided scalar value of 1.0 is broadcast to all nn values of the Jacobian matrix.

Dimensioned Inputs and Outputs

The above example assumes all inputs and outputs are scalar at each node. Sometimes the user may encounter a situation in which the inputs and/or outputs are vectors, matrices, or tensors at each node. In this case the dimension of the variable is num_nodes, with the dimension of the variable at a single node filling out the remaining indices. A 3-vector is thus dimensioned (num_nodes, 3), while a 3 x 3 matrix would be sized (num_nodes, 3, 3).

Non-Vector Inputs

Declaring inputs as vectors means that they have the potential to be used either as parameters or as dynamic controls which can assume a different value at each node. For some quantities, such as gravitational acceleration in the Brachistochrone example, we can assume that the value will never need to be dynamic. To accommodate this, parameters can be declared static with the argument dynamic=False. This prevents Dymos from "fanning out" the static value to the n nodes in the ODE system.

Providing the ODE to the Phase

Phases in Dymos are instantiated with both the ode_class and the transcription to be used. Internally, Dymos needs to instantiate the ODE multiple times. This instantiation takes the form:

ode_instance = ode_class(num_nodes=<int>, **ode_init_kwargs)

This allows an OpenMDAO ExecComp to be used as an ODE via a lambda. For instance, the brachistochrone ODE can be written as:

ode = lambda num_nodes: om.ExecComp(['vdot = g * cos(theta)',
                                     'xdot = v * sin(theta)',
                                     'ydot = -v * cos(theta)'],
                                    g={'value': 9.80665, 'units': 'm/s**2'},
                                    v={'shape': (num_nodes,), 'units': 'm/s'},
                                    theta={'shape': (num_nodes,), 'units': 'rad'},
                                    vdot={'shape': (num_nodes,),
                                          'units': 'm/s**2',
                                          'tags': ['state_rate_source:v']},
                                    xdot={'shape': (num_nodes,),
                                          'units': 'm/s',
                                          'tags': ['state_rate_source:v']},
                                    ydot={'shape': (num_nodes,),
                                          'units': 'm/s',
                                          'tags': ['state_rate_source:v']},
                                    has_diag_partials=True)

phase = dm.Phase(ode_class=ode, transcription=t)

Note the use of has_diag_partials=True to provide more efficient graph coloring for the derivatives.

In theory, this also means you can implement Python's __call__ method for an ODE. The following code will return a copy of the brachistochrone ODE with the appropriate number of nodes. Note that the implementation below does not deal with any options provided via the ode_init_kwargs.

class CallableBrachistochroneODE(om.ExplicitComponent):

   def initialize(self):
       self.options.declare('num_nodes', types=int)

   def __call__(self, num_nodes, **kwargs):
       from copy import deepcopy
       ret = deepcopy(self)
       ret.options['num_nodes'] = num_nodes
       return ret

   def setup(self):
       nn = self.options['num_nodes']

       # Inputs
       self.add_input('v', val=np.zeros(nn), desc='velocity', units='m/s')

       self.add_input('g', val=9.80665, desc='grav. acceleration', units='m/s/s')

       self.add_input('theta', val=np.ones(nn), desc='angle of wire', units='rad')

       self.add_output('xdot', val=np.zeros(nn), desc='velocity component in x', units='m/s',
                       tags=['state_rate_source:x', 'state_units:m'])

       self.add_output('ydot', val=np.zeros(nn), desc='velocity component in y', units='m/s',
                       tags=['state_rate_source:y', 'state_units:m'])

       self.add_output('vdot', val=np.zeros(nn), desc='acceleration magnitude', units='m/s**2',
                       tags=['state_rate_source:v', 'state_units:m/s'])

       self.declare_partials(of='*', wrt='*', method='cs')
       self.declare_coloring(wrt='*', tol=1.0E-12)

An instance of the above ODE can then be provided to the phase upon instantiation.

ode = CallableBrachistochroneODE(num_nodes=1)
phase = dm.Phase(ode_class=ode, transcription=t)

This can potentially lead to unintended behavior if multiple copeis of the ODE are intended to share data. See the Python docs for some of the limitations of deepcopy.