Это обсуждение интеграции ODE, и вы выбрали довольно сложную тему, потому что 1) существует множество различных способов выполнить интеграцию (например, явный Эйлер, RK4, BDF ...) 2) перенос аналитических производныхчерез интеграцию времени очень сложно.
для # 2, корень сложности - именно та проблема, которую вы определили.Вы не можете использовать простую структуру for-loop внутри одного компонента, если вы также хотите создать свой ODE из набора различных компонентов, организованных в группу.
Хорошей новостью является то, что уже есть библиотека, которая обрабатывает все время интеграции для вас: Dymos .По состоянию на апрель 2019 года эта библиотека активно разрабатывается самой командой OpenMDAO и все еще подвергается пересмотру некоторых режимов API и добавлению функций.Несмотря на довольно плавные API, я бы порекомендовал вам взглянуть на примеры там.Для сложных задач это ваш лучший выбор.
Тем не менее, вы можете добиться простой интеграции euler с пошаговым временем без дополнительной библиотеки.Хитрость в том, что вы выделяете один экземпляр вашего ODE для каждого временного шага и передаете новое состояние по линии от одного экземпляра к следующему.Вот мертвый простой пример падающего объекта в условиях постоянной силы тяжести:
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. Dymos внутренне работает по-разному, работая с компонентами, которые векторизованы во времени для гораздо большей эффективности.Тем не менее, если вы действительно заинтересованы в самой мертвой простой схеме временной интеграции в OpenMDAO ... вот оно.