Dymos: как я могу настроить объект для минимизации значения функции переменных состояния в конце времени моделирования

Используя стандартную задачу о брахистохроне, я сделал игрушечную модификацию целевой функции, чтобы показать вам, как добиться того, о чем вы спрашиваете. Новая целевая функция, которую я определил, не является сверхзначимой, но она передает суть.

Ключевая идея здесь заключается в том, что вам разрешено добавлять в модель дополнительные компоненты в дополнение к группе traj. Они могут быть до или после группы traj, но в вашем случае они нужны после. Затем вы можете подключать вещи к другим компонентам и группе traj.

Здесь я сделал очень простой ExecComp, который суммировал начальные и конечные условия переменной состояния V и установил это в качестве цели. Для этого я подключил к timeseries выходам фазы. Очевидно, вы могли бы построить другое уравнение, которое работает с окончательным значением состояния времени. Операторы подключения в этом примере показывают, как это сделать.

import numpy as np
import openmdao.api as om
import dymos as dm
import matplotlib.pyplot as plt

# First define a system which computes the equations of motion
class BrachistochroneEOM(om.ExplicitComponent):
    def initialize(self):
        self.options.declare('num_nodes', types=int)

    def setup(self):
        nn = self.options['num_nodes']

        # Inputs
        self.add_input('v', val=np.zeros(nn), units='m/s', desc='velocity')
        self.add_input('theta', val=np.zeros(nn), units='rad', desc='angle of wire')
        self.add_output('xdot', val=np.zeros(nn), units='m/s', desc='x rate of change')
        self.add_output('ydot', val=np.zeros(nn), units='m/s', desc='y rate of change')
        self.add_output('vdot', val=np.zeros(nn), units='m/s**2', desc='v rate of change')

        # Ask OpenMDAO to compute the partial derivatives using complex-step
        # with a partial coloring algorithm for improved performance
        self.declare_partials(of='*', wrt='*', method='cs')
        self.declare_coloring(wrt='*', method='cs', show_summary=True)

    def compute(self, inputs, outputs):
        v, theta = inputs.values()
        outputs['vdot'] = 9.80665 * np.cos(theta)
        outputs['xdot'] = v * np.sin(theta)
        outputs['ydot'] = -v * np.cos(theta)

p = om.Problem()

# Define a Trajectory object
traj = p.model.add_subsystem('traj', dm.Trajectory())

# Define a Dymos Phase object with GaussLobatto Transcription
tx = dm.GaussLobatto(num_segments=10, order=3)
phase = dm.Phase(ode_class=BrachistochroneEOM, transcription=tx)

traj.add_phase(name='phase0', phase=phase)

# Set the time options
phase.set_time_options(fix_initial=True,
                       duration_bounds=(0.5, 10.0))
# Set the state options
phase.add_state('x', rate_source='xdot',
                fix_initial=True, fix_final=True)
phase.add_state('y', rate_source='ydot',
                fix_initial=True, fix_final=True)
phase.add_state('v', rate_source='vdot',
                fix_initial=True, fix_final=False)
# Define theta as a control.
phase.add_control(name='theta', units='rad',
                  lower=0, upper=np.pi)


######################################################
# Post trajectory calculations for composite objective
######################################################

p.model.add_subsystem('other_calcs', om.ExecComp("f=v[0]+v[1]", v = {"shape":(2,)}))

p.model.connect('traj.phase0.timeseries.states:v', 'other_calcs.v', src_indices=[0,-1])

###################################################################
# Standard Dymos API for objective function - minimize final time.
###################################################################

# phase.add_objective('time', loc='final')

###################################################################
# OpenMDAO API for objective functions that require post-dymos calculations
###################################################################

p.model.add_objective('other_calcs.f')

# Set the driver.
p.driver = om.ScipyOptimizeDriver()

# Allow OpenMDAO to automatically determine total
# derivative sparsity pattern.
# This works in conjunction with partial derivative
# coloring to give a large speedup
p.driver.declare_coloring()

# Setup the problem
p.setup()

# Now that the OpenMDAO problem is setup, we can guess the
# values of time, states, and controls.
p.set_val('traj.phase0.t_duration', 2.0)

# States and controls here use a linearly interpolated
# initial guess along the trajectory.
p.set_val('traj.phase0.states:x',
          phase.interp('x', ys=[0, 10]),
          units='m')
p.set_val('traj.phase0.states:y',
          phase.interp('y', ys=[10, 5]),
          units='m')
p.set_val('traj.phase0.states:v',
          phase.interp('v', ys=[0, 5]),
          units='m/s')
# constant initial guess for control
p.set_val('traj.phase0.controls:theta', 90, units='deg')

# Run the driver to solve the problem and generate default plots of
# state and control values vs time
dm.run_problem(p, make_plots=True, simulate=True)

(1) всегда ли индекс [-1] переменных timeseries является последним моментом времени в транскрипции?

Да, это так.

(2) Могу ли я добавить к узлам один или несколько конкретных моментов времени? Возможно ли, чтобы узлы гарантированно находились в определенных фиксированных точках времени, которые могут не совпадать с точками коллокации?

Да...

Есть «правильный» способ сделать это, и вы обязательно должны попробовать его в первую очередь. Тогда есть «неправильный» способ, который я упомяну для полноты, и он может работать для вас, но также, вероятно, имеет некоторые числовые проблемы, о которых я не буду вдаваться в подробности, но это причина того, что это не правильный путь.

--- Правильный путь ---

Если у вас есть какие-то промежуточные моменты времени, когда вам нужно что-то измерить/ограничить, вам следует разбить траекторию на несколько фаз, и вы можете установить точки останова в ключевые моменты времени, которые вам необходимо учитывать. Итак, если у вас есть траектория, продолжительность которой неизвестна, но она будет составлять не менее 10 секунд, и вы знаете, что вам нужно выполнить некоторую постобработку производительности при t = 10, тогда сделайте ее двухфазной. Одна фаза от t=0 до t=10 и вторая фаза от t=10 до t=t_final. Даже если t_final равно 10, все равно все в порядке. Dymos позволяет использовать 0 фаз продолжительности.

--- Другой способ ---

Если поставить свою задачу с фиксированной продолжительностью, то это точно возможно. Возможно, вам повезет, и точки коллокации совпадут с нужным вам временем. Если вам не повезет, вы можете использовать другую технику, чтобы заставить Dymos интерполировать временную историю во вторичную сетку для вас. Эта сетка может иметь другую транскрипцию и количество сегментов, чем основная сетка словосочетаний (но будет иметь то же время начала и окончания). Вы можете выбрать транскрипцию Radau с порядком 1, что дает вам простой фиксированный временной шаг. Тогда количество сегментов в основном равно количеству временных шагов, и теоретически вы можете гарантировать, что у вас есть интерполированное значение в любое время, которое вы можете достичь с помощью целочисленного деления вашей общей продолжительности. Если вам нужно два раза, которые не совпадают с одним целочисленным делением, вы можете добавить два дополнительных временных ряда с разным количеством сегментов.

Я упоминаю этот метод только потому, что есть некоторые проблемы, когда вы не просто хотите провести измерение в один момент времени, но хотите, чтобы оно было в некотором равномерном распределении времени (возможно, вы хотите выполнить БПФ данных). .). В этих случаях вторичные временные ряды являются хорошим вариантом.

Если желаемое время не соответствует целочисленному делению, то определенно не используйте этот подход. Кроме того, если ваша фаза имеет нефиксированную продолжительность, то ДЕЙСТВИТЕЛЬНО определенно не используйте этот подход. Это просто не сработает, поскольку вы не можете удерживать фиксированные точки времени при изменении продолжительности.

(3.1) Можете ли вы иметь несколько целей?

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

Только одна цель. Вызовите ЛИБО API-интерфейс dymos или OpenMDAO. Не оба.

(3.2), Как можно настроить постобработку с интегрированным значением?

Dymos по своей сути является интегратором. Если вам нужно интегрированное значение, вы просто добавляете его как новое состояние к фазе. Вам вообще не нужно редактировать ODE. Dymos понимает, что вы пытаетесь сделать, если вы передаете существующее состояние в качестве источника скорости для нового состояния. Таким образом, для проблемы брахистохроны, если вы хотите отслеживать общее интегрированное движение по оси y, вы должны добавить новое состояние (мы назовем его y_tot, а state_rate будет самим y). Тогда последний момент времени во временном ряду для y_tot будет иметь интегрированное значение.

Вот сценарий, демонстрирующий приемы, которые я описал для (2) и (3). Добавление дополнительных фаз с точками останова в определенное время описано в этом примере документации Dymos.

import numpy as np
import openmdao.api as om
import dymos as dm
import matplotlib.pyplot as plt

# First define a system which computes the equations of motion
class BrachistochroneEOM(om.ExplicitComponent):
    def initialize(self):
        self.options.declare('num_nodes', types=int)

    def setup(self):
        nn = self.options['num_nodes']

        # Inputs
        self.add_input('v', val=np.zeros(nn), units='m/s', desc='velocity')
        self.add_input('theta', val=np.zeros(nn), units='rad', desc='angle of wire')
        self.add_output('xdot', val=np.zeros(nn), units='m/s', desc='x rate of change')
        self.add_output('ydot', val=np.zeros(nn), units='m/s', desc='y rate of change')
        self.add_output('vdot', val=np.zeros(nn), units='m/s**2', desc='v rate of change')

        # Ask OpenMDAO to compute the partial derivatives using complex-step
        # with a partial coloring algorithm for improved performance
        self.declare_partials(of='*', wrt='*', method='cs')
        self.declare_coloring(wrt='*', method='cs', show_summary=True)

    def compute(self, inputs, outputs):
        v, theta = inputs.values()
        outputs['vdot'] = 9.80665 * np.cos(theta)
        outputs['xdot'] = v * np.sin(theta)
        outputs['ydot'] = -v * np.cos(theta)

p = om.Problem()

# Define a Trajectory object
traj = p.model.add_subsystem('traj', dm.Trajectory())

# Define a Dymos Phase object with GaussLobatto Transcription
tx = dm.GaussLobatto(num_segments=10, order=3)
phase = dm.Phase(ode_class=BrachistochroneEOM, transcription=tx)

traj.add_phase(name='phase0', phase=phase)

# Set the time options
phase.set_time_options(fix_initial=True,
                       duration_bounds=(0.5, 10.0))
# Set the state options
phase.add_state('x', rate_source='xdot',
                fix_initial=True, fix_final=True)
phase.add_state('y', rate_source='ydot',
                fix_initial=True, fix_final=True)
phase.add_state('y_tot', rate_source='y',
                fix_initial=True, fix_final=True)
phase.add_state('v', rate_source='vdot',
                fix_initial=True, fix_final=False)
# Define theta as a control.
phase.add_control(name='theta', units='rad',
                  lower=0, upper=np.pi)


# compressed transcription so we don't get duplicate nodes between each segment
phase.add_timeseries('fixed_step', dm.Radau(num_segments=10, order=1))
phase.add_timeseries_output('v', timeseries = "fixed_step")

######################################################
# Post trajectory calculations for composite objective
######################################################

p.model.add_subsystem('other_calcs', om.ExecComp("f=v[0]+v[1]", v = {"shape":(2,)}))
p.model.connect('traj.phase0.timeseries.states:v', 'other_calcs.v', src_indices=[0,-1])

###################################################################
# Standard Dymos API for objective function - minimize final time.
###################################################################

phase.add_objective('time', loc='final')

###################################################################
# OpenMDAO API for objective functions that require post-dymos calculations
###################################################################

# p.model.add_objective('other_calcs.f')

# Set the driver.
p.driver = om.ScipyOptimizeDriver()

# Allow OpenMDAO to automatically determine total
# derivative sparsity pattern.
# This works in conjunction with partial derivative
# coloring to give a large speedup
p.driver.declare_coloring()

# Setup the problem
p.setup()

# Now that the OpenMDAO problem is setup, we can guess the
# values of time, states, and controls.
p.set_val('traj.phase0.t_duration', 2.0)

# States and controls here use a linearly interpolated
# initial guess along the trajectory.
p.set_val('traj.phase0.states:x',
          phase.interp('x', ys=[0, 10]),
          units='m')
p.set_val('traj.phase0.states:y',
          phase.interp('y', ys=[10, 5]),
          units='m')
p.set_val('traj.phase0.states:y_tot',
          phase.interp('y', ys=[0, 0]),
          units='m*s')
p.set_val('traj.phase0.states:v',
          phase.interp('v', ys=[0, 5]),
          units='m/s')
# constant initial guess for control
p.set_val('traj.phase0.controls:theta', 90, units='deg')

# Run the driver to solve the problem and generate default plots of
# state and control values vs time
dm.run_problem(p, make_plots=True, simulate=True)


print(p.get_val("traj.phase0.fixed_step.states:v"))
time = p.get_val("traj.phase0.fixed_step.time")[:,0]
dt = time[1:] - time[0:-1]
print("time", time)
# note that dymos keeps duplicates of values between each segment, so you see double values
# which gives 0 time steps between each segment
print("dt: ", dt)
print(p.get_val("traj.phase0.timeseries.states:y_tot"))