Как интегрировать оценку дискретного времени в компоненты OpenMDAO?

Это обсуждение интеграции ODE, и вы выбрали довольно сложную тему, потому что 1) существует множество различных способов выполнить интеграцию (например, явный Эйлер, RK4, BDF...) 2) проведение аналитических производных через интегрирование по времени очень сложно.

для № 2 корень трудности именно в той проблеме, которую вы определили. Вы не можете использовать простую структуру цикла for внутри одного компонента, если вы также хотите построить ОДУ из набора различных компонентов, организованных в группу.

Хорошей новостью является то, что уже написана библиотека, которая выполняет всю интеграцию по времени за вас: Димос. По состоянию на апрель 2019 года эта библиотека находится в стадии активной разработки самой командой OpenMDAO и все еще подвергается изменениям API некоторых режимов и добавлению функций. Несмотря на несколько плавные API, я бы порекомендовал вам взглянуть на приведенные там примеры. Для сложных задач это лучший выбор.

Однако вы можете выполнить простую интеграцию Эйлера с временным шагом без дополнительной библиотеки. Хитрость заключается в том, что вы штампуете один экземпляр вашего ОДУ для каждого временного шага и передаете новое состояние по цепочке от одного экземпляра к другому. Вот очень простой пример падающего объекта под действием постоянной гравитации:

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()

Это даст временную структуру для вашей модели (сгенерированную с использованием Встроенный в OpenMDAO просмотрщик N2).

И вы можете видеть, что вы получаете ожидаемую квадратичную позицию по времени.

Вы можете сделать более сложную интеграцию, закодировав другую схему ODE в этот компонент (например, RK4). Вы также можете написать более сложную группу, состоящую из нескольких компонентов, которая будет служить временным шагом, а затем отпечатать несколько копий этой группы.

Я хочу подчеркнуть, что, хотя приведенный выше пример хорош для понимания, это не очень эффективный способ выполнить интеграцию времени в OpenMDAO. Димос внутренне делает вещи совсем по-другому, работая с компонентами, которые векторизованы во времени для гораздо большей эффективности. Тем не менее, если вас действительно интересует самая простая схема временной интеграции в OpenMDAO... вот она.