Лучший способ справиться с сингулярностями в уравнениях движения Даймоса?

Привет всем, я работаю над оптимизацией траектории ракеты с Dymos. Уравнения движения моей модели относятся к сферической невращающейся Земле.

где v — скорость ракеты, gamma — угол траектории полета к местному горизонту, h — высота над поверхностью, тета — угол по большому кругу траектории, m — масса, эпсилон — угол атаки, T — доверие, D — сопротивление, L подъемная сила. Они очень похожи на уравнения движения в примере с водяной ракетой Dymos. Я использовал check_partials, и мои аналитические производные кажутся правильными.

Как видите, скорость изменения гаммы делится на v, поэтому у вас есть сингулярность, где v=0 в начале траектории. Вы можете обойти это, начав v с некоторого небольшого числа, которое составляет небольшую часть предполагаемой конечной скорости, например, от 1 до м/с для орбитальной ракеты.

Однако чем ниже я устанавливаю начальную скорость, тем труднее я вижу, что оптимизатор и уточнение сетки должны работать, чтобы получить конвергентное решение. Например,

установив конечные условия h = 200 км, v = 7668 м / с и gamma = 0 (круговая орбита), с начальной v = 10 м / с, я должен начать с 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 м/с, те же конечные условия и начальную сетку из 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 секунды.

Учитывая то, как точки сетки концентрируются в начале траектории, я думаю, что сингулярность в gamma_dot, вызванная v, является проблемой.

Что я могу сделать, чтобы улучшить производительность?

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')

Это код ЕОМ

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)