How to integrate discrete time estimation into OpenMDAO Components?

Question

I am trying to familiarise myself with OpenMDAO. One thing I am finding hard to understand is how integrated/state variables work in separate components of OpenMDAO. I think it's quite crucial and I want to understand the basics.

So let's say I have a rocket with constant thrust and changing mass and I want to model the first 10 seconds of flight with discrete time and a timestep of 0.1 second. Let's say the changing mass is a function of time as well, but also pretend it is a difficult variable dependent on many things, so we want to have it calculated in a different component and simply be an input to this component.

-How does one ensure that the Mass is updated in the discrete timestep computation of the Velocity, while being in a different component?

The code below is a mediocre attempt at what I would try to solve this example problem, but I am rather sure now Mass is a simple static input, while it should change in some arbitrary relation to time. (Let's assume the Thrust is constant)

from openmdao.api import ExplicitComponent

class VelocityComponent(ExplicitComponent):

    def setup(self):
        self.add_input('T', desc='Propulsion Thrust')
        self.add_input('M', desc='Instanteneous Mass')\

        self.add_output('v', desc='Satellite Velocity')

        self.declare_partials('*','*')

    def compute(self, inputs, outputs)

        v = 10                          #some initial velocity value
        t = 0                           #initial time value
        tstep = 0.1
        tend = 10

        for i in range(0,tend/tstep):
            a = inputs['T']/inputs['M'] #calculate acceleration  
            v += a                      #update velocity
            t += tstep                  #next time step 

        outputs['v'] = v 

velocity, v, should be integrated with time dependent acceleration, a, not a constant acceleration.

PS:As I am rather new to all of this, but willing to learn, any tips for resources that could help a beginner like me with OpenMDAO are much appreciated.

PSPS: I have read the Beginner and Advanced User Guide of the OpenMDAO documentation, but couldn't find an example with integrated variables. The old documentation has an example of an engine and transmission system, and its engine component does include state variables and some discrete integration steps, but it uses an older OpenMDAO version and I don't know how that would work in the newer version (if I even understand the old one correctly)

Answer 1

This is a discussion about ODE integration, and you've picked a pretty complex topic because 1) there are lots of different ways to accomplish the integration (e.g. explicit Euler, RK4, BDF ... ) 2) carrying analytic derivatives through the time-integration is extra-tricky.

for #2, the root of the difficulty is exactly the problem you have identified. You can not use a simple for-loop structure inside a single component, if you also want to construct your ODE out of a set of different components organized into a group.

The good news is that there is a library already written that handles all the time-integration for you: Dymos. As of April 2019, this library is under active development by the OpenMDAO team itself and is still undergoing some modes API revisions and feature additions. Despite the somewhat fluid APIs, I would recommend you take a look at the examples there. For complex problems, this is your best bet.

However, you can achieve simple time-stepping euler integration without an additional library. The trick is that you stamp out one instance of your ODE for each time-step and pass the new state down the line from one instance to the next. Here is a dead simple example of a falling object under constant gravity:

from openmdao.api import IndepVarComp, Problem, ExplicitComponent


class Cannonball(ExplicitComponent): 

    def initialize(self): 

        self.options.declare('delta_t', default=0.1)

    def setup(self): 

        self.add_input('Yi', units='m',   desc='position at the start of the time-step')
        self.add_input('Vi', units='m/s', desc='velocity at the start of the time-step')

        self.add_output('Ye', units='m',   desc='position at the end of the time-step')
        self.add_output('Ve', units='m/s', desc='velocity at the end of the time-step')

        self.declare_partials(of='*', wrt='*', method='cs')


    def compute(self, inputs, outputs): 

        dt = self.options['delta_t']

        outputs['Ve'] = 9.81 * dt + inputs['Vi']
        outputs['Ye'] = 0.5 * 9.81 * dt**2 + inputs['Vi'] * dt + inputs['Yi']


if __name__ == "__main__": 
    import numpy as np
    import matplotlib.pylab as plt

    N_TIMES = 10

    p = Problem()

    ivc = p.model.add_subsystem('init_conditions', IndepVarComp(), promotes=['*'])
    ivc.add_output('Y0', 100., units='m')
    ivc.add_output('V0', 0, units='m/s')

    p.model.connect('Y0', 't_0.Yi')
    p.model.connect('V0', 't_0.Vi')

    for i in range(N_TIMES): 
        p.model.add_subsystem(f't_{i}', Cannonball())


    for i in range(N_TIMES-1): 
        p.model.connect(f't_{i}.Ye', f't_{i+1}.Yi')
        p.model.connect(f't_{i}.Ve', f't_{i+1}.Vi')

    p.setup()

    p.run_model()


    # collect the data into an array for plotting
    Y = [p['Y0'],]
    V = [p['V0'],]
    for i in range(N_TIMES): 
        Y.append(p[f't_{i}.Ye'])
        V.append(p[f't_{i}.Ve'])

    times = np.arange(N_TIMES+1) * .01 # delta_t

    fig, ax = plt.subplots()
    ax.plot(times, Y)
    ax.set_ylabel('velocity (m/s')
    ax.set_xlabel('time (s)')

    plt.show()

This will give a feed-foward in time structure to your model (generated using OpenMDAO's built-in N2 viewer).

And you can see that you get the expected quadratic position with respect to time.

You could make a more complex integration by coding a different ODE scheme into this component (RK4 for example). You could also write a more complex group, made of of multiple components, to serve as your time-step and then stamp out multiple copies of that group.

I want to stress that, while the above example is good for generating an understanding this is not a very efficient way to actually do time-integration in OpenMDAO. Dymos does things very differently internally, working with components that are vectorized in time for much greater efficiency. None the less, if you really are interested in the most dead simple time-integration scheme within OpenMDAO... this is it.