Dymos导出导数/雅可比矩阵时KeyError问题求助

我尝试从evtol_dymos_vectorized.py文件中导出梯度信息，参考相关方法并按OpenMDAO的指引记录导数，却遇到错误：KeyError: 'traj.phases.phase0->final_boundary_constraint->x'。

我的实现代码 #

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

import sys

sys.path.insert(0, "../ode")

from evtol_dynamics_comp_vectorized import Dynamics as DynamicsVectorized

import verify_data

if __name__ == '__main__':
    input_arg_1 = 0.0

    input_arg_2 = 'ns'

    # Some specifications
    prop_rad = 0.75
    wing_S = 9.
    wing_span = 6.
    num_blades = 3.
    blade_chord = 0.1
    num_props = 8

    # User-specified input dictionary
    input_dict = {'T_guess': 9.8 * 725 * 1.2,  # initial thrust guess
                  'x_dot_initial': 0.,  # initial horizontal speed
                  'y_dot_initial': 0.01,  # initial vertical speed
                  'y_initial': 0.01,  # initial vertical displacement
                  'A_disk': np.pi * prop_rad ** 2 * num_props,  # total propeller disk area
                  'AR': wing_span ** 2 / (0.5 * wing_S),  # aspect ratio of each wing
                  'e': 0.68,  # span efficiency factor of each wing
                  't_over_c': 0.12,  # airfoil thickness-to-chord ratio
                  'S': wing_S,  # total wing reference area
                  'CD0': 0.35 / wing_S,  # coefficient of drag of the fuselage, gear, etc.
                  'm': 725.,  # mass of aircraft
                  'a0': 5.9,  # airfoil lift-curve slope
                  'alpha_stall': 15. / 180. * np.pi,  # wing stall angle
                  'rho': 1.225,  # air density
                  'induced_velocity_factor': int(input_arg_1) / 100.,  # induced-velocity factor
                  'stall_option': input_arg_2,  # stall option: 's' allows stall, 'ns' does not
                  'R': prop_rad,  # propeller radius
                  'solidity': num_blades * blade_chord / np.pi / prop_rad,  # solidity
                  'omega': 136. / prop_rad,  # angular rotation rate
                  'prop_CD0': 0.012,  # CD0 for prop profile power
                  'k_elec': 0.9,  # electrical and mechanical losses factor
                  'k_ind': 1.2,  # induced-losses factor
                  'nB': num_blades,  # number of blades per propeller
                  'bc': blade_chord,  # representative blade chord
                  'n_props': num_props  # number of propellers
                  }

    p = om.Problem()

    traj = dm.Trajectory()
    p.model.add_subsystem('traj', traj)

    phase = dm.Phase(transcription=dm.GaussLobatto(num_segments=10, order=3, solve_segments=False,
                                                   compressed=False),
                     ode_class=DynamicsVectorized,
                     ode_init_kwargs={'input_dict': input_dict})

    traj.add_phase('phase0', phase)

    phase.set_time_options(fix_initial=True, duration_bounds=(5, 60), duration_ref=30)
    phase.add_state('x', fix_initial=True, rate_source='x_dot', ref0=0, ref=900, defect_ref=100)
    phase.add_state('y', fix_initial=True, rate_source='y_dot', ref0=0, ref=300, defect_ref=300)
    phase.add_state('vx', fix_initial=True, rate_source='a_x', ref0=0, ref=10)
    phase.add_state('vy', fix_initial=True, rate_source='a_y', ref0=0, ref=10)
    phase.add_state('energy', fix_initial=True, rate_source='energy_dot', ref0=0, ref=1E7, defect_ref=1E5)

    phase.add_control('power', lower=1e3, upper=311000, ref0=1e3, ref=311000, rate_continuity=False)
    phase.add_control('theta', lower=0., upper=3 * np.pi / 4, ref0=0, ref=3 * np.pi / 4,
                      rate_continuity=False)

    phase.add_timeseries_output(['CL', 'CD'])

    # Objective
    phase.add_objective('energy', loc='final', ref0=0, ref=1E7)

    # Boundary Constraints
    phase.add_boundary_constraint('y', loc='final', lower=305,
                                  ref=100)  # Constraint for the final vertical displacement
    phase.add_boundary_constraint('x', loc='final', equals=900,
                                  ref=100)  # Constraint for the final horizontal displacement
    phase.add_boundary_constraint('x_dot', loc='final', equals=67.,
                                  ref=100)  # Constraint for the final horizontal speed

    # Path Constraints
    phase.add_path_constraint('y', lower=0., upper=305,# was 305
                              ref=300)  # Constraint for the minimum vertical displacement
    phase.add_path_constraint('acc', upper=0.3,
                              ref=1.0)  # Constraint for the acceleration magnitude
    phase.add_path_constraint('aoa', lower=-np.radians(15), upper=np.radians(15), ref0=-np.radians(15),
                              ref=np.radians(15))  # Constraint for the angle of attack
    phase.add_path_constraint('thrust', lower=10, ref0=10,
                              ref=100)  # Constraint for the thrust magnitude

    # # Setup the driver
    driver = p.driver = om.pyOptSparseDriver()

    # p.driver.options['optimizer'] = 'SNOPT'
    # p.driver.opt_settings['Major optimality tolerance'] = 1e-4
    # p.driver.opt_settings['Major feasibility tolerance'] = 1e-6
    # p.driver.opt_settings['Major iterations limit'] = 1000
    # p.driver.opt_settings['Minor iterations limit'] = 100_000_000
    # p.driver.opt_settings['iSumm'] = 6

    driver.options['optimizer'] = 'IPOPT'
    driver.opt_settings['max_iter'] = 1000
    driver.opt_settings['alpha_for_y'] = 'safer-min-dual-infeas'
    driver.opt_settings['print_level'] = 5
    driver.opt_settings['nlp_scaling_method'] = 'gradient-based'
    driver.opt_settings['tol'] = 5.0E-5

    driver.declare_coloring(tol=1.0E-8)
    
    ## Added this to record the derivatives (gradient info)
    driver.recording_options['record_derivatives'] = True
    driver.add_recorder(om.SqliteRecorder('cases.sql'))

    p.setup()
    p.set_solver_print(0)

    p.set_val('traj.phase0.t_initial', 0.0)
    p.set_val('traj.phase0.t_duration', 30)
    p.set_val('traj.phase0.states:x', phase.interpolate(ys=[0, 900], nodes='state_input'))
    p.set_val('traj.phase0.states:y', phase.interpolate(ys=[0.01, 300], nodes='state_input'))
    p.set_val('traj.phase0.states:vx', phase.interpolate(ys=[0, 60], nodes='state_input'))
    p.set_val('traj.phase0.states:vy', phase.interpolate(ys=[0.01, 10], nodes='state_input'))
    p.set_val('traj.phase0.states:energy', phase.interpolate(ys=[0, 1E7], nodes='state_input'))

    p.set_val('traj.phase0.controls:power', phase.interpolate(xs=np.linspace(0, 28.368, 500),
                                                              ys=verify_data.powers.ravel(),
                                                              nodes='control_input'))
    p.set_val('traj.phase0.controls:theta', phase.interpolate(xs=np.linspace(0, 28.368, 500),
                                                              ys=verify_data.thetas.ravel(),
                                                              nodes='control_input'))

    p.set_val('traj.phase0.controls:power', 200000.0)
    p.set_val('traj.phase0.controls:theta', phase.interpolate(ys=[0.001, np.radians(85)], nodes='control_input'))


    dm.run_problem(p, run_driver=True, simulate=True)
    
    p.cleanup()

    
## Added this section to access the recorded derivaives
cr = om.CaseReader('cases.sql')

# Get derivatives associated with the last iteration.
deriv_1 = cr.get_case(-1).derivatives

print("deriv_1 = ",deriv_1)

报错输出 #

File "/home/sps79w/anaconda3/envs/Aero_Optimize/lib/python3.9/site-packages/spyder_kernels/py3compat.py", line 356, in compat_exec
    exec(code, globals, locals)

  File "/home/sps79w/Documents/PhD_Code/RevHack2020-master/problems/evtol_trajectory/evtol_dymos/evtol_dymos_vectorized.py", line 153, in <module>
    cr = om.CaseReader('cases.sql')

  File "/home/sps79w/anaconda3/envs/Aero_Optimize/lib/python3.9/site-packages/openmdao/recorders/case_reader.py", line 26, in CaseReader
    return SqliteCaseReader(filename, pre_load, metadata_filename)

  File "/home/sps79w/anaconda3/envs/Aero_Optimize/lib/python3.9/site-packages/openmdao/recorders/sqlite_reader.py", line 168, in __init__
    self._load_cases()

  File "/home/sps79w/anaconda3/envs/Aero_Optimize/lib/python3.9/site-packages/openmdao/recorders/sqlite_reader.py", line 360, in _load_cases
    self._driver_cases._load_cases()

  File "/home/sps79w/anaconda3/envs/Aero_Optimize/lib/python3.9/site-packages/openmdao/recorders/sqlite_reader.py", line 1259, in _load_cases
    for case in self.cases(cache=True):

  File "/home/sps79w/anaconda3/envs/Aero_Optimize/lib/python3.9/site-packages/openmdao/recorders/sqlite_reader.py", line 1420, in cases
    case = Case('driver', row, self._prom2abs, self._abs2prom, self._abs2meta,

  File "/home/sps79w/anaconda3/envs/Aero_Optimize/lib/python3.9/site-packages/openmdao/recorders/case.py", line 191, in __init__
    self.derivatives = PromAbsDict(jacobian, prom2abs['output'], abs2prom['output'],

  File "/home/sps79w/anaconda3/envs/Aero_Optimize/lib/python3.9/site-packages/openmdao/recorders/case.py", line 938, in __init__
    abs_keys, prom_key = self._deriv_keys(key)

  File "/home/sps79w/anaconda3/envs/Aero_Optimize/lib/python3.9/site-packages/openmdao/recorders/case.py", line 984, in _deriv_keys
    abs_of = [of] if of in abs2prom else prom2abs[of]

KeyError: 'traj.phases.phase0->final_boundary_constraint->x'

解决方案 #

这个错误源于Dymos生成的边界约束内部名称与你尝试访问的路径不匹配，或是记录器配置时机不当导致元数据映射错误，可按以下步骤修复：

1. 调整记录器配置时机 #

将记录器的配置代码移到p.setup()之后执行，确保问题所有变量和约束的路径已正确生成：

p.setup()
p.set_solver_print(0)

# 移动记录器配置到setup之后
driver.recording_options['record_derivatives'] = True
driver.add_recorder(om.SqliteRecorder('cases.sql'))

2. 确认约束的实际绝对路径 #

在p.setup()之后打印所有约束的绝对路径，获取正确的访问键名：

p.setup()
# 打印所有约束的绝对路径
for con in p.model.list_constraints(return_format='dict').keys():
    print(f"Constraint absolute path: {con}")

3. 避免提前执行cleanup #

p.cleanup()会清除问题的元数据，导致CaseReader无法正确解析记录的导数，直接移除该行代码即可。

4. 版本兼容性检查 #

确保OpenMDAO和Dymos版本兼容（例如OpenMDAO 3.x搭配Dymos 1.x），版本不匹配也可能导致此类路径映射问题。

修改后的核心代码片段：

# ... 其他代码保持不变 ...

p.setup()
p.set_solver_print(0)

# 配置记录器
driver.recording_options['record_derivatives'] = True
driver.add_recorder(om.SqliteRecorder('cases.sql'))

# 打印约束路径（可选）
for con in p.model.list_constraints(return_format='dict').keys():
    print(f"Constraint path: {con}")

# ... 变量初始化代码保持不变 ...

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

# 移除p.cleanup()

# 访问导数
cr = om.CaseReader('cases.sql')
deriv_1 = cr.get_case(-1).derivatives
print("deriv_1 = ", deriv_1)

内容的提问来源于stack exchange，提问作者MSTGrad