Как лучше всего представить Матрицу переходов состояний и как использовать ее для поиска периодических гало-орбит?

Question

Как лучше всего представить Матрицу переходов состояний и как использовать ее для поиска периодических гало-орбит?

ооо

Сначала я изложу свой математический вопрос о матрицах распространения состояния и перехода состояния , а затем покажу вам простую задачу, для которой я хотел бы использовать эти концепции для создания семейства плотно расположенных гало-орбит.

Я также начну с утверждения, что я ищу Aha! напечатайте ответ. Я не надеюсь на объяснение, пока это превосходное интуитивное объяснение кватернионов . Мне не нужно, чтобы все было проработано, просто какое-то объяснение того, как можно понять, получить и использовать Матрицу переходов состояний в этом контексте.

Следующее довольно стандартно, я цитирую из статьи, которая оказалась у меня под рукой, Хуан Сенент, Сезар Окампо и Антонио Капелла; Переходы с малой тягой и переменным импульсом и наведение на нестабильные периодические орбиты. Журнал руководства, контроля и динамики, 28 (2) марта – апреля 2005 г .:

Для динамической системы

$\dot{Икс} знак равно ф (Икс)$ $\mathbf{\dot x} = \mathbf{f}(\mathbf{x})$

оценивается из $t_0=0$ некоторым $t=t_f$ , дифференциал конечного состояния при $t_f$ дан кем-то

$д {Икс}_{ф} знак равно Φ (т_{ф}, т_{0}) дельта {Икс}_{0} + {\dot{Икс}}_{ф} д т_{ф}$ $\text{d} \mathbf{x}_f = \mathbf{\Phi}(t_f, t_0) \delta \mathbf{x}_0 + \mathbf{\dot x}_f \text{d} t_f$

где матрица перехода состояний удовлетворяет

$\dot{Φ} (т, т_{0}) знак равно Ф (Икс (т)) Φ (т, т_{0})$ $\mathbf{\dot \Phi} (t,t_0) = \mathbf{F}(\mathbf{x}(t)) \mathbf{\Phi}(t, t_0)$

а также

$Φ (т_{0}, т_{0}) знак равно я_{6 \times 6}$ $\mathbf{\Phi} (t_0, t_0) = \mathbf{I}_{6 \times 6}$

а также $\mathbf{F}$ - якобиан векторного поля, используемого в качестве матрицы распространения состояния,

$Ф (Икс (т)) знак равно \frac{\partial ф (Икс)}{\partial Икс}$ $\mathbf{F}(\mathbf{x}(t)) = \frac{\partial\mathbf{f}(\mathbf{x})}{\partial \mathbf{x}}$

Я начал с классической статьи , написанной Кэтлин Коннор Хауэлл «Трехмерные , периодические «гало» орбиты небесной механики 32 (1984) 53-71». В нем описывается метод поиска решений для гало-орбит в круговой ограниченной задаче трех тел (CR3BP), очень близкий к методу, описанному Бреквеллом, Дж. В. и Брауном, Дж. В.: 1979, Семейство трехмерных периодических орбит «Гало». в ограниченной задаче трех тел Земля-Луна Целест. мех. 20 , 389.

Howell, 1984, подробно описывает пошаговую процедуру нахождения членов семейства гало-орбит относительно коллинеарных точек либрации Лагранжа, обладающих симметрией относительно плоскости xz, используя тот факт, что для этой группы орбит три из шести компонент вектора состояния должны сходиться к нулю в точке пересечения орбиты с плоскостью.

В статье приведены шесть примеров гало-орбит, и с приведенными там числами я могу проинтегрировать векторы состояния, убедиться, что три компонента вектора состояния $y, v_x, v_z$ действительно проходят через ноль в середине орбиты и составляют хороший график.

Что я хотел бы сделать, так это интуитивно понять, что такое вектор распространения состояния и вектор перехода состояния, и как использовать их для более быстрой сходимости на новом члене семейства гало-орбит, чем если бы я только начал снимать орбиты в кластере. вокруг начальной точки и использовал что-то простое, например, наискорейший спуск, чтобы найти следующую орбиту с $y, v_x, v_z$ все равно нулю.

\ddot{Икс} знак равно Икс + 2 \dot{у} - \frac{(1 - мю) (Икс + мю)}{р_{1}^{3}} - \frac{мю (Икс - 1 + мю)}{р_{2}^{3}}

$\ddot{x}=x+2\dot{y}-\frac{(1-\mu)(x+\mu)}{r_1^3}-\frac{\mu(x-1+\mu)}{r_2^3}$

\ddot{у} знак равно - 2 \dot{Икс} + у (1 - \frac{1 - мю}{р_{1}^{3}} - \frac{мю}{р_{2}^{3}})

$\ddot{y}=-2\dot{x}+y\left( 1-\frac{1-\mu}{r_1^3} -\frac{\mu}{r_2^3}\right)$

\ddot{г} знак равно - г (\frac{1 - мю}{р_{1}^{3}} + \frac{мю}{р_{2}^{3}})

$\ddot{z}=-z\left( \frac{1-\mu}{r_1^3} + \frac{\mu}{r_2^3} \right)$

куда

р_{1} знак равно \sqrt{(Икс + мю)^{2} + у^{2} + г^{2}}

$r_1=\sqrt{(x+\mu)^2 + y^2 + z^2}$

р_{2} знак равно \sqrt{(Икс - 1 + мю)^{2} + у^{2} + г^{2}}

$r_2=\sqrt{(x-1+\mu)^2 + y^2 + z^2}$

ПРИМЕЧАНИЕ! Я считаю, что метки для позиций L ${}_1$ и я ${}_2$ в GIF и скрипт переставлены (неправильные метки/названия). Я скоро обновлю изображение.

def deriv(X, t):
    x, y, z, xdot, ydot, zdot = X
    r1 = np.sqrt((x      + mu)**2 + y**2 + z**2)
    r2 = np.sqrt((x - 1. + mu)**2 + y**2 + z**2)

    term_1 = x + 2. * ydot
    term_2 = -(1.-mu) * (x + mu) / r1**3
    term_3 =     -mu  * (x - 1. + mu) / r2**3
    xddot  = term_1 + term_2 + term_3

    term_1 = -2. * xdot
    term_2 = 1. - (1.-mu)/r1**3 - mu/r2**3 
    yddot  = term_1 + y * term_2

    term_1 = (1. - mu)/r1**3 + mu/r2**3  # should be plus???
    zddot  = -z * term_1

    return np.array([xdot, ydot, zdot, xddot, yddot, zddot])


class Sat(object):
    def __init__(self, X0, T0, nu12):
        self.X0 = X0
        self.pos0 = X0[:3]
        self.v0   = X0[3:]
        self.T0 = T0
        self.nu1, self.nu2 = nu12       


import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint
from mpl_toolkits.mplot3d import Axes3D

# From "Three-Dimensional, Periodic 'Halo' Orbits,
# Kathleen Connor Howell, Celestial Mechanics 32 (1984) 53-71 

pi, twopi = np.pi, 2*np.pi
mu = 0.04

# starting points:
x0     =   [0.723268, 0.729988, 0.753700, 0.777413, 0.801125, 0.817724]
y0     = 6*[0.0]
z0     =   [0.040000, 0.215589, 0.267595, 0.284268, 0.299382, 0.313788]
xdot0  = 6*[0.0]
ydot0  =   [0.198019, 0.397259, 0.399909, 0.361870, 0.312474, 0.271306]
zdot0  = 6*[0.0]

# X0s    = np.array(zip(x0, y0, z0, xdot0, ydot0, zdot0)) Nope! 
X0s    = np.array(list(zip(x0, y0, z0, xdot0, ydot0, zdot0)))

Thalf0s = [1.300177, 1.348532, 1.211253, 1.101099, 1.017241, 0.978653]
T0s     = [2.0*x for x in Thalf0s]

nu1s    = [1181.69,    51.07839,  4.95816,  1.101843,  0.94834,  1.10361]
nu2s    = [   0.98095, -0.90203, -0.40587, -0.420200, -1.58429, -2.09182]
nu12s   = zip(nu1s, nu2s)

n_half  = 200
fractional_times  = np.linspace(0.0, 1.0, 2*n_half+1)

rtol, atol = 1E-12, 1E-12

sats   = []
for X0, T0, nu12 in zip(X0s, T0s, nu12s):
    sat = Sat(X0, T0, nu12)
    sat.n_half  = n_half
    sat.t = sat.T0 * fractional_times
    sat.rtol, sat.atol = rtol, atol    
    sats.append(sat)

for sat in sats:
    answer, info = ODEint(deriv, sat.X0, sat.t,
                          rtol=sat.rtol, atol=sat.atol,
                          full_output = True )
    sat.answer   = answer
    sat.mid    = answer[sat.n_half]
    sat.mid    = answer[sat.n_half]
    sat.info     = info

if 1 == 1:
    xL2, xL1 = 0.74091, 1.21643  # lazy!
    fig = plt.figure(figsize=[10, 8])
    ax = fig.add_subplot(1, 1, 1, projection='3d')

    for sat in sats:
        x,  y,  z  = sat.answer.T[:3]
        ax.plot(x, y, z)

    ax.plot([0.0-mu], [0], [0], 'ob', markersize=20)
    ax.plot([1.0-mu], [0], [0], 'og', markersize=12)
    ax.plot([xL2], [0], [0], 'ok', markersize=8)
    ax.plot([xL1], [0], [0], 'ok', markersize=8)

    ax.set_xlim(0.7, 1.25)
    ax.set_ylim(-0.225, 0.225)
    ax.set_zlim(-0.15, 0.40)
    ax.text(xL1, 0, -0.05, "L1", fontsize=14, horizontalalignment='center')
    ax.text(xL2, 0, -0.05, "L2", fontsize=14, horizontalalignment='center')

    nplot    = 80
    thetas   = np.linspace(0, twopi, nplot+1)[:-1]
    azimuths = -90 + 10.0 * np.cos(thetas)

    fnames = []
    for i, azim in enumerate(azimuths):
        fname = "haloz_3D_" + str(10000+i)[1:]
        ax.elev, ax.azim = 0, azim
        plt.savefig(fname)
        fnames.append(fname)

    # tight cropping
    for i in range(len(fnames)):
        fname_in  = "haloz_3D_" + str(10000+i)[1:]
        fname_out = "haloz_3D_crop_" + str(10000+i)[1:] + ".png"
        img = plt.imread(fname_in + ".png")
        plt.imsave(fname_out, img[200:-175, 240:-190])

пользователь20022

Я знаю, что написание этого самого комментария лишает меня права на "Ага!" момент, но я боюсь, что я не экстрасенс ;-) Является ли ваша проблема объяснением

δ x_{f} = Φ (t_{f}, t_{0}) δ x_{0} + {\dot{x}}_{f} δ t

$\delta x_f = \Phi(t_f, t_0)\delta x_0 + \dot{x}_f \delta t$ ? Или это приложение, которое Хауэлл или Брейквелл и Браун делают из этого?

ооо

@LucJ.Bourhis, я знаю, это сложный вопрос. В эту секунду у меня нет под рукой B&B, но насколько я помню, это написано в стиле учебника, в том смысле, что оно ясное, методичное и продолжается от начала до конца, но не имело большого количества интуитивных помощников, те короткие пути, которые помогают некоторым людям некоторое время просто «получить это» все сразу. Что такое матрица переходов состояний? Матрица какого перехода состояния орбиты это? Когда у меня есть это чувство, мне намного легче пройти через вывод или применение. Ответ может состоять из нескольких предложений.

Ответы (3)

Как лучше всего представить Матрицу переходов состояний и как использовать ее для поиска периодических гало-орбит?

Я знаю, что написание этого самого комментария лишает меня права на "Ага!" момент, но я боюсь, что я не экстрасенс ;-) Является ли ваша проблема объяснением $\delta x_f = \Phi(t_f, t_0)\delta x_0 + \dot{x}_f \delta t$ ? Или это приложение, которое Хауэлл или Брейквелл и Браун делают из этого?
@LucJ.Bourhis, я знаю, это сложный вопрос. В эту секунду у меня нет под рукой B&B, но насколько я помню, это написано в стиле учебника, в том смысле, что оно ясное, методичное и продолжается от начала до конца, но не имело большого количества интуитивных помощников, те короткие пути, которые помогают некоторым людям некоторое время просто «получить это» все сразу. Что такое матрица переходов состояний? Матрица какого перехода состояния орбиты это? Когда у меня есть это чувство, мне намного легче пройти через вывод или применение. Ответ может состоять из нескольких предложений.

КрисР · Answer 1

Матрица перехода состояний (STM)

СТМ представляет собой процедуру линеаризации динамической системы. Он может использоваться для любой нелинейной динамической системы и используется для аппроксимации динамики системы за короткий период времени. В астродинамике он используется, в частности, для статистического определения орбиты (stat OD) и круговой ограниченной задачи третьего тела (CRTBP).

Вычисление STM для статистической OD подробно объясняется в «Статистическом определении орбиты» Tapley, Schultz, Born, Elsevier 2004. В частности, в разделах 1.2.5 и 4.2.1. С этого момента эта ссылка будет обозначаться как «(1)».

Системная динамика

Позволять $\boldsymbol{X}$ быть состоянием вашей системы в декартовой системе координат. В следующих, $\bf{r}$ а также $\bf{v}$ соответственно соответствуют положению и скорости КА; $\dot\gamma$ соответствует производной по времени $\gamma$ переменная. Выбор положения и скорости часто используется для задач начального уровня. Если вы делаете более серьезную статистическую OD, вы также захотите добавить гравитационный параметр, положение ваших наземных станций и т. д., но важно отметить, что изменение вашего вектора состояния также изменит STM и матрицу A (см. ниже).

Икс знак равно [\begin{matrix} р \\ в \end{matrix}] знак равно [\begin{matrix} Икс \\ у \\ г \\ \dot{Икс} \\ \dot{у} \\ \dot{г} \end{matrix}]

$\boldsymbol{X}=\left[\begin{array}{c} \boldsymbol{r}\\ \boldsymbol{v} \end{array}\right]=\left[\begin{array}{c} x\\ y\\ z\\ \dot{x}\\ \dot{y}\\ \dot{z} \end{array}\right]$

Затем мы можем выразить производную по времени состояния $\boldsymbol{X}$ следующим образом:

\dot{Икс} знак равно [\begin{matrix} \dot{р} \\ \dot{в} \end{matrix}] знак равно [\begin{matrix} \dot{Икс} \\ \dot{у} \\ \dot{г} \\ \ddot{Икс} \\ \ddot{у} \\ \ddot{г} \end{matrix}] знак равно Ф (Икс, т)

$\boldsymbol{\dot{X}}=\left[\begin{array}{c} \boldsymbol{\dot{r}}\\ \boldsymbol{\dot{v}} \end{array}\right]=\left[\begin{array}{c} \dot{x}\\ \dot{y}\\ \dot{z}\\ \ddot{x}\\ \ddot{y}\\ \ddot{z} \end{array}\right]=F\left(\boldsymbol{X}, t\right)$

В этой формулировке $F$ функция соответствует полной динамике системы: эта функция интегрируется по времени, если вы вычисляете реальную динамику, т.е. это представление уравнений движения. Предполагая задачу двух тел, $\boldsymbol{\dot{v}}$ есть ускорение, обусловленное только первичным телом, т.е. $-\frac{\mu}{r^3}\boldsymbol{X}$ . Если моделировать более сложную динамику, $F$ функция также будет включать их.

Назначение СТМ

Как было сказано выше, STM — это линеаризация вашей динамики. Итак, мы начнем с дискретизации времени и предположения , что система ведет себя линейно в течение этого времени. Это очень полезное приближение. На самом деле, это позволяет упростить моделирование: вместо того, чтобы распространять вашу динамику (т.е. $F$ функция) за заданное время интегрирования вам просто нужно умножить состояние $X_{i-1}$ с СТМ $\Phi$ чтобы получить $X_i$ . Кроме того, согласно (1), STM имеет следующие свойства (номер раздела и страницы указан в первой строке для справки):

Вычисление СТМ

Итак, на данный момент мы знаем, что СТМ — это линеаризация динамической системы, которая позволяет нам рассматривать ее как линейную систему в течение короткого периода времени. Итак, нам нужно линеаризовать динамику системы вокруг заданного состояния, здесь ссылка . Эта ссылка основывается на времени и обновляется через STM. Другими словами, мы вычисляем начальный STM, вычисляем состояние в следующий раз, а затем повторно вычисляем STM вокруг этого нового состояния.

Ниже приводится отрывок из лекции доктора МакМэхона. То, что отмечено звездочкой, соответствует эталонному состоянию.

Мы можем ясно видеть здесь, что мы просто вычисляем ряд Тейлора $F$ работать с первого раза! Итак, математически это просто. Однако на практике это соответствует производной ускорения, так что это немного утомительно для вычислений (но Mathematica или Sage Math (теперь CoCalc) могут помочь с их символьными производными, это может помочь ). Во всяком случае, это частичное обычно упоминается как $A$ матрица (по крайней мере, по моему опыту).

Связь между матрицей A и STM из «Анализ лагранжевой среды Солнце-Земля для New Worlds Observer (NWO)», Deccia 2017 ( ссылка )

Я думаю, что хорошим примером является рассмотрение того, как это можно сделать в коде (это из моей библиотеки астродинамики, которая находится на Голанге, извините... я думаю/надеюсь, что это все еще относительно читабельно). Во- первых, вычисление матрицы А с рядом возможных возмущений на основе конфигурации миссии. Во- вторых, серия тестовых случаев . Помимо прочего, тест проверяет, что норма разницы между предыдущим состоянием и новым состоянием (рассчитанная с помощью СТМ) находится в пределах $0.1$ (это несколько произвольно, но состояние имеет положение и скорость космического корабля LEO, так что это небольшая разница). В-третьих, вы можете проверить исходный код GMAT (который я сделал доступным на Github для удобства — проверьте их репозиторий sourceforge на наличие последних обновлений).

Гало-орбиты и STM

Судя по вашему вопросу, вы уже знаете орбиты гало, поэтому я не буду углубляться в них (я в любом случае не эксперт в них, поэтому могу сказать неправильные вещи). Короче говоря, гало вращается по квазипериодическим орбитам вокруг точек либрации (они периодические в CRTPB). Точки либрации — это точки равновесия между двумя массивными телами. По сути, орбита будет периодической в течение заданного времени. $T$ (и, следовательно, быть гало-орбитой) тогда и только тогда, когда на половине ее периода движение (т. е. скорость) космического корабля равно нулю во всех направлениях, кроме одного. Этот раздаточный материал доктора Дэвиса (из CCAR в CU Boulder) по поиску орбит гало из первоначального предположения подробно описывает, как это запрограммировать. Добавлю следующие уточнения:

Все вычисления выполняются после нормализации между обоими телами.
Это решает гало-орбиты только в круговой ограниченной задаче трех тел. В других условиях проблемы этот метод может не применяться как таковой или вообще не применяться.
$T/2$ соответствует полупериоду времени
STM интегрирован между нулевым временем и временем $T/2$ : это весь период дискретизации. (Если исходить из фона OD stat, это время намного больше, чем вы бы использовали).
Метод одиночной стрельбы позволяет найти орбиты, имеющие хотя бы один период. Гало-орбиты нестабильны по своей природе, поэтому вполне вероятно, что распространение «конечной» гало-орбиты заставляет ее расходиться после более чем одного витка (см. Рисунок ниже).

Отвечая на ваш вопрос (надеюсь)

Почему вы хотите использовать STM для поиска орбит Halo вместо того, чтобы перебирать все подряд?

Грубая сила редко бывает хорошей идеей. Он медленный, потому что ищет все возможные решения. Это полностью зависит от вашей дикретизации пространства решений. Представьте, что вы установили размер шага равным 0,5 в позиции нормализованного кадра, но на самом деле решение имеет приращение 0,2, тогда ваш метод никогда не сойдется.
СТМ позволяет выполнять несколько итераций, которые все ближе и ближе приближаются к орбите Гало. Вы должны ожидать, что алгоритм сойдется менее чем за 5-6 итераций (это ничто по сравнению с грубой силой).
Вы говорите о самом крутом спуске. Я считаю, что это будет включать метод градиентного спуска для поиска глобальных решений задач оптимизации. К СТМ можно применить градиентный спуск, но он не может работать с полной динамикой (система нелинейна). Кроме того, метод градиентного спуска применим к выпуклым задачам, но ваша задача не обязательно выпуклая (честно говоря, я не думаю, что это вообще): вы можете не найти решения. Таким образом, вам нужно найти выпуклую двойственную задачу и решить ее. Но переход к двойной задаче будет очень сложным, учитывая, что у вас нелинейная система. Наконец, что более важно, чем все вышеперечисленные математические вещи, какую функцию стоимости вы минимизируете? Где оптимальная задача?

Код?

Отказ от ответственности: я не проверял этот код Matlab. Он может содержать ошибки, иметь пограничные случаи, ломаться в определенных случаях и т. д. и т. д. Но это может помочь получить представление о том, как это реализовать: непроверенный код . (Я думаю, что я включил все файлы, необходимые для запуска этого, но если это не так, дайте мне знать в комментариях, и я добавлю его - у меня нет проблем с тем, чтобы поделиться своим кодом, как раз наоборот)

Я сосредоточил этот ответ в основном на STM, потому что мне показалось, что вы хорошо разбираетесь в орбитах Halo. Возможно, я забыл ряд важных «связей» между абзацами в моем ответе, которые убирают «Ага!» момент, который вы ищете. Если это так, пожалуйста, дайте мне знать здесь, и я улучшу его. Более того, я попросил друга, который исследует орбиты гало, проверить мой ответ, поэтому я мог бы внести некоторые изменения.
Ага! Я чувствовал, что грядет что-то хорошее. Хорошо, этот заслуживает свежего кофе, закопается завтра утром. Спасибо!!
Ваше первое предложение уже прояснило мое замешательство, я могу сказать, что это будет действительно хорошо, спасибо! Сегодня я « запутался » с другим ответом, но в эти выходные я испытаю его.
Там, где я упоминаю самый крутой (или градиентный) спуск, я на самом деле думаю только о том, чтобы начать с чего-то близкого и сойтись на гораздо более близком решении или «выйти» из одного ореола в соседний, чтобы найти семью. Это не лучший способ сделать это, но это будет лучше, чем грубая сила. Кроме того, используется полупериод $T/2$ просто экономия процессорного времени? Если бы кто-то расширился до более широкого разнообразия орбит точек либрации, имеющих непрерывную причудливую форму, он бы просто переключился на полный период, или в этом есть что-то особенное? $T/2$ ?
@uhoh, извините за задержку с ответом. Насколько я знаю, гало-орбиты не обязательно «близки» друг к другу. Даже с процедурой, изложенной в раздаточном материале Дэвиса, вам нужно иметь разумное предположение о действительной орбите Гало. Поэтому попытка угадать их, вероятно, займет некоторое время. Я подозреваю, но не уверен, что в некоторых публикациях есть целый список действительных гало-орбит. (1/н)
@uhoh, использование полупериода T/2 для физики орбиты. На полпути по орбите вы достигнете точки, где скорость космического корабля только в одном направлении по отношению к вращающейся раме обоих тел. Думайте об этом, как когда космический корабль находится в апогее в системе отсчета ECEF, где отсчет находится точно в надире: компоненты скорости меняют знак.
@uhoh, Так что да, в этом есть что-то особенное $T/2$ . Это время, когда вы ожидаете «увидеть» периодичность, и поэтому вы корректируете первоначальную оценку гало-орбиты (через СТМ), чтобы в это время она имела только одну ненулевую компоненту скорости. (3/n; n=3).
Я имею в виду $T$ вместо $T/2$ , и причина, по которой я спрашиваю, заключается в том, что существует целый зоопарк орбит CR3BP (рис. 40) , и, вероятно, некоторые из них не имеют подходящей симметрии для использования. $T/2$ (например, головастики или что-либо, связанное с треугольными точками либрации рис. 8 ), так что в этих случаях вам придется использовать $T$ . Поэтому мне любопытно, можно ли просто использовать $T$ везде, или если есть какое-то математическое преимущество помимо процессорного времени, которое $T/2$ обеспечивает более $T$ в случаях, когда вы знаете или хотите потребовать симметрии заранее.
Рисунки из Элементарные периодические орбиты, связанные с точками либрации в круговой ограниченной задаче трех тел , из EJ DOEDEL et al, Int. Дж. Раздвоенный хаос 17, 2625 (2007). doi.org/10.1142/S0218127407018671 (прямая ссылка worldscientific.com/doi/abs/10.1142/S0218127407018671 ). Или то, что я называю орбитальной бумагой для спагетти .
И то, что я имею в виду под «близко расположенными», будет, например, двумя соседними орбитами на одном из этих графиков спагетти, которые будут приращением вдоль непрерывного участка одной из линий, показанных на графиках бифуркации. Таким образом, вместо шести траекторий, показанных в моем вопросе, может быть шестьсот, более близко расположенных друг к другу. Если он у вас есть, вы можете найти смежное решение, произведя произвольное смещение, а затем повторно сойдясь. Вы не пойдете на старое решение, вы найдете новое, близкое соседнее решение. Суть моего вопроса!! "...создать семейство плотно расположенных гало-орбит".
Баунти снова голоден? Я добавил один в Как найти почти прямолинейные гало-орбиты (NHRO)? поиск ответов, чтобы дополнить простой ответ, который я оставил до сих пор. Возможно, простое объяснение линеаризации или как выбрать начальную точку?
@uhoh Я попрошу коллегу создать несколько графиков. В настоящее время мы делаем некоторую работу, финансируемую NRHO ;-)

пользователь20022 · Answer 2

Давайте попробуем! Для простоты я буду рассматривать одномерное уравнение движения

\begin{matrix} (1) & м \ddot{Икс (т)} знак равно а (т) Икс (т) + б (т) \dot{Икс} (т) \end{matrix}

$m \ddot{x(t)} = a(t) x(t) + b(t) \dot{x}(t) \tag{1}$

Приложения к гало-орбите на самом деле проще, потому что коэффициенты $a(t)$ а также $b(t)$ не будет зависеть от времени.

Теория линейных дифференциальных уравнений сообщает нам два важных результата:

Первоначальные условия $x(0)=x_0,\ \dot{x}(0)=\dot{x}_0$ полностью исправить решение;
Любая линейная комбинация двух решений также является решением.

Первый результат означает, что должна существовать функция, отображающая $(x_0,\dot{x}_0)$ на $x(t)$ . Второй результат гарантирует, что эта функция является линейной, т.е.

Икс (т) знак равно α (т) {Икс}_{0} + β (т) {\dot{Икс}}_{0}

$x(t) = \alpha(t)x_0 + \beta(t)\dot{x}_0$

Но тогда скорость имеет тот же вид

\dot{Икс} (т) знак равно \dot{α} (т) {Икс}_{0} + \dot{β} (т) {\dot{Икс}}_{0}

$\dot{x}(t) = \dot{\alpha}(t)x_0 + \dot{\beta}(t)\dot{x}_0$

и поэтому мы можем собрать все вместе

\begin{matrix} (2) & (\begin{matrix} Икс (т) \\ \dot{Икс} (т) \end{matrix}) знак равно \underset{Φ (т, т_{0})}{\underset{⏟}{(\begin{matrix} Φ_{11} (т, т_{0}) & Φ_{12} (т, т_{0}) \\ Φ_{21} (т, т_{0}) & Φ_{22} (т, т_{0}) \end{matrix})}} (\begin{matrix} {Икс}_{0} \\ {\dot{Икс}}_{0} \end{matrix}) \end{matrix}

$\begin{pmatrix} x(t) \\ \dot{x}(t) \end{pmatrix} = \underbrace{\begin{pmatrix} \Phi_{11}(t,t_0) & \Phi_{12}(t,t_0) \\ \Phi_{21}(t,t_0) & \Phi_{22}(t,t_0) \end{pmatrix}}_{\displaystyle\Phi(t,t_0)} \begin{pmatrix} x_0 \\ \dot{x}_0 \end{pmatrix} \tag{2}$

А также $\Phi(t,t_0)$ называется матрицей перехода от времени $t_0$ ко времени $t$ .

Из этого уравнения, поскольку $x(t)$ удовлетворяет дифференциальному уравнению (1), с которого мы начали, мы можем разумно ожидать $\Phi(t,t_0)$ чтобы удовлетворить один слишком. Чтобы найти его, нам просто нужно продифференцировать (2)

\begin{matrix} (3а) & (\begin{matrix} \dot{Икс} (т) \\ \ddot{Икс} (т) \end{matrix}) знак равно \dot{Φ} (т, т_{0}) (\begin{matrix} {Икс}_{0} \\ {\dot{Икс}}_{0} \end{matrix}) \end{matrix}

$\begin{pmatrix} \dot{x}(t) \\ \ddot{x}(t) \end{pmatrix} = \dot{\Phi}(t,t_0)\begin{pmatrix} x_0 \\ \dot{x}_0 \end{pmatrix}\tag{3a}$

куда $\dot{\Phi}(t,t_0)$ обозначает дифференцирование по $t$ , сохраняя $t_0$ постоянный. Но тогда левая сторона читается

(\begin{matrix} \dot{Икс} (т) \\ \ddot{Икс} (т) \end{matrix}) знак равно \underset{А (т)}{\underset{⏟}{(\begin{matrix} 0 & 1 \\ \frac{1}{м} а (т) & \frac{1}{м} б (т) \end{matrix})}} (\begin{matrix} Икс (т) \\ \dot{Икс} (т) \end{matrix})

$\begin{pmatrix} \dot{x}(t) \\ \ddot{x}(t) \end{pmatrix} = \underbrace{\begin{pmatrix} 0 & 1\\ \frac{1}{m}a(t) & \frac{1}{m}b(t) \end{pmatrix}}_{A(t)}\begin{pmatrix} x(t) \\ \dot{x}(t) \end{pmatrix}$ Затем мы используем (2) для замены

(\begin{matrix} x (t) \\ \dot{x} (t) \end{matrix})

$\begin{pmatrix} x(t) \\ \dot{x}(t) \end{pmatrix}$ с правой стороны.

\begin{matrix} (3б) & (\begin{matrix} \dot{Икс} (т) \\ \ddot{Икс} (т) \end{matrix}) знак равно А (т) Φ (т, т_{0}) (\begin{matrix} {Икс}_{0} \\ {\dot{Икс}}_{0} \end{matrix}) \end{matrix}

$\begin{pmatrix} \dot{x}(t) \\ \ddot{x}(t) \end{pmatrix} =A(t)\Phi(t,t_0)\begin{pmatrix} x_0 \\ \dot{x}_0 \end{pmatrix} \tag{3b}$

Приравнивая правые части (3а) и (3б), получаем

\dot{Φ} (т, т_{0}) (\begin{matrix} {Икс}_{0} \\ {\dot{Икс}}_{0} \end{matrix}) знак равно А (т) Φ (т, т_{0}) (\begin{matrix} {Икс}_{0} \\ {\dot{Икс}}_{0} \end{matrix})

$\dot{\Phi}(t,t_0)\begin{pmatrix} x_0 \\ \dot{x}_0 \end{pmatrix} = A(t)\Phi(t,t_0)\begin{pmatrix} x_0 \\ \dot{x}_0 \end{pmatrix}$

Это равенство должно выполняться для любого $x_0$ и любой $\dot{x}_0$ . Таким образом, матрицы, действующие на $\begin{pmatrix} x_0 \\ \dot{x}_0 \end{pmatrix}$ в обеих частях уравнения должны быть равны, и мы получаем искомое дифференциальное уравнение,

\begin{matrix} (4) & \dot{Φ} (т, т_{0}) знак равно А (т) Φ (т, т_{0}) . \end{matrix}

$\dot{\Phi}(t,t_0) = A(t)\Phi(t,t_0). \tag{4}$

Написав все это, я чувствую, что должен объяснить последний трюк в статье Хауэлла. Итак, у нас есть $x(t)$ и мы хотим понять, что может немного изменить его. $x(t)$ зависит от $t$ , такие разные $t$ по $\delta t$ индуцирует изменение, пропорциональное производной: $\dot{x}(t)\delta t$ . Но $x(t)$ также зависит от $x_0$ а также $\dot{x}_0$ и эта зависимость дается выражением (2). Вторая строка матрицы, если быть точным, и вариация $\Phi_{21}(t,t_0)\delta x_0 + \Phi_{22}(t,t_0)\delta \dot{x}_0$ . Затем, если мы рассмотрим только небольшие вариации, мы можем просто суммировать эти два вклада и получить:

дельта \dot{Икс} (т) знак равно Φ_{21} (т, т_{0}) дельта {Икс}_{0} + Φ_{22} (т, т_{0}) дельта {\dot{Икс}}_{0} + \dot{Икс} (т) дельта т

$\delta \dot{x}(t) = \Phi_{21}(t,t_0)\delta x_0 + \Phi_{22}(t,t_0)\delta \dot{x}_0 + \dot{x}(t)\delta t$

В интересующей вас проблеме, $t$ полупериод $T/2$ , а вариация $\delta \dot{x}(T/2)$ происходит либо от небольшой вариации $T/2$ , для тех же начальных условий, или от небольшого изменения начальных условий, для того же полупериода.

Я надеюсь, что это принесет некоторое просветление, и я желаю вам всего наилучшего для вашего проекта!

о боже, теперь я начал смотреть. Можете ли вы сказать мне, что $\dot{\mathbf{\Phi}}(t, t_0)$ в кавычках в моем вопросе имеется в виду? Это так же просто, как частная производная по отношению к $t$ и будет $\frac{\partial}{\partial t}$ применяется по месту к каждому элементу ${\mathbf{\Phi}(t, t_0})$ ?
@uhoh, к сожалению, это не так просто, как производные по времени. Уравнения движения не являются явными функциями времени, поэтому вы должны численно интегрировать их в свой код Python. Вы должны накапливать STM с течением времени, начиная с личности. Ключом является матрица F, представляющая собой частные производные уравнений движения.
@DuffBeerBaron, тогда что означает точка сверху $\mathbf{\Phi}$ значит если нет $\partial / \partial t$ поэлементно? Теперь я начинаю вспоминать, что именно эта точка остановила меня на моем пути.
В целом, STM — это «просто» выравнивание динамики за короткий период времени. Таким образом, вы берете части динамики ваших систем по отношению к каждой переменной, которую вы используете (обычно 3 для положения и еще 3 для скорости). Матрица А — это вторая производная вашего STM. Тэпли, Шульц, Борн довольно хорошо объясняют это в своей книге «Орбитальное определение».
В защиту Люка скажу, что проблема с STM заключается в том, что поначалу их трудно понять, но как только вы начинаете их использовать, они становятся очень простыми для понимания и, следовательно, их становится труднее объяснить. знак равно
@uhoh да, я извиняюсь за эту фразу. (Бывший) Академики невыносимы, не так ли? Я отредактировал свой ответ, чтобы попытаться лучше объяснить каждый шаг. Последний оставшийся трюк, который использовал Коннор Хауэлл, — это уравнение внизу страницы 56. После того, как оно объяснено, вы должны быть готовы к своему проекту.
Ниже уравнения 1 есть $x(t) = \alpha(t)x_0 + \beta(t)\dot{x}_0$ . Если $\alpha$ а также $\beta$ независимы от времени, как вы утверждаете, это просто $x(t) = \alpha x_0 + \beta \dot{x}_0$ Это правильно? Этот ответ «Как лучше всего придумать…» превращается в пошаговый вывод, где объясняются математические шаги. Я считаю, что это работает, поэтому не нужно убеждать, это не нужно доказывать.
« У меня нет B&B… но насколько я помню, это написано в стиле учебника, в том смысле, что оно ясное, методичное и продолжается от начала до конца, но не имеет большого количества интуитивных помощников, те короткие пути, которые помогают некоторым людям некоторое время просто "понимать" все сразу. Что такое матрица перехода состояния? Какого перехода состояния орбиты она представляет собой матрицу? Как только у меня появляется это чувство, мне гораздо проще пройтись по выведению или применению. Ответ мог бы состоять из нескольких предложений » .
Тем не менее (и бесстыдно и беспричинно повторенный), вывод и пошаговое объяснение здесь чрезвычайно полезны. Иногда для того, чтобы расколоть кокос, требуется атака со всех сторон. Спасибо!
@uhoh Да, нет проблем, на самом деле Хауэлл был профессором в моем старшекурснике. Раз уж я привлек ваше внимание, не могли бы вы объяснить мне последний «трюк»? Я не понимаю, почему есть термин 𝑥˙(𝑡)𝛿𝑡. Я думал, что STM сопоставляет состояние в момент t1 с состоянием в момент t2, так что не последний ли термин выполняет интеграцию по методу Эйлера поверх STM?
@Николай Вам повезло! :-) Я посмотрю, но не думаю, что смогу ответить на ваш вопрос, и я думаю, что есть еще несколько человек, которые могут. Я рекомендую вам просто опубликовать это как новый вопрос и позволить экспертам насладиться им.
Думаю, я понял это. Это похоже на прямой член Эйлера, но только внешне. 𝛿𝑡 действительно должно быть 𝛿(Точка). Конечное дельта-состояние может меняться не только в зависимости от начального положения, но и в зависимости от периода. Я сам не занимался математикой, но это может сработать, если вы добавите Period к вектору состояния и вычислите STM с этим 7-элементным вектором состояния. Вот статья одного из студентов Хауэлл, которая также объясняет это (многие студенты объясняли это в разных магистерских диссертациях) engineering.purdue.edu/people/kathleen.howell.1/Publications/…

ДаффПивоБарон · Answer 3

Сначала я постараюсь ответить на ваши два вопроса просто. Если эти ответы слишком просты или не соответствуют действительности, дайте мне знать, и я отредактирую ответ.

1) Что такое вектор распространения состояния и матрица перехода состояния (STM)?

Вектор распространения состояния — это просто положение и скорость в данный момент времени.

STM представляет собой матрицу, которая фиксирует чувствительность распространения к начальному состоянию. Таким образом, он отвечает на вопрос: «Если я изменю свою начальную координату x на 5 метров, насколько изменится мое конечное положение и скорость?»

2) Как я могу использовать STM для улучшения конвергенции на новых гало-орбитах?

Вы можете использовать СТМ для достижения более быстрой сходимости на новых гало-орбитах, отображая необходимые изменения на пересечении оси Y обратно в начальное состояние. (Например, если вы доберетесь до перекрестка со скоростью +2 по оси Z, вы можете использовать СТМ для расчета другого начального состояния, в котором скорость по оси Z будет уменьшена примерно на 2. (с учетом ошибок линеаризации) Д-р Дэвис из CU Boulder ( CCAR) предоставляет следующий раздаточный материал в курсе для выпускников «Проектирование межпланетных миссий», который она преподает:

http://ccar.colorado.edu/imd/2015/documents/SingleShootingHandout.pdf

Кроме того, вот краткое изложение проекта по гало-орбитам, которое включает ряд полезных цифр: http://ccar.colorado.edu/asen5050/projects/projects_2012/dowling/introduction.html

Классная ссылка, спасибо, посмотрю!
Однострочники особенно полезны, спасибо! Это также отличная ссылка; Я посмотрю.
@uhoh, я только что предложил еще один источник CCAR от профессора, который преподает IMD. Если это изменение не будет одобрено, вот ссылка: ccar.colorado.edu/imd/2015/documents/SingleShootingHandout.pdf . Предыдущая ссылка из студенческого проекта.
@uhoh, это соответствует разнице между состоянием, которое вы ожидаете, и тем, к которому приводит ваш пропагатор в момент времени T/2 (т.е. за полпериода). Из начального состояния в нормализованной вращающейся системе отсчета вы распространяете орбиту до тех пор, пока не достигнете ожидаемой точки пересечения (например, y=0). Остановите распространение (используя odeeventsв Matlab) и вычислите разницу. Затем выполните математику, чтобы найти поправку, необходимую для начального состояния.
@uhoh, это вычитание векторов между обоими состояниями. Состояние на половине периода и желаемое состояние (где одна составляющая положения равна нулю, а одна составляющая скорости равна нулю).
@uhoh, конечно, я могу сделать это за несколько часов и, в конце концов, также написать пример, похожий на тот, который Дэвис заставила класс сделать домашнюю работу несколько месяцев назад (я брал ее уроки).
@uhoh, если бы кто-то с историей редактирования мог это исправить, было бы здорово. Сейчас я разговариваю по телефону и собираюсь лечь спать, так что я исправлю это через несколько часов. Простите за это!
Ссылка на раздаточный материал по одиночной стрельбе не работает. Я получил файл и разместил его на STORJ . 230 тысяч загрузок, так что все должно быть в порядке.

Как лучше всего представить Матрицу переходов состояний и как использовать ее для поиска периодических гало-орбит?

ооо

пользователь20022

ооо

Ответы (3)

КрисР

Матрица перехода состояний (STM)

Системная динамика

Назначение СТМ

Вычисление СТМ

Гало-орбиты и STM

Отвечая на ваш вопрос (надеюсь)

Код?

КрисР

ооо

ооо

ооо

КрисР

КрисР

КрисР

ооо

ооо

ооо

ооо

КрисР

пользователь20022

ооо

ДаффПивоБарон

ооо

КрисР

КрисР

пользователь20022

ооо

ооо

ооо

Николай

ооо

Николай

ДаффПивоБарон

ооо

ооо

КрисР

КрисР

КрисР

КрисР

КрисР

Николай