<jupyter_start><jupyter_text>Defining ODEsOptimal 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 OptionsODE 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 ProblemDymos 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. TimeThe following time options can be set via the `set_time_options` method of Phase:<jupyter_code>from IPython.display import display, HTML, IFrame, Code
from dymos.phase.options import TimeOptionsDictionary
obj = TimeOptionsDictionary()
display(HTML(obj.to_table(fmt='html')))<jupyter_output><empty_output><jupyter_text>StatesStates have the following options set via the `add_state` and `set_state_options` methods of Phase:<jupyter_code>from dymos.phase.options import StateOptionsDictionary
obj = StateOptionsDictionary()
display(HTML(obj.to_table(fmt='html')))<jupyter_output><empty_output><jupyter_text>ControlsInputs 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:<jupyter_code>from dymos.phase.options import ControlOptionsDictionary
obj = ControlOptionsDictionary()
display(HTML(obj.to_table(fmt='html')))<jupyter_output><empty_output><jupyter_text>ParametersInputs 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:<jupyter_code>from dymos.phase.options import ParameterOptionsDictionary
obj = ParameterOptionsDictionary()
display(HTML(obj.to_table(fmt='html')))<jupyter_output><empty_output><jupyter_text>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 VariablesDymos will also automatically find and add any states that have been declared in components in the ODE. The syntaxfor declaring them is as follows. ```pythonself.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.| 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}'.                                      || 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.   | Automatic target detectionDymos will attempt to automatically connect variables to targets in the same name if those targetsexist at the top-level of the ODE and have the same name as the state, control, parameter, or 'time'. Vectorizing the ODEIn addition to specifying the ODE Options, a system used as an ODE is required to have a metadataentry 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 variablevalue at a discrete time in the trajectory.The nodes are discretization or collocation locations in the polynomials which representeach segment.  The number of nodes in a given phase (to be evaluated by the ODE system) is determinedby the number of segments in the phase and the polynomial order in each segment.  When Dymos instantiatesthe 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 somethinglike this:<jupyter_code>class MyODESystem(om.ExplicitComponent):

    def initialize(self):
        self.metadata.declare('num_nodes', types=int)<jupyter_output><empty_output><jupyter_text>The inputs and outputs of the system are expected to provide a scalar or dimensionedvalue *at each node*.  Vectorization of the component via numpy adds a significant performance increasecompared to using a for loop to cycle through calculations at each node.  It's important to rememberthat 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 valuesof the input variables at time `t`.  When the output variables are scalar at any given time, thisresults in components whose Jacobian matrices are diagonal.  This large degree of sparsity leadsto computational advantages when using sparse-aware optimizers like SNOPT.  Users should declarethe 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 thesetup, compute, and compute partials methods might look like this:<jupyter_code>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<jupyter_output><empty_output><jupyter_text>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 intime (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.  Sincethe 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, inrow-major order.In this example, the partial of `y` with respect to `b` is linear, so we can simply provide it inthe `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 OutputsThe above example assumes all inputs and outputs are scalar at each node.  Sometimes the user mayencounter a situation in which the inputs and/or outputs are vectors, matrices, or tensors ateach node.  In this case the dimension of the variable is `num_nodes`, with the dimension of thevariable 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 InputsDeclaring 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 PhasePhases 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: ```pythonode_instance = ode_class(num_nodes=, **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:<jupyter_code>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)<jupyter_output><empty_output><jupyter_text>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`.<jupyter_code>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)<jupyter_output><empty_output><jupyter_text>An instance of the above ODE can then be provided to the phase upon instantiation.<jupyter_code>ode = CallableBrachistochroneODE(num_nodes=1)
t = dm.Radau(num_segments=10, order=3)
phase = dm.Phase(ode_class=ode, transcription=t)<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.09673820305834482, 'fraction_numerical': 0.0024939292511649276, 'mean_word_length': 3.242204899777283, 'pattern_counts': {'":': 0, '<': 37, '<?xml version=': 0, '>': 37, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 38, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '141628', 'license': 'Apache-2.0', 'n_tokens_mistral': 4010, 'n_tokens_neox': 3880, 'n_words': 1707, 'path': 'docs/getting_started/defining_odes.ipynb', 'repo_name': 'cashmesh/dymos', 'size': '15237'}