Dymos/OpenMDAO相位联动解决读取错误问题

Question

我正在运行并研究此处提供的 Dymos/OpenMDAO 示例： SSTO 地球发射

它运行良好，具有不同的边界条件和参数。

我的目标是使轨迹分为两个阶段，这样我就可以模拟两级火箭，其中第一级将分离并产生质量“跳跃”。

我正在建立一个逐步模型。因此，首先是设置两个阶段，相同的边界条件（初始和最终）并使用以下方法链接它们：

traj.link_phases(phases=['phase0', 'phase1'], vars=['*'], connected=True)

它按预期工作，相同的轨迹，但有两个阶段。

我正在准备模型以减轻烧毁的第一级的重量。但现在只是将所有变量的联系改为显式（控制变量除外）：

traj.link_phases(phases=['phase0', 'phase1'], vars=['time','x','y','vx','vy','m'])

该案例有效，但在加载绘图解决方案时出现错误：

Exception has occurred: KeyError
'traj.linkages.phase0:m_final|phase1:m_initial'
KeyError: 'traj.linkages.phase0:m_final|phase1:m_initial'

During handling of the above exception, another exception occurred:

  File "/dymos/ssto/main_phases_works.py", line 113, in <module>
    sol = om.CaseReader('dymos_solution.db').get_case('final')
KeyError: 'traj.linkages.phase0:m_final|phase1:m_initial'

显然，它不是在读取链接信息进行绘图？有谁知道为什么吗？

谢谢你 豪尔赫

Answer 1

修改后的 SSTO 代码仅适用于一个阶段（仅用于正确比较）

import matplotlib.pyplot as plt
import openmdao.api as om
import dymos as dm
from dymos.examples.plotting import plot_results

p = om.Problem(model=om.Group())
p.driver = om.pyOptSparseDriver()

p.driver.options['optimizer'] = 'IPOPT'

p.driver.opt_settings['derivative_test'] = 'first-order'
p.driver.opt_settings['mu_strategy'] = 'adaptive'
p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'
p.driver.opt_settings['bound_mult_init_method'] = 'mu-based'
p.driver.opt_settings['mu_init'] = 0.01
p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'
p.driver.opt_settings['max_iter']= 1000
p.driver.opt_settings['tol']= 1e-2
p.driver.declare_coloring()

from launch_vehicle_ode import LaunchVehicleODE

traj = dm.Trajectory()

phase0 = dm.Phase(ode_class=LaunchVehicleODE, ode_init_kwargs={'CD': 1.0, 'S': 0.2},
                 transcription=dm.Radau(num_segments=10, order=3, compressed=False))

traj.add_phase('phase0', phase0)

p.model.add_subsystem('traj', traj)
#
phase0.set_time_options(    fix_initial=True, duration_bounds=(1, 500))
phase0.add_state('x',       fix_initial=True, rate_source='xdot', fix_final=True, lower = 0.0)
phase0.add_state('y',       fix_initial=True,rate_source='ydot',  fix_final=True, lower = 0.0)
phase0.add_state('vx',      fix_initial=True,rate_source='vxdot', fix_final=True, lower = 0.0)
phase0.add_state('vy',      fix_initial=True, rate_source='vydot',fix_final=True, lower = 0.0)
phase0.add_state('m',       fix_final=True,rate_source='mdot', lower = 500.0)
phase0.add_control('theta', fix_initial=True, units='deg', lower=-90.0, upper=90.0, targets=['theta'])
phase0.add_parameter('thrust', val = 25000.0, opt=False, targets=['thrust'])
#
phase0.add_objective('time', loc='final')

p.model.linear_solver = om.DirectSolver()
#
p.setup(check=True)

p.set_val('traj.phase0.t_initial', 0.0)
p.set_val('traj.phase0.t_duration', 150.0)
p.set_val('traj.phase0.states:x', phase0.interp('x', [0, 80.0E3]))
p.set_val('traj.phase0.states:y', phase0.interp('y', [0, 50.0E3]))
p.set_val('traj.phase0.states:vx', phase0.interp('vx', [0.0, 2000.0]))
p.set_val('traj.phase0.states:vy', phase0.interp('vy', [10.0, 0.0]))
p.set_val('traj.phase0.states:m', phase0.interp('m', [1500.0, 500.0]))
p.set_val('traj.phase0.controls:theta', phase0.interp('theta', [89.0, 0.0]))
p.set_val('traj.phase0.parameters:thrust', 25000, units='N')
#
phase0.set_refine_options(tol=1.0E-2)

dm.run_problem(p, simulate=True, run_driver=True, refine_iteration_limit=1)

#POST
sol = om.CaseReader('dymos_solution.db').get_case('final')
sim = om.CaseReader('dymos_simulation.db').get_case('final')

plot_results([('traj.phase0.timeseries.x', 'traj.phase0.timeseries.y',
               'range (m)', 'altitude (m)'),
              ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.y',
               'time (s)', 'altitude (m)'),
              ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.vx',
               'time (s)', 'xspeed (m/s)'),
              ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.vy',
               'time (s)', 'yspeed (m/s)'),
              ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.theta',
               'time (s)', 'theta (deg)'),
              ('traj.phase0.timeseries.time', 'traj.phase0.timeseries.m',
               'time (s)', 'm (kg)'),
               ],
             title='Supersonic Minimum Time-to-Climb Solution',
             p_sol=sol, p_sim=sim)

plt.show()

两个阶段的代码只是将其分成两部分并链接它们：

import matplotlib.pyplot as plt
import openmdao.api as om
import dymos as dm
from dymos.examples.plotting import plot_results

p = om.Problem(model=om.Group())
p.driver = om.pyOptSparseDriver()

p.driver.options['optimizer'] = 'IPOPT'

p.driver.opt_settings['derivative_test'] = 'first-order'
p.driver.opt_settings['mu_strategy'] = 'adaptive'
p.driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'
p.driver.opt_settings['bound_mult_init_method'] = 'mu-based'
p.driver.opt_settings['mu_init'] = 0.01
p.driver.opt_settings['nlp_scaling_method'] = 'gradient-based'
p.driver.opt_settings['max_iter']= 1000
p.driver.opt_settings['tol']= 1e-2
p.driver.declare_coloring()

from launch_vehicle_ode import LaunchVehicleODE

traj = dm.Trajectory()

phase0 = dm.Phase(ode_class=LaunchVehicleODE, ode_init_kwargs={'CD': 1.0, 'S': 0.2},
                 transcription=dm.Radau(num_segments=10, order=3, compressed=False))
phase1 = dm.Phase(ode_class=LaunchVehicleODE, ode_init_kwargs={'CD': 1.0, 'S': 0.2},
                 transcription=dm.Radau(num_segments=10, order=3, compressed=False))

traj.add_phase('phase0', phase0)
traj.add_phase('phase1', phase1)

p.model.add_subsystem('traj', traj)

phase0.set_time_options(fix_initial=True, duration_bounds=(1, 500))
phase0.add_state('x',  fix_initial=True, rate_source='xdot', lower = 0.0)
phase0.add_state('y',  fix_initial=True,rate_source='ydot', lower = 0.0)
phase0.add_state('vx', fix_initial=True,rate_source='vxdot', lower = 0.0)
phase0.add_state('vy', fix_initial=True, rate_source='vydot', lower = 0.0)
phase0.add_state('m',  rate_source='mdot', lower=1100.0)
phase0.add_control('theta', fix_initial=True, units='deg', lower=-90.0, upper=90.0, targets=['theta'])
phase0.add_parameter('thrust', val = 25000.0, opt=False, targets=['thrust'])

phase1.set_time_options(duration_bounds=(1, 500))
phase1.add_state('x', rate_source='xdot',  fix_final=True, lower = 0.0)
phase1.add_state('y', rate_source='ydot',  fix_final=True, lower = 0.0)
phase1.add_state('vx',rate_source='vxdot', fix_final=True, lower = 0.0)
phase1.add_state('vy',rate_source='vydot', fix_final=True, lower = 0.0)
phase1.add_state('m', rate_source='mdot',   fix_final=True, lower = 500.0)
phase1.add_control('theta', units='deg', lower=-90.0, upper=90.0, targets=['theta'])
phase1.add_parameter('thrust', val = 25000.0, opt=False, targets=['thrust'])
#
# Set the options for our constraints and objective
#

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

p.model.linear_solver = om.DirectSolver()

p.setup(check=True)

p.set_val('traj.phase0.t_initial', 0.0)
p.set_val('traj.phase0.t_duration', 150.0)
p.set_val('traj.phase0.states:x', phase0.interp('x',   [0, 50.0E3]))
p.set_val('traj.phase0.states:y', phase0.interp('y',   [0, 45.0E3]))
p.set_val('traj.phase0.states:vx', phase0.interp('vx', [0.0, 1000.0]))
p.set_val('traj.phase0.states:vy', phase0.interp('vy', [10.0, 0.0]))
p.set_val('traj.phase0.states:m', phase0.interp('m',   [1500.0, 1100.0]))
p.set_val('traj.phase0.controls:theta', phase0.interp('theta', [89.0, 0.0]))
p.set_val('traj.phase0.parameters:thrust', 25000, units='N')

p.set_val('traj.phase1.t_initial', 100.0)
p.set_val('traj.phase1.t_duration', 150.0)
p.set_val('traj.phase1.states:x', phase1.interp('x',   [50.0E3, 80.0E3]))
p.set_val('traj.phase1.states:y', phase1.interp('y',   [45.0E3, 50.0E3]))
p.set_val('traj.phase1.states:vx', phase1.interp('vx', [1000.0, 2000.0]))
p.set_val('traj.phase1.states:vy', phase1.interp('vy', [0.0, 0.0]))
p.set_val('traj.phase1.states:m', phase1.interp('m',   [1100.0, 500.0]))
p.set_val('traj.phase1.controls:theta', phase1.interp('theta', [0.0, 0.0]))
p.set_val('traj.phase1.parameters:thrust', 25000, units='N')

traj.link_phases(phases=['phase0', 'phase1'], vars=['*'], connected=True)

phase0.set_refine_options(tol=1.0E-2)
phase1.set_refine_options(tol=1.0E-2)
dm.run_problem(p, simulate=False, run_driver=True, refine_iteration_limit=1, 
               simulate_kwargs={'times_per_seg' : 100, 'rtol': 1.0E-2})

p.driver.scaling_report()

sol = om.CaseReader('dymos_solution.db').get_case('final')

fig, axes = plt.subplots(nrows=1, ncols=1, figsize=(10, 6))

time_imp = {'phase0': p.get_val('traj.phase0.timeseries.time'),
            'phase1': p.get_val('traj.phase1.timeseries.time')}
r_imp = {'phase0': p.get_val('traj.phase0.timeseries.x'),
         'phase1': p.get_val('traj.phase1.timeseries.x')}
h_imp = {'phase0': p.get_val('traj.phase0.timeseries.y'),
         'phase1': p.get_val('traj.phase1.timeseries.y')}

axes.plot(r_imp['phase0'], h_imp['phase0'], 'bo')
axes.plot(r_imp['phase1'], h_imp['phase1'], 'ro')

axes.set_xlabel('range (m)')
axes.set_ylabel('altitude (m)')
axes.grid(True)

fig, axes = plt.subplots(nrows=6, ncols=1, figsize=(10, 6))
states = ['x', 'y', 'vx', 'vy', 'm','theta']
for i, state in enumerate(states):
    x_imp = {'phase0': sol.get_val(f'traj.phase0.timeseries.{state}'),
             'phase1': sol.get_val(f'traj.phase1.timeseries.{state}')}
    axes[i].set_ylabel(state)
    axes[i].grid(True)
    axes[i].plot(time_imp['phase0'], x_imp['phase0'], 'bo')
    axes[i].plot(time_imp['phase1'], x_imp['phase1'], 'ro')

plt.show()

用connected=True链接阶段就可以了，我得到了收敛。 尽管如此，我的目的是质量“跳跃”，以模拟两级火箭中燃烧的第一级的下落。所以我正在努力：

traj.link_phases(phases=['phase0', 'phase1'], vars=['time','x','y','vx','vy'], connected=True)
traj.add_linkage_constraint('phase0', 'phase1', 'm', 'm', loc_a='final', loc_b='initial',
                            mult_a=1.0, mult_b=-1.0, equals = 200.0, units='kg', connected=True)

虽然它似乎解释了正确的质量“跳跃”，但我无法让它收敛。我想删除质量中的“连接”设置（它现在按预期完全连续），但要么我在第二阶段没有收敛和/或负时间，要么出现绘图错误。

我很困惑为什么这个问题看起来如此困难。有线索吗？

谢谢你