处理 Dymos 运动方程中的奇点的最佳方法？
Best way to handle singularaties in Dymos equations of motion?
大家好，我正在使用 Dymos 进行火箭弹道优化。我模型的运动方程适用于球形非旋转地球
其中 v 是火箭的速度，gamma 是与当地水平面的飞行路径角度，h 是表面上方的高度，theta 是沿轨迹大圆的角度，m 是质量，epsilon 是攻角，T是信任，D 是阻力，L 是升力。这些与 Dymos 水火箭示例中的运动方程非常相似。我使用了 check_partials，我的解析导数似乎是正确的。
如您所见，伽马的变化率除以 v，因此在轨迹开始处有一个奇点，其中 v=0。您可以通过以某个较小的数字开始 v 来解决这个问题，该数字是预期最终速度的一小部分 - 例如轨道火箭的 1 到 m/s。
但是，我设置的初始速度越低，我就越难看到优化器和网格细化需要工作才能获得收敛的解决方案。例如，
设置 h=200 km，v= 7668 m/s 和 gamma=0（圆形轨道）的最终条件，v initial = 10 m/s，我必须从 20 段开始6阶才能达到收敛，最终的网格是这样的
段数 = 25
段结束 = [-1.0, -0.9833, -0.9667, -0.95, -0.9333, -0.9, -0.85, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0 , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
段顺序 = [9, 9, 9, 9, 9, 9, 9, 9, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6]
此聚合网格的优化时间为 45 秒
如果我使用 v_initial = 1000 m/s，相同的最终条件，以及 20 个 3 阶段的起始网格，最终网格是
段结束 = [-1.0, -0.95, -0.9, -0.8, -0.7, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0]
段顺序 = [6, 6, 6, 6, 6, 3, 3, 3, 3, 3, 3, 3, 3]
收敛网格的优化时间为 3.18 秒。
考虑到网格点集中到轨迹起点的方式，我认为由 v 引起的 gamma_dot 中的奇异性是问题所在。
我可以做些什么来提高性能？
import matplotlib.pyplot as plt
import openmdao.api as om
from openmdao.utils.assert_utils import assert_near_equal
import dymos as dm
import numpy as np
from dymos.examples.plotting import plot_results
import sys
from rocket_ode import rocketODE
#
# Instantiate the problem and configure the optimization driver
#
p = om.Problem(model=om.Group())
p.driver = om.pyOptSparseDriver()
p.driver.options['optimizer'] = 'IPOPT'
p.driver.declare_coloring()
h_target = 200E3
v_circ = (3.986004418E14/(h_target + 6378E3))**0.5
#
# Instantiate the trajectory and phase
#
traj = dm.Trajectory()
seg_ends = [-1.0, -0.9833, -0.9667, -0.95, -0.9333, -0.9, -0.85, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
seg_order = [9, 9, 9, 9, 9, 9, 9, 9, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6]
phase = dm.Phase(ode_class=rocketODE,
transcription=dm.Radau(num_segments=25, order=seg_order,
segment_ends=seg_ends, compressed=False))
traj.add_phase('phase0', phase)
p.model.add_subsystem('traj', traj)
#
# Set the options on the optimization variables
# Note the use of explicit state units here since much of the ODE uses imperial units
# and we prefer to solve this problem using metric units.
# #
phase.set_time_options(fix_initial=True, duration_bounds=(50, 1000))
# phase.set_time_options(fix_initial=True, fix_duration=True, duration_ref=100.0)
phase.add_state('v', fix_initial=True, lower=10.0, units='m/s',
ref=1.0E3, targets=['eom.v'], rate_source='eom.v_dot')
phase.add_state('h', fix_initial=True, lower=0, upper=1.0E6, units='m',
ref=1.0E5, targets=['eom.h'], rate_source='eom.h_dot')
phase.add_state('gamma', fix_initial=True, lower=0, upper=np.pi, units='rad',
ref=1.0, targets=['eom.gamma'], rate_source='eom.gamma_dot')
phase.add_state('lambda', fix_initial=True, lower=0., upper=np.pi, units='rad',
ref=0.01, targets=['eom.lambda'], rate_source='eom.lambda_dot')
phase.add_state('alpha', fix_initial=True, lower=-np.pi / 2, upper=np.pi / 2, units='rad',
ref=1, defect_ref=1, targets=['eom.alpha'],
rate_source='eom.alpha_dot')
phase.add_state('m', fix_initial=True, lower=100.0, upper=1.0E5, units='kg',
ref=1.0E4, targets=['eom.m'], rate_source='eom.m_dot')
phase.add_control('alpha_rate', fix_initial=True, units='rad/s', lower=-0.1, upper=0.1, scaler=1.0, targets=['eom.alpha_rate'])
# phase.add_control('alpha', fix_initial=True, units='deg', lower=-89.0, upper=89.0, scaler=1.0, targets=['eom.alpha'])
# phase.add_control('thrust', units='N', lower=0, upper=1.35E5, ref=1E5, targets=['eom.thrust'])
phase.add_parameter('Isp', val=500.0, units='s', opt=False, targets=['eom.Isp'])
phase.add_parameter('thrust', val=1.3E5, opt=False, targets=['eom.thrust'], units='N')
phase.add_boundary_constraint('h', loc='final', ref=1.0E5, equals=h_target)
phase.add_boundary_constraint('v', loc='final', ref=1.0E3, equals=v_circ)
phase.add_boundary_constraint('gamma', loc='final', equals=0.0, units='rad')
# Minimize time at the end of the phase
phase.add_objective('time', loc='final', scaler=1.0)
p.model.linear_solver = om.DirectSolver()
#
# Setup the problem and set the initial guess
#
p.setup(check=True)
p['traj.phase0.t_initial'] = 0.0
p['traj.phase0.t_duration'] = 300
p['traj.phase0.states:v'] = phase.interpolate(ys=[10.0, v_circ], nodes='state_input')
p['traj.phase0.states:h'] = phase.interpolate(ys=[0.0, h_target], nodes='state_input')
p['traj.phase0.states:gamma'] = phase.interpolate(ys=[np.pi/2, 0.0], nodes='state_input')
p['traj.phase0.states:lambda'] = phase.interpolate(ys=[0.0, 0.02], nodes='state_input')
p['traj.phase0.states:m'] = phase.interpolate(ys=[10000.0, 1000.0], nodes='state_input')
p['traj.phase0.states:alpha'] = phase.interpolate(ys=[0.0, 0.0], nodes='state_input')
p['traj.phase0.controls:alpha_rate'] = phase.interpolate(ys=[0.0, 0.0], nodes='control_input')
# Solve for the optimal trajectory
#
# dm.run_problem(p,refine_method='hp',refine_iteration_limit=20)
dm.run_problem(p)
#
# Get the explicitly simulated solution and plot the results
#
exp_out = traj.simulate(times_per_seg=100,method='RK45')
这是 EOM 代码
import numpy as np
import openmdao.api as om
class rocketEOM(om.ExplicitComponent):
def initialize(self):
self.options.declare('num_nodes', types=int)
self.options.declare('CD', types=float, default=0.5,
desc='coefficient of drag')
self.options.declare('S', types=float, default=7.069,
desc='aerodynamic reference area (m**2)')
self.options.declare('r0', types=float, default=6.3781E6,
desc='Earth radius m')
self.options.declare('g0', types=float, default=9.80665,
desc='acceleration at earth surface')
self.options.declare('omega', types=float, default=7.292115E-5,
desc='rate of Earth rotation, rad/s')
def setup(self):
nn = self.options['num_nodes']
# Inputs
self.add_input('v',
val=np.zeros(nn),
desc='relative velocity',
units='m/s')
self.add_input('h',
val=np.zeros(nn),
desc='altitude',
units='m')
self.add_input('gamma',
val=np.zeros(nn),
desc='flight path angle',
units='rad')
self.add_input('lambda',
val=np.zeros(nn),
desc='longitude',
units='rad')
self.add_input('alpha',
val=np.zeros(nn),
desc='angle of attack',
units='rad')
self.add_input('alpha_rate',
val=np.zeros(nn),
desc='angle of attack rate',
units='rad/s')
self.add_input('m',
val=np.zeros(nn),
desc='mass',
units='kg')
self.add_input('rho',
val=np.zeros(nn),
desc='density',
units='kg/m**3')
self.add_input('thrust',
val=1.4E5 * np.ones(nn),
desc='thrust',
units='N')
self.add_input('Isp',
val=400.0 * np.ones(nn),
desc='specific impulse',
units='s')
# Outputs
self.add_output('v_dot',
val=np.zeros(nn),
desc='acceleraton',
units='m/s**2')
self.add_output('gamma_dot',
val=np.zeros(nn),
desc='flight path angle rate of change',
units='rad/s')
self.add_output('h_dot',
val=np.zeros(nn),
desc='height rate of change',
units='m/s')
self.add_output('lambda_dot',
val=np.zeros(nn),
desc='longitude rate of change',
units='rad/s')
self.add_output('alpha_dot',
val=np.zeros(nn),
desc='angle of attack rate of change',
units='rad/s')
self.add_output('m_dot',
val=np.zeros(nn),
desc='mass rate of change',
units='kg/s')
# Setup partials
ar = np.arange(self.options['num_nodes'])
self.declare_partials(of='v_dot', wrt='h', rows=ar, cols=ar)
self.declare_partials(of='v_dot', wrt='gamma', rows=ar, cols=ar)
self.declare_partials(of='v_dot', wrt='m', rows=ar, cols=ar)
self.declare_partials(of='v_dot', wrt='alpha', rows=ar, cols=ar)
self.declare_partials(of='v_dot', wrt='thrust', rows=ar, cols=ar)
self.declare_partials(of='gamma_dot', wrt='v', rows=ar, cols=ar)
self.declare_partials(of='gamma_dot', wrt='h', rows=ar, cols=ar)
self.declare_partials(of='gamma_dot', wrt='gamma', rows=ar, cols=ar)
self.declare_partials(of='gamma_dot', wrt='m', rows=ar, cols=ar)
self.declare_partials(of='gamma_dot', wrt='alpha', rows=ar, cols=ar)
self.declare_partials(of='gamma_dot', wrt='thrust', rows=ar, cols=ar)
self.declare_partials(of='h_dot', wrt='v', rows=ar, cols=ar)
self.declare_partials(of='h_dot', wrt='gamma', rows=ar, cols=ar)
self.declare_partials(of='lambda_dot', wrt='v', rows=ar, cols=ar)
self.declare_partials(of='lambda_dot', wrt='h', rows=ar, cols=ar)
self.declare_partials(of='lambda_dot', wrt='gamma', rows=ar, cols=ar)
self.declare_partials(of='m_dot', wrt='Isp', rows=ar, cols=ar)
self.declare_partials(of='m_dot', wrt='thrust', rows=ar, cols=ar)
self.declare_partials(of='alpha_dot', wrt='alpha_rate', rows=ar, cols=ar, val=1.0)
def compute(self, inputs, outputs):
alpha = inputs['alpha']
v = inputs['v']
h = inputs['h']
gamma = inputs['gamma']
m = inputs['m']
T = inputs['thrust']
ISP = inputs['Isp']
g = self.g_func(h)
r0 = self.options['r0']
D = 0 # calc drag here
L = 0 # calc lift here
g = self.g_func(h)
r = r0 + h
sin_gamma = np.sin(gamma)
cos_gamma = np.cos(gamma)
v_dot = (T * np.cos(alpha) - D) / m - g * sin_gamma
gamma_dot = ((T * np.sin(alpha) + L) / m - (g - v**2 / r) * cos_gamma) / v
h_dot = v * sin_gamma
lambda_dot = v * cos_gamma / r
m_dot = -T / (9.80665 * ISP)
outputs['v_dot'] = v_dot
outputs['gamma_dot'] = gamma_dot
outputs['h_dot'] = h_dot
outputs['lambda_dot'] = lambda_dot
outputs['m_dot'] = m_dot
outputs['alpha_dot'] = inputs['alpha_rate']
def g_func(self, h):
r0 = self.options['r0']
g0 = self.options['g0']
return g0 * (r0 / (r0 + h))**2
def compute_partials(self, inputs, jacobian):
r0 = self.options['r0']
g0 = self.options['g0']
alpha = inputs['alpha']
v = inputs['v']
h = inputs['h']
gamma = inputs['gamma']
m = inputs['m']
T = inputs['thrust']
ISP = inputs['Isp']
g = self.g_func(h)
r0 = self.options['r0']
D = 0 # calc drag here
L = 0 # calc lift here
# g = self.g_func(h)
r = r0 + h
sin_gamma = np.sin(gamma)
cos_gamma = np.cos(gamma)
sin_alpha = np.sin(alpha)
cos_alpha = np.cos(alpha)
jacobian['v_dot', 'h'] = 2 * g0 * r0**2 * sin_gamma / r**3
jacobian['v_dot', 'gamma'] = -g0 * r0**2 * cos_gamma / r**2
jacobian['v_dot', 'm'] = (D - T * cos_alpha) / m**2
jacobian['v_dot', 'alpha'] = -T * sin_alpha / m
jacobian['v_dot', 'thrust'] = cos_alpha / m
jacobian['gamma_dot', 'v'] = cos_gamma * (g / v**2 + 1 / r) + (L - T * sin_alpha) / (m * v**2)
jacobian['gamma_dot', 'h'] = (2 * g / (r * v) - v / r**2) * cos_gamma
jacobian['gamma_dot', 'gamma'] = ((g0 * r0**2 - v**2 * r) * sin_gamma) / (r**2 * v)
jacobian['gamma_dot', 'm'] = -(-L + T * sin_alpha) / (m**2 * v)
jacobian['gamma_dot', 'alpha'] = (T * cos_alpha) / (m * v)
jacobian['gamma_dot', 'thrust'] = sin_alpha / (m * v)
jacobian['h_dot', 'v'] = sin_gamma
jacobian['h_dot', 'gamma'] = v * cos_gamma
jacobian['lambda_dot', 'v'] = cos_gamma / r
jacobian['lambda_dot', 'h'] = -v * cos_gamma / r**2
jacobian['lambda_dot', 'gamma'] = -v * sin_gamma / r
jacobian['m_dot', 'Isp'] = T / (9.80665 * ISP**2)
jacobian['m_dot', 'thrust'] = - 1 / (9.80665 * ISP)
我同意这里的奇异性是个问题，您将其设置为一些小的初始值的方法通常是一个很好的方法。
虽然它们通常非常合理，但当 EOM 开始变得奇异时，当前在 Dymos 中的网格细化算法可能会出现问题，我想你在这里看到了。
你还能尝试什么？
-您可以将问题分解为多个阶段。第一阶段将涵盖俯仰开始前的垂直上升。通过为轨迹的这个阶段使用一组非奇异的 EOM，您可以让车辆在过渡到第二阶段之前建立一些速度，使用您当前的 EOM。
-您可以重新表述问题以使用在您的操作条件附近没有奇点的 EOM。飞行路径 EOM 通常良好且良性，但它们在 v=0 时存在奇点（在 3D 情况下 gamma=+/-90 度）。
EOM 的笛卡尔公式，其中车辆惯性位置和速度为 [x, y, vx, vy]，将提供仅在行星中心奇异的 EOM。参数化初始状态有点棘手（尤其是在 3D 情况下），因为初始速度是固定在旋转地球上的发射台的初始速度。
您也可以尝试将飞行路径角度作为第一阶段的固定参数，将其固定为 90 度，直到建立足够的速度。在这种情况下，为了获得更好的行为，您将接受动态中的一些小错误。
大家好，我正在使用 Dymos 进行火箭弹道优化。我模型的运动方程适用于球形非旋转地球
其中 v 是火箭的速度，gamma 是与当地水平面的飞行路径角度，h 是表面上方的高度，theta 是沿轨迹大圆的角度，m 是质量，epsilon 是攻角，T是信任，D 是阻力，L 是升力。这些与 Dymos 水火箭示例中的运动方程非常相似。我使用了 check_partials，我的解析导数似乎是正确的。
如您所见，伽马的变化率除以 v，因此在轨迹开始处有一个奇点，其中 v=0。您可以通过以某个较小的数字开始 v 来解决这个问题，该数字是预期最终速度的一小部分 - 例如轨道火箭的 1 到 m/s。
但是，我设置的初始速度越低，我就越难看到优化器和网格细化需要工作才能获得收敛的解决方案。例如，
设置 h=200 km，v= 7668 m/s 和 gamma=0（圆形轨道）的最终条件，v initial = 10 m/s，我必须从 20 段开始6阶才能达到收敛，最终的网格是这样的
段数 = 25 段结束 = [-1.0, -0.9833, -0.9667, -0.95, -0.9333, -0.9, -0.85, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0 , 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0] 段顺序 = [9, 9, 9, 9, 9, 9, 9, 9, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6]
此聚合网格的优化时间为 45 秒
如果我使用 v_initial = 1000 m/s，相同的最终条件，以及 20 个 3 阶段的起始网格，最终网格是
段结束 = [-1.0, -0.95, -0.9, -0.8, -0.7, -0.6, -0.4, -0.2, 0.0, 0.2, 0.4, 0.6, 0.8, 1.0] 段顺序 = [6, 6, 6, 6, 6, 3, 3, 3, 3, 3, 3, 3, 3]
收敛网格的优化时间为 3.18 秒。
考虑到网格点集中到轨迹起点的方式，我认为由 v 引起的 gamma_dot 中的奇异性是问题所在。
我可以做些什么来提高性能？
import matplotlib.pyplot as plt
import openmdao.api as om
from openmdao.utils.assert_utils import assert_near_equal
import dymos as dm
import numpy as np
from dymos.examples.plotting import plot_results
import sys
from rocket_ode import rocketODE
#
# Instantiate the problem and configure the optimization driver
#
p = om.Problem(model=om.Group())
p.driver = om.pyOptSparseDriver()
p.driver.options['optimizer'] = 'IPOPT'
p.driver.declare_coloring()
h_target = 200E3
v_circ = (3.986004418E14/(h_target + 6378E3))**0.5
#
# Instantiate the trajectory and phase
#
traj = dm.Trajectory()
seg_ends = [-1.0, -0.9833, -0.9667, -0.95, -0.9333, -0.9, -0.85, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0]
seg_order = [9, 9, 9, 9, 9, 9, 9, 9, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6, 6]
phase = dm.Phase(ode_class=rocketODE,
transcription=dm.Radau(num_segments=25, order=seg_order,
segment_ends=seg_ends, compressed=False))
traj.add_phase('phase0', phase)
p.model.add_subsystem('traj', traj)
#
# Set the options on the optimization variables
# Note the use of explicit state units here since much of the ODE uses imperial units
# and we prefer to solve this problem using metric units.
# #
phase.set_time_options(fix_initial=True, duration_bounds=(50, 1000))
# phase.set_time_options(fix_initial=True, fix_duration=True, duration_ref=100.0)
phase.add_state('v', fix_initial=True, lower=10.0, units='m/s',
ref=1.0E3, targets=['eom.v'], rate_source='eom.v_dot')
phase.add_state('h', fix_initial=True, lower=0, upper=1.0E6, units='m',
ref=1.0E5, targets=['eom.h'], rate_source='eom.h_dot')
phase.add_state('gamma', fix_initial=True, lower=0, upper=np.pi, units='rad',
ref=1.0, targets=['eom.gamma'], rate_source='eom.gamma_dot')
phase.add_state('lambda', fix_initial=True, lower=0., upper=np.pi, units='rad',
ref=0.01, targets=['eom.lambda'], rate_source='eom.lambda_dot')
phase.add_state('alpha', fix_initial=True, lower=-np.pi / 2, upper=np.pi / 2, units='rad',
ref=1, defect_ref=1, targets=['eom.alpha'],
rate_source='eom.alpha_dot')
phase.add_state('m', fix_initial=True, lower=100.0, upper=1.0E5, units='kg',
ref=1.0E4, targets=['eom.m'], rate_source='eom.m_dot')
phase.add_control('alpha_rate', fix_initial=True, units='rad/s', lower=-0.1, upper=0.1, scaler=1.0, targets=['eom.alpha_rate'])
# phase.add_control('alpha', fix_initial=True, units='deg', lower=-89.0, upper=89.0, scaler=1.0, targets=['eom.alpha'])
# phase.add_control('thrust', units='N', lower=0, upper=1.35E5, ref=1E5, targets=['eom.thrust'])
phase.add_parameter('Isp', val=500.0, units='s', opt=False, targets=['eom.Isp'])
phase.add_parameter('thrust', val=1.3E5, opt=False, targets=['eom.thrust'], units='N')
phase.add_boundary_constraint('h', loc='final', ref=1.0E5, equals=h_target)
phase.add_boundary_constraint('v', loc='final', ref=1.0E3, equals=v_circ)
phase.add_boundary_constraint('gamma', loc='final', equals=0.0, units='rad')
# Minimize time at the end of the phase
phase.add_objective('time', loc='final', scaler=1.0)
p.model.linear_solver = om.DirectSolver()
#
# Setup the problem and set the initial guess
#
p.setup(check=True)
p['traj.phase0.t_initial'] = 0.0
p['traj.phase0.t_duration'] = 300
p['traj.phase0.states:v'] = phase.interpolate(ys=[10.0, v_circ], nodes='state_input')
p['traj.phase0.states:h'] = phase.interpolate(ys=[0.0, h_target], nodes='state_input')
p['traj.phase0.states:gamma'] = phase.interpolate(ys=[np.pi/2, 0.0], nodes='state_input')
p['traj.phase0.states:lambda'] = phase.interpolate(ys=[0.0, 0.02], nodes='state_input')
p['traj.phase0.states:m'] = phase.interpolate(ys=[10000.0, 1000.0], nodes='state_input')
p['traj.phase0.states:alpha'] = phase.interpolate(ys=[0.0, 0.0], nodes='state_input')
p['traj.phase0.controls:alpha_rate'] = phase.interpolate(ys=[0.0, 0.0], nodes='control_input')
# Solve for the optimal trajectory
#
# dm.run_problem(p,refine_method='hp',refine_iteration_limit=20)
dm.run_problem(p)
#
# Get the explicitly simulated solution and plot the results
#
exp_out = traj.simulate(times_per_seg=100,method='RK45')
这是 EOM 代码
import numpy as np
import openmdao.api as om
class rocketEOM(om.ExplicitComponent):
def initialize(self):
self.options.declare('num_nodes', types=int)
self.options.declare('CD', types=float, default=0.5,
desc='coefficient of drag')
self.options.declare('S', types=float, default=7.069,
desc='aerodynamic reference area (m**2)')
self.options.declare('r0', types=float, default=6.3781E6,
desc='Earth radius m')
self.options.declare('g0', types=float, default=9.80665,
desc='acceleration at earth surface')
self.options.declare('omega', types=float, default=7.292115E-5,
desc='rate of Earth rotation, rad/s')
def setup(self):
nn = self.options['num_nodes']
# Inputs
self.add_input('v',
val=np.zeros(nn),
desc='relative velocity',
units='m/s')
self.add_input('h',
val=np.zeros(nn),
desc='altitude',
units='m')
self.add_input('gamma',
val=np.zeros(nn),
desc='flight path angle',
units='rad')
self.add_input('lambda',
val=np.zeros(nn),
desc='longitude',
units='rad')
self.add_input('alpha',
val=np.zeros(nn),
desc='angle of attack',
units='rad')
self.add_input('alpha_rate',
val=np.zeros(nn),
desc='angle of attack rate',
units='rad/s')
self.add_input('m',
val=np.zeros(nn),
desc='mass',
units='kg')
self.add_input('rho',
val=np.zeros(nn),
desc='density',
units='kg/m**3')
self.add_input('thrust',
val=1.4E5 * np.ones(nn),
desc='thrust',
units='N')
self.add_input('Isp',
val=400.0 * np.ones(nn),
desc='specific impulse',
units='s')
# Outputs
self.add_output('v_dot',
val=np.zeros(nn),
desc='acceleraton',
units='m/s**2')
self.add_output('gamma_dot',
val=np.zeros(nn),
desc='flight path angle rate of change',
units='rad/s')
self.add_output('h_dot',
val=np.zeros(nn),
desc='height rate of change',
units='m/s')
self.add_output('lambda_dot',
val=np.zeros(nn),
desc='longitude rate of change',
units='rad/s')
self.add_output('alpha_dot',
val=np.zeros(nn),
desc='angle of attack rate of change',
units='rad/s')
self.add_output('m_dot',
val=np.zeros(nn),
desc='mass rate of change',
units='kg/s')
# Setup partials
ar = np.arange(self.options['num_nodes'])
self.declare_partials(of='v_dot', wrt='h', rows=ar, cols=ar)
self.declare_partials(of='v_dot', wrt='gamma', rows=ar, cols=ar)
self.declare_partials(of='v_dot', wrt='m', rows=ar, cols=ar)
self.declare_partials(of='v_dot', wrt='alpha', rows=ar, cols=ar)
self.declare_partials(of='v_dot', wrt='thrust', rows=ar, cols=ar)
self.declare_partials(of='gamma_dot', wrt='v', rows=ar, cols=ar)
self.declare_partials(of='gamma_dot', wrt='h', rows=ar, cols=ar)
self.declare_partials(of='gamma_dot', wrt='gamma', rows=ar, cols=ar)
self.declare_partials(of='gamma_dot', wrt='m', rows=ar, cols=ar)
self.declare_partials(of='gamma_dot', wrt='alpha', rows=ar, cols=ar)
self.declare_partials(of='gamma_dot', wrt='thrust', rows=ar, cols=ar)
self.declare_partials(of='h_dot', wrt='v', rows=ar, cols=ar)
self.declare_partials(of='h_dot', wrt='gamma', rows=ar, cols=ar)
self.declare_partials(of='lambda_dot', wrt='v', rows=ar, cols=ar)
self.declare_partials(of='lambda_dot', wrt='h', rows=ar, cols=ar)
self.declare_partials(of='lambda_dot', wrt='gamma', rows=ar, cols=ar)
self.declare_partials(of='m_dot', wrt='Isp', rows=ar, cols=ar)
self.declare_partials(of='m_dot', wrt='thrust', rows=ar, cols=ar)
self.declare_partials(of='alpha_dot', wrt='alpha_rate', rows=ar, cols=ar, val=1.0)
def compute(self, inputs, outputs):
alpha = inputs['alpha']
v = inputs['v']
h = inputs['h']
gamma = inputs['gamma']
m = inputs['m']
T = inputs['thrust']
ISP = inputs['Isp']
g = self.g_func(h)
r0 = self.options['r0']
D = 0 # calc drag here
L = 0 # calc lift here
g = self.g_func(h)
r = r0 + h
sin_gamma = np.sin(gamma)
cos_gamma = np.cos(gamma)
v_dot = (T * np.cos(alpha) - D) / m - g * sin_gamma
gamma_dot = ((T * np.sin(alpha) + L) / m - (g - v**2 / r) * cos_gamma) / v
h_dot = v * sin_gamma
lambda_dot = v * cos_gamma / r
m_dot = -T / (9.80665 * ISP)
outputs['v_dot'] = v_dot
outputs['gamma_dot'] = gamma_dot
outputs['h_dot'] = h_dot
outputs['lambda_dot'] = lambda_dot
outputs['m_dot'] = m_dot
outputs['alpha_dot'] = inputs['alpha_rate']
def g_func(self, h):
r0 = self.options['r0']
g0 = self.options['g0']
return g0 * (r0 / (r0 + h))**2
def compute_partials(self, inputs, jacobian):
r0 = self.options['r0']
g0 = self.options['g0']
alpha = inputs['alpha']
v = inputs['v']
h = inputs['h']
gamma = inputs['gamma']
m = inputs['m']
T = inputs['thrust']
ISP = inputs['Isp']
g = self.g_func(h)
r0 = self.options['r0']
D = 0 # calc drag here
L = 0 # calc lift here
# g = self.g_func(h)
r = r0 + h
sin_gamma = np.sin(gamma)
cos_gamma = np.cos(gamma)
sin_alpha = np.sin(alpha)
cos_alpha = np.cos(alpha)
jacobian['v_dot', 'h'] = 2 * g0 * r0**2 * sin_gamma / r**3
jacobian['v_dot', 'gamma'] = -g0 * r0**2 * cos_gamma / r**2
jacobian['v_dot', 'm'] = (D - T * cos_alpha) / m**2
jacobian['v_dot', 'alpha'] = -T * sin_alpha / m
jacobian['v_dot', 'thrust'] = cos_alpha / m
jacobian['gamma_dot', 'v'] = cos_gamma * (g / v**2 + 1 / r) + (L - T * sin_alpha) / (m * v**2)
jacobian['gamma_dot', 'h'] = (2 * g / (r * v) - v / r**2) * cos_gamma
jacobian['gamma_dot', 'gamma'] = ((g0 * r0**2 - v**2 * r) * sin_gamma) / (r**2 * v)
jacobian['gamma_dot', 'm'] = -(-L + T * sin_alpha) / (m**2 * v)
jacobian['gamma_dot', 'alpha'] = (T * cos_alpha) / (m * v)
jacobian['gamma_dot', 'thrust'] = sin_alpha / (m * v)
jacobian['h_dot', 'v'] = sin_gamma
jacobian['h_dot', 'gamma'] = v * cos_gamma
jacobian['lambda_dot', 'v'] = cos_gamma / r
jacobian['lambda_dot', 'h'] = -v * cos_gamma / r**2
jacobian['lambda_dot', 'gamma'] = -v * sin_gamma / r
jacobian['m_dot', 'Isp'] = T / (9.80665 * ISP**2)
jacobian['m_dot', 'thrust'] = - 1 / (9.80665 * ISP)
我同意这里的奇异性是个问题，您将其设置为一些小的初始值的方法通常是一个很好的方法。
虽然它们通常非常合理，但当 EOM 开始变得奇异时，当前在 Dymos 中的网格细化算法可能会出现问题，我想你在这里看到了。
你还能尝试什么？
-您可以将问题分解为多个阶段。第一阶段将涵盖俯仰开始前的垂直上升。通过为轨迹的这个阶段使用一组非奇异的 EOM，您可以让车辆在过渡到第二阶段之前建立一些速度，使用您当前的 EOM。
-您可以重新表述问题以使用在您的操作条件附近没有奇点的 EOM。飞行路径 EOM 通常良好且良性，但它们在 v=0 时存在奇点（在 3D 情况下 gamma=+/-90 度）。
EOM 的笛卡尔公式，其中车辆惯性位置和速度为 [x, y, vx, vy]，将提供仅在行星中心奇异的 EOM。参数化初始状态有点棘手（尤其是在 3D 情况下），因为初始速度是固定在旋转地球上的发射台的初始速度。
您也可以尝试将飞行路径角度作为第一阶段的固定参数，将其固定为 90 度，直到建立足够的速度。在这种情况下，为了获得更好的行为，您将接受动态中的一些小错误。