openmdao - 动力 : How can I set object to minimize a function value of state variables at the end of simulation time

Question

我想用Dymos解决最优控制问题:

enter image description here

服从动力系统:

enter image description here

我有两个问题:

(1)如何设置V(X(T))的目标函数(object1)，其中最终时间点的状态变量值是函数V的输入。

我觉得设置

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

将最小化速率为 V(X) 的变量从 t=0 到 t=T 的积分，对吗？

(2) 我想我可以设置一个时间变化率 L(X, u) 的输出 objectval2 然后设置

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 和 顺序，例如num_segments=2 和 order=3 或 num_segments=20 和 order=5, v[- 1] 将始终是 t=2000 时 v 的值，而 v[0] 将始终是 v 在 t=0?

(2) 我可以在节点上添加一个或多个特定的时间点吗？例如，如果我解决一个从 t=0 到 t=2000 的控制问题，我想最小化 t=500 时 v 的值之和，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')

优化器是否最小化从t=0到t=T和v(T)<的x积分之和？如果不是，我如何设置优化器以最小化从 t=0 到 t=T 和 v 的 x 积分总和(T)?

Answer 1

使用标准的 Brachistochrone 问题，我对目标函数进行了玩具修改，以向您展示如何实现您的要求。我定义的新目标函数不是一个非常有意义的目标函数，但它明白了这一点。

这里的关键思想是，除了 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)