Defining ODEs

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.

(4)#\[\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 to 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 systems 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:

Option	Default	Acceptable Values	Acceptable Types	Description
duration_adder	N/A	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	N/A	N/A	['Number']	Unit-reference value for the duration of the integration variable across the phase.
duration_ref0	N/A	N/A	['Number']	Zero-reference value for the duration of the integration variable across the phase.
duration_scaler	N/A	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	N/A	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	N/A	N/A	['Number']	Unit-reference value for the initial value of the integration variable.
initial_ref0	N/A	N/A	['Number']	Zero-reference value for the initial value of the integration variable.
initial_scaler	N/A	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.
name	time	N/A	['str']	Name of the integraiton variable in the phase.
t_duration_balance_options	{}	N/A	N/A	options dictionary for the duration residual
t_duration_targets	[]	N/A	N/A	targets in the ODE to which the total duration of the phase is connected
t_initial_targets	[]	N/A	N/A	targets in the ODE to which the initial time of the phase is connected
targets	unspecified	N/A	N/A	targets in the ODE to which the integration variable is connected
time_phase_targets	unspecified	N/A	N/A	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:

Option	Default	Acceptable Values	Acceptable Types	Description	Deprecation
adder	N/A	N/A	['Iterable', 'Number']	Adder of the state variable at the discretization nodes. This option is invalid if opt=False.	N/A
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.	State option 'connected_initial' is deprecated. Please use 'input_initial' instead.
continuity	True	N/A	['bool', 'dict']	Enforce continuity of state values at segment boundaries. This option is invalid if opt=False.	N/A
continuity_ref	N/A	N/A	['Number']	Reference unit value for continuity at segment boundaries instead of scaler.This option is invalid if opt=False.	N/A
continuity_scaler	N/A	N/A	['Number']	Scaler for continuity at segment boundaries. This option is invalid if opt=False.	N/A
defect_ref	N/A	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.	N/A
defect_scaler	N/A	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.	N/A
desc		N/A	['str']	description of the state variable	N/A
final_bounds	N/A	N/A	['Iterable']	Bounds on the value of the state at the end of the phase. This option is invalid if opt=False.	N/A
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.	N/A
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.	N/A
initial_bounds	N/A	N/A	['Iterable']	Bounds on the value of the state at the start of the phase. This option is invalid if opt=False.	N/A
input_final	False	[True, False]	['bool']	Whether the final value of the state is expected to be connected to an exterior value.	N/A
input_initial	False	[True, False]	['bool']	Whether the initial value of the state is expected to be connected to an exterior value.	N/A
lower	N/A	N/A	['Iterable', 'Number']	Lower bound of the state variable at the discretization nodes. This option is invalid if opt=False.	N/A
name	Required	N/A	['str']	name of ODE state variable	N/A
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.	N/A
rate_source	N/A	N/A	['str']	ODE-path or phase variable providing the derivative of the state variable	N/A
ref	N/A	N/A	['Iterable', 'Number']	Unit-reference value of the state variable at the discretization nodes. This option is invalid if opt=False.	N/A
ref0	N/A	N/A	['Iterable', 'Number']	Zero-reference value of the state variable at the discretization nodes. This option is invalid if opt=False.	N/A
scaler	N/A	N/A	['Iterable', 'Number']	Scaler of the state variable at the discretization nodes. This option is invalid if opt=False.	N/A
shape	N/A	N/A	['Iterable']	shape of the state variable, as determined by introspection	N/A
solve_segments	N/A	[False, 'forward', 'backward']	N/A	If '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.	N/A
source	N/A	N/A	['str']	RHS-path or phase variable providing value of the state variable, for Analytic transcription only!	N/A
targets	unspecified	N/A	N/A	Targets in the ODE to which the state is connected	N/A
units	unspecified	N/A	N/A	units in which the state variable is defined	N/A
upper	N/A	N/A	['Iterable', 'Number']	Upper bound of the state variable at the discretization nodes. This option is invalid if opt=False.	N/A
val	1.0	N/A	['Iterable', 'Number']	Default value of the state variable at the discretization nodes	N/A

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:

Option	Default	Acceptable Values	Acceptable Types	Description
adder	N/A	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.
continuity_ref	N/A	N/A	['Number']	Reference unit value for continuity at segment boundaries instead of scaler.This option is invalid if opt=False.
continuity_scaler	N/A	N/A	['Number']	Scaler for continuity at segment boundaries. This option is invalid if opt=False.
desc		N/A	['str']	The description of the control variable.
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	N/A	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_ref	N/A	N/A	['Number']	Reference unit value for rate2 continuity at segment boundaries instead of scaler.This option is invalid if opt=False.
rate2_continuity_scaler	N/A	N/A	['Number']	Scaler of the rate2 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_ref	N/A	N/A	['Number']	Reference unit value for rate continuity at segment boundaries instead of scaler.This option is invalid if opt=False.
rate_continuity_scaler	N/A	N/A	['Number']	Scaler of the 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	N/A	N/A	['Iterable', 'Number']	The unit-reference value of the control variable at the nodes. This option is invalid if opt=False.
ref0	N/A	N/A	['Iterable', 'Number']	The zero-reference value of the control variable at the nodes. This option is invalid if opt=False.
scaler	N/A	N/A	['Iterable', 'Number']	The scaler of the control variable at the nodes. This option is invalid if opt=False.
shape	N/A	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	N/A	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:

Option	Default	Acceptable Values	Acceptable Types	Description	Deprecation
adder	N/A	N/A	['Iterable', 'Number']	The adder of the parameter. This option is invalid if opt=False.	N/A
desc		N/A	['str']	The description of the parameter.	N/A
dynamic	unspecified	[True, False, unspecified]	N/A	True if this parameter can be used as a dynamic control, else False.If _unspecified, attempt to determine through introspection.	Option dynamic has been replaced by option 'static_target' and will be removed in Dymos 2.0.0. Note that 'static_target' has the opposite meaning of option 'dynamic', so parameters with option 'dynamic' set to False should now use 'static_target' set to True.
include_timeseries	N/A	[True, False]	['bool']	True if the static parameters should be included in output timeseries, else False.If None (default) set the value based on Phase.timeseries_options["include_parameters"]	N/A
lower	N/A	N/A	['Iterable', 'Number']	The lower bound of the parameter. This option is invalid if opt=False.	N/A
name	Required	N/A	['str']	The name of ODE system parameter to be set via parameter.	N/A
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.	N/A
ref	N/A	N/A	['Iterable', 'Number']	The unit-reference value of the parameter. This option is invalid if opt=False.	N/A
ref0	N/A	N/A	['Iterable', 'Number']	The zero-reference value of the parameter. This option is invalid if opt=False.	N/A
scaler	N/A	N/A	['Iterable', 'Number']	The scaler of the parameter. This option is invalid if opt=False.	N/A
shape	unspecified	N/A	N/A	The shape of the parameter.	N/A
static_target	unspecified	N/A	N/A	True if the target of this parameter does NOT have a unique value at each node in the ODE.If _unspecified, attempt to determine through introspection.	Use option `static_targets` to specify whether all targets are static (static_targets=True), none are static (static_targets=False), static_targets are determined via introspection (static_targets=_unspecified), or give an explicit sequence of the static targets.
static_targets	unspecified	N/A	N/A	If a boolean, specifies whether all targets are static (True), or no targets are static (False). Otherwise, provide a list of the static targets within the ODE. If left unspecified, static targets will be determined by finding inptus tagged with 'dymos.static_target'.	N/A
targets	unspecified	N/A	N/A	Targets in the ODE to which the state is connected	N/A
units	unspecified	N/A	N/A	The units in which the parameter is defined.	N/A
upper	N/A	N/A	['Iterable', 'Number']	The upper bound of the parameter. This option is invalid if opt=False.	N/A
val	0.0	N/A	['Iterable', 'ndarray', 'Number']	The default value of the parameter in the phase.	N/A

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 not sized as num_nodes in the first dimension), then the user may specify option static_target = True to override this behavior.

Tagging Variables#

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=['dymos.state_rate_source:v', 'dymos.state_units:m/s'])

The state is defined by adding tags to the state rate’s output. The tag ‘dymos.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 ‘dymos.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.

Inputs which are sized as ‘static’ in the ODE (the first dimension is not the number of nodes), must be made known to Dymos when using them as parameters. This is normally accomplished with the ‘static_target=True’ option when calling add_parameter. To avoid repetitive use of this option, users can add the tag ‘dymos.static_target’ to any input which should always be considered static. When connecting to multiple targets, a ‘dymos.static_target’ tag on one of them will result in all of them being considered static. Inputs tagged with ‘dymos.static_target’ cannot be used as controls or polynomial controls.

Similarly, the tag ‘dymos.static_output’ is used to convey that an output of the ODE is not sized with num_nodes as its first dimension. These outputs may not be used as timeseries outputs, nor may they be used as boundary or path constraints. The general OpenMDAO add_constraint method can still be applied to them by using the full path to the output in the ODE.

Tag	Applies to	Description
dymos.state_rate_source:{state_name}	Outputs	Specifies that the tagged output provides the rate used to integrate the state of the given name.
dymos.state_units:{units}	Outputs	Specifies the default units to be used for the state whose rate is provided by the tagged output. May only be used on outputs also tagged with ‘dymos.state_rate_source:{state_name}’ or ‘dymos.state_source:{state_name}’.
dymos.state_source:{state_name}	Outputs	For AnalyticPhase only, this tag specifies that the output provides the analytic solution to a state.
dymos.static_target	Inputs	Specifies that the input is not shaped with num_nodes as the first dimension. This effects the way in which parameters are connected to the input. Inputs tagged with ‘dymos.static_target’ cannot be used as controls.
dymos.static_output	Outputs	Specifies that the output is not shaped with num_nodes as the first dimension, and may not be used as a timeseries output, boundary constraint, or path constraint.

Automatic target detection#

Dymos will attempt to automatically connect variables to targets in the same name if those targets exist at the top-level of the ODE and have the same name as the state, control, parameter, or ‘time’.

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:

class MyODESystem(om.ExplicitComponent):

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

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.

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 static_target=True. 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)

t = dm.Radau(num_segments=10, order=3)
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')

        # Outputs
        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'])

        # Derivatives
        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)
t = dm.Radau(num_segments=10, order=3)
phase = dm.Phase(ode_class=ode, transcription=t)

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