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

Я хочу использовать Dymos для решения задачи оптимального управления:

введите здесь описание изображения

с учетом динамической системы:

введите здесь описание изображения

У меня есть два вопроса:

(1) Как установить объектную функцию (object1) V (X (T)), где значения переменных состояния в конечный момент времени являются входными данными функции V.

Я думаю, что настройка

phase.add_objective('V(X)', loc='final', scale=1)

минимизирует интеграл от t=0 до t=T переменной, скорость которой равна V(X), верно?

(2) Я думаю, что могу настроить выход objectval2 со скоростью изменения времени L (X, u), а затем установить

phase.add_objective('objectval2', loc='final', scale=1)

тогда как мне сообщить оптимизатору, что я хочу минимизировать сумму objectval1 + objectval2 вместо objectval1*objectval2 или любой другой функции objectval1 и objectval2.

Я прочитал учебник и не нашел примера, похожего на мой случай.

#################

Комментарии: Спасибо. Я не могу добавить следующие уточняющие вопросы в комментариях к ответу, он сказал, что они слишком длинные. Поэтому я добавил их сюда.

Я думаю, что ключевой код здесь:

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

Я видел это в учебнике, и у меня есть 3 вопроса:

(1) Для ввода ключевого слова src_indices я думаю, что src_indices=[-1] передает последнее значение v(v[-1]) в a, а src_indices=[0] передает первое значение v(v[0]) в b. Я думаю, что момент времени в v[-1] зависит от количества узлов, которое определяется параметрами num_segments и order. Так всегда ли функция транскрипции устанавливает позицию последнего узла в последний момент времени? Например, если я решаю задачу управления от t=0 до t=2000, независимо от того, как я устанавливаю num_segments и order, т.е. num_segments=2 и order=3 или num_segments=20 и order=5, v[-1] всегда будет значением v в t=2000, а v[0] всегда будет значением v в t=0?

(2) Могу ли я добавить к узлам один или несколько конкретных моментов времени? Например, если я решаю задачу управления от t = 0 до t = 2000, я хочу минимизировать сумму значений v в t=500, t=1000 и t=2000, могу ли я потребовать, чтобы узлы включали моменты времени в t=500, t=1000 и t=2000, независимо от того, как я установил num_segments и order?

(3) Последний вопрос касается настройки целевой функции, которую я задал в первый раз. В примере вы закомментировали # phase.add_objective('time', loc='final'). Что сделает оптимизатор, если я установлю код:

phase.add_objective('x', loc='final')
p.model.add_objective('other_calcs.f')

Минимизирует ли оптимизатор сумму интеграла x от t=0 до t=T и v(T)? Если нет, то как настроить оптимизатор на минимизацию суммы интеграла x от t=0 до t=T и v(T)?

Пример Джастина довольно подробный, и здесь есть описание того, как это сделать: openmdao.github.io/dymos/faq/downstream_analysis.html

Rob Falck 19.11.2022 23:48

Спасибо. Я видел в учебнике пример с использованием om.ExecComp и src_indices, но мне все еще не ясно, и я добавил свои уточняющие вопросы к исходным вопросам.

Wei 20.11.2022 01:45
Шаблоны Angular PrimeNg
Шаблоны Angular PrimeNg
Как привнести проверку типов в наши шаблоны Angular, использующие компоненты библиотеки PrimeNg, и настроить их отображение с помощью встроенной...
Создайте ползком, похожим на звездные войны, с помощью CSS и Javascript
Создайте ползком, похожим на звездные войны, с помощью CSS и Javascript
Если вы веб-разработчик (или хотите им стать), то вы наверняка гик и вам нравятся "Звездные войны". А как бы вы хотели, чтобы фоном для вашего...
Документирование API с помощью Swagger на Springboot
Документирование API с помощью Swagger на Springboot
В предыдущей статье мы уже узнали, как создать Rest API с помощью Springboot и MySql .
Начала с розового дизайна
Начала с розового дизайна
Pink Design - это система дизайна Appwrite с открытым исходным кодом для создания последовательных и многократно используемых пользовательских...
Шлюз в PHP
Шлюз в PHP
API-шлюз (AG) - это сервер, который действует как единая точка входа для набора микросервисов.
14 Задание: Типы данных и структуры данных Python для DevOps
14 Задание: Типы данных и структуры данных Python для DevOps
проверить тип данных используемой переменной, мы можем просто написать: your_variable=100
2
2
99
2
Перейти к ответу Данный вопрос помечен как решенный

Ответы 2

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

Ключевая идея здесь заключается в том, что вам разрешено добавлять в модель дополнительные компоненты в дополнение к группе 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)

Спасибо. Я видел в учебнике пример с использованием om.ExecComp и src_indices, но мне все еще не ясно, и я добавил свои уточняющие вопросы к исходным вопросам.

Wei 20.11.2022 01:45
Ответ принят как подходящий

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

Большое спасибо. Я установил продолжительность t от 0 до 10 и использовал dm.Radau(num_segments=10, order=1), я видел на каждой итерации «время = [ 0. 1. 1. 2. 2. 3. 3. 4. 4. 5. 5. 6. 6. 7. 7. 8. 8. 9. 9. 10.]' с дублированием моментов времени между каждым сегментом в подсистеме. Я могу настроить num_segments, чтобы получить нужные мне моменты времени, или интерполировать их в подсистеме. Таким образом, после того, как я установил num_segments и order на «транскрипцию», временные точки не меняются алгоритмом на каждой итерации? Или только временные точки в phase.add_timeseries не изменятся?

Wei 20.11.2022 22:00

Использование order=1 в транскрипциях может привести к неожиданному поведению. Мое предложение состояло бы в том, чтобы использовать вторичный временной ряд, который использует транскрипцию GaussLobatto порядка 3. Если мы установим временной ряд для вывода в узлах ввода состояния, это должно фиксировать только моменты времени окончания сегмента. См. этот пример для настройки вторичных выходных данных временных рядов. После tx1 = dm.GaussLobatto(num_segments=10, order=3, compressed=False) добавьте временной ряд как phase0.add_timeseries('timeseries2', transcription=tx1, subset='state_input')

Rob Falck 20.11.2022 22:58

Спасибо. Я не совсем понимаю ваше высказывание: «Если мы установим временной ряд для вывода в узлах ввода состояния, это должно фиксировать только моменты времени окончания сегмента». Сначала установите tx1 = dm.GaussLobatto (num_segments = 10, order = 3). Когда я использую «phase0.add_timeseries («timeseries2», транскрипция = tx1)», точки времени, вводимые в подсистему, равны [0. , 0,5, 1, 1, 1,5, 2, 2, 2,5, 3, 3, 3,5, 4, 4, 4,5, 5, 5, 5,5, 6, 6, 6,5 , 7. , 7. , 7.5, 8. , 8. , 8.5, 9. , 9. , 9.5, 10. ].

Wei 21.11.2022 03:49

Когда я использую 'phase0.add_timeseries('timeseries2', транскрипция = tx1, subset = 'state_input')', точки времени, вводимые в подсистему, равны [ 0., 1., 1., 2., 2., 3., 3., 4., 4., 5., 5., 6., 6., 7., 7., 8., 8., 9., 9., 10.]. На 10 моментов времени меньше по сравнению с тем, чтобы не использовать subset='state_input'. Что происходит, когда я устанавливаю аргументы ключевых слов subset='state_input? Существует ли уравнение для расчета разницы количества введенных моментов времени между установкой и отсутствием установки subset='state_input из параметра num_segments и order?

Wei 21.11.2022 03:58

Для определения многочлена 3-го порядка требуется 4 элемента информации. В GaussLobatto 3-го порядка это значения состояния и скорости состояния в конечных точках каждого сегмента. Здесь вводятся состояния (узлы ввода состояния). Ограничения «дефекта» оцениваются в средних точках сегмента. Это узлы, которые не отображаются, когда вы указываете, что вам нужны узлы ввода состояния. Здесь немного больше описания: openmdao.github.io/dymos/getting_started/collocation.html

Rob Falck 21.11.2022 20:46

Понятно. Спасибо. У меня простой вопрос по алгоритму. После решения проблемы есть ли информация от dm или p, чтобы проверить, успешно ли завершилась оптимизация? И сообщить финал collocation defect? Иногда решение и симуляция хорошо подходят (я думаю, что collocation defect в этом случае малы), а иногда нет. Но они печатают одну и ту же информацию EXIT: Maximum Number of Iterations Exceeded или EXIT: Restoration Failed!. Даже я поставил p.driver.opt_settings['max_iter']=int(1e5), там тоже написано EXIT:Maximum Number of Iterations Exceeded.

Wei 22.11.2022 01:21

Похоже, вы используете SLSQP. Я настоятельно рекомендую вместо этого переключиться на IPOPT. Вы можете получить его как часть pyoptsparse (помощник по установке для Linux: github.com/OpenMDAO/build_pyoptsparse)

Justin Gray 22.11.2022 15:42

Я использую IPOPT. У меня работает только IPOPT. Другие оптимизаторы дают мне ошибку. Моя настройка optimizer = 'IPOPT', p.driver.options['optimizer'] = optimizer, p.driver.opt_settings['max_iter'] = int(5e4). Есть ли информация в dm или p, которая говорит мне об успешности окончательной оптимизации или нет, или какая-то окончательная оценка collocation defect. Он печатает EXIT: Maximum Number of Iterations Exceeded или ВЫХОД: Ошибка восстановления!. When I plot, sometimes the collocations are good sometimes are bad. But the printed info are the same, ВЫХОД: Превышено максимальное количество итераций`.

Wei 22.11.2022 18:53

Другие вопросы по теме