# 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*.

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_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. | 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 |

dymos.static_output |
Outputs |
Specifies that the output is |

## 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.