¿Cómo pensar mejor en la Matriz de Transición de Estado y cómo usarla para encontrar órbitas periódicas de Halo?

Question

¿Cómo pensar mejor en la Matriz de Transición de Estado y cómo usarla para encontrar órbitas periódicas de Halo?

UH oh

Plantearé mi pregunta matemática sobre las matrices de transición de estado y propagación de estado primero, luego le mostraré un problema simple para el cual me gustaría usar estos conceptos para generar una familia densamente espaciada de órbitas de halo.

También comenzaré con la declaración de que estoy buscando un ¡Ajá! escriba respuesta. No espero una explicación tan larga como esta excelente e intuitiva explicación de los cuaterniones . No necesito que se resuelva todo, solo alguna explicación de cómo se podría entender, obtener y usar la Matriz de Transición del Estado en este contexto.

Lo siguiente es bastante estándar, estoy citando un artículo que tengo a mano en este momento, Juan Senent, César Ocampo y Antonio Capella; Transferencias de impulso específico variable de bajo empuje y guía para órbitas periódicas inestables. Journal of Guidance, Control, and Dynamics, 28 (2) de marzo a abril de 2005:

Para el sistema dinámico

$\dot{X} = F (X)$ $\mathbf{\dot x} = \mathbf{f}(\mathbf{x})$

evaluado desde $t_0=0$ Para algo $t=t_f$ , el diferencial de estado final en $t_f$ es dado por

$d X_{F} = Φ (t_{F}, t_{0}) d X_{0} + {\dot{X}}_{F} d t_{F}$ $\text{d} \mathbf{x}_f = \mathbf{\Phi}(t_f, t_0) \delta \mathbf{x}_0 + \mathbf{\dot x}_f \text{d} t_f$

donde la matriz de transición de estado satisface

$\dot{Φ} (t, t_{0}) = F (X (t)) Φ (t, t_{0})$ $\mathbf{\dot \Phi} (t,t_0) = \mathbf{F}(\mathbf{x}(t)) \mathbf{\Phi}(t, t_0)$

y

$Φ (t_{0}, t_{0}) = {yo}_{6 \times 6}$ $\mathbf{\Phi} (t_0, t_0) = \mathbf{I}_{6 \times 6}$

y $\mathbf{F}$ es el jacobiano del campo vectorial utilizado como matriz de propagación de estado,

$F (X (t)) = \frac{\partial F (X)}{\partial X}$ $\mathbf{F}(\mathbf{x}(t)) = \frac{\partial\mathbf{f}(\mathbf{x})}{\partial \mathbf{x}}$

Empecé con el artículo clásico escrito por Kathleen Connor Howell Three-Dimensional, Periodic 'Halo' Orbits Celestial Mechanics 32 (1984) 53-71. Describe una técnica para encontrar soluciones para las órbitas de halo en el problema circular restringido de 3 cuerpos (CR3BP), siguiendo de cerca una técnica descrita por Breakwell, JV y Brown, JV: 1979, The "Halo" Family of 3-Dimensional Periodic Orbits en el Problema de 3 Cuerpos Restringido Tierra-Luna Celest. mecánico 20 , 389.

Howell 1984 describe en detalle un procedimiento paso a paso para encontrar miembros de una familia de órbitas de halo alrededor de los puntos de libración colineales de Lagrange que tienen simetría alrededor del plano xz, aprovechando el hecho de que para este grupo de órbitas tres de los seis componentes del vector de estado debe converger a cero en el punto donde la órbita interseca al plano.

El documento tabula seis ejemplos de órbitas de halo, y con los números dados allí puedo integrar los vectores de estado, verificar que los tres componentes del vector de estado $y, v_x, v_z$ de hecho, pase por cero en el punto medio de las órbitas y haga un buen gráfico.

Lo que me gustaría hacer es entender intuitivamente qué es un vector de propagación de estado y un vector de transición de estado, y cómo usarlos para converger más rápido en un nuevo miembro de la familia de órbitas de halo que si hubiera comenzado a disparar órbitas en un cúmulo. alrededor de un punto de partida y usó algo simple como el descenso más pronunciado para encontrar la siguiente órbita con $y, v_x, v_z$ todos iguales a cero.

\ddot{X} = X + 2 \dot{y} - \frac{(1 - m) (X + m)}{r_{1}^{3}} - \frac{m (X - 1 + m)}{r_{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{y} = - 2 \dot{X} + y (1 - \frac{1 - m}{r_{1}^{3}} - \frac{m}{r_{2}^{3}})

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

\ddot{z} = - z (\frac{1 - m}{r_{1}^{3}} + \frac{m}{r_{2}^{3}})

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

dónde

r_{1} = \sqrt{(X + m)^{2} + y^{2} + z^{2}}

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

r_{2} = \sqrt{(X - 1 + m)^{2} + y^{2} + z^{2}}

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

¡NOTA! Creo que las etiquetas para las posiciones de L ${}_1$ y yo ${}_2$ en el GIF y el guión están transpuestos (etiquetas/nombres incorrectos). Actualizaré la imagen pronto.

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

usuario20022

Sé que escribir este mismo comentario me descalifica para un "¡Ajá!" momento pero me temo que no soy psíquico ;-) ¿Es su problema la explicación de

δ 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$ ? ¿O es la aplicación que hacen Howell, o Breakwell y Brown?

UH oh

@ LucJ.Bourhis ya es una pregunta difícil, lo sé. No tengo B&B a mano en este momento, pero por lo que recuerdo es que está escrito al estilo de un libro de texto, es claro, metódico y procede de principio a fin, pero no tenía muchos ayudantes intuitivos. esos atajos que ayudan a algunas personas a "entenderlo" todo a la vez. ¿Qué es una matriz de transición de estado? ¿De qué transición del estado de una órbita es la matriz? Una vez que tengo ese sentido, es mucho más fácil para mí pasar por la derivación o aplicación. La respuesta podría ser unas pocas frases.

Respuestas (3)

¿Cómo pensar mejor en la Matriz de Transición de Estado y cómo usarla para encontrar órbitas periódicas de Halo?

Sé que escribir este mismo comentario me descalifica para un "¡Ajá!" momento pero me temo que no soy psíquico ;-) ¿Es su problema la explicación de $\delta x_f = \Phi(t_f, t_0)\delta x_0 + \dot{x}_f \delta t$ ? ¿O es la aplicación que hacen Howell, o Breakwell y Brown?
@ LucJ.Bourhis ya es una pregunta difícil, lo sé. No tengo B&B a mano en este momento, pero por lo que recuerdo es que está escrito al estilo de un libro de texto, es claro, metódico y procede de principio a fin, pero no tenía muchos ayudantes intuitivos. esos atajos que ayudan a algunas personas a "entenderlo" todo a la vez. ¿Qué es una matriz de transición de estado? ¿De qué transición del estado de una órbita es la matriz? Una vez que tengo ese sentido, es mucho más fácil para mí pasar por la derivación o aplicación. La respuesta podría ser unas pocas frases.

CrisR · Answer 1

La Matriz de Transición del Estado (STM)

El STM es un procedimiento de linealización de un sistema dinámico. Se puede usar para cualquier sistema dinámico no lineal y se usa para aproximar la dinámica de un sistema en períodos cortos de tiempo. En astrodinámica, se utiliza especialmente para la determinación estadística de la órbita (stat OD) y el problema circular restringido del tercer cuerpo (CRTBP).

El cálculo del STM para la OD estadística se explica en profundidad en "Determinación estadística de la órbita" de Tapley, Schultz, Born, Elsevier 2004. Específicamente, las secciones 1.2.5 y 4.2.1. En adelante, esta referencia se denominará "(1)".

Sistemas dinámicos

Dejar $\boldsymbol{X}$ sea el estado de su sistema en un marco cartesiano. En el siguiente, $\bf{r}$ y $\bf{v}$ corresponden respectivamente a la posición y la velocidad de la nave espacial; $\dot\gamma$ corresponde a la derivada temporal de la $\gamma$ variable. Elegir la posición y la velocidad es a menudo lo que usará para los problemas de nivel de entrada. Si está haciendo una estadística OD más seria , también querrá agregar el parámetro gravitacional, la posición de sus estaciones terrestres, etc. pero es importante tener en cuenta que cambiar su vector de estado también cambiará el STM y la matriz A (cf. abajo).

X = [\begin{matrix} r \\ v \end{matrix}] = [\begin{matrix} X \\ y \\ z \\ \dot{X} \\ \dot{y} \\ \dot{z} \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]$

Entonces podemos expresar la derivada temporal del estado $\boldsymbol{X}$ como sigue:

\dot{X} = [\begin{matrix} \dot{r} \\ \dot{v} \end{matrix}] = [\begin{matrix} \dot{X} \\ \dot{y} \\ \dot{z} \\ \ddot{X} \\ \ddot{y} \\ \ddot{z} \end{matrix}] = F (X, t)

$\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)$

En esta formulación, el $F$ La función corresponde a la dinámica completa del sistema: esta función se integra durante un período de tiempo si está calculando la dinámica real, es decir, es una representación de las ecuaciones de movimiento. Suponiendo el problema de los dos cuerpos, $\boldsymbol{\dot{v}}$ es la aceleración debida únicamente al cuerpo primario, es decir $-\frac{\mu}{r^3}\boldsymbol{X}$ . Si se modelan dinámicas más complejas, el $F$ función también incluirá estos.

Propósito del STM

Como se dijo anteriormente, el STM es una linealización de su dinámica. Entonces comenzamos discretizando el tiempo y asumiendo que el sistema se comporta de manera lineal durante ese tiempo. Esta es una aproximación muy útil. De hecho, permite simplificar la simulación: en lugar de tener que propagar su dinámica (es decir, la $F$ función) durante un tiempo de integración dado, simplemente necesita multiplicar el estado $X_{i-1}$ con el STM $\Phi$ para obtener $X_i$ . Además, según (1), el STM tiene las siguientes propiedades (la sección y el número de página se muestran en la primera línea como referencia):

Cálculo del STM

Entonces, a partir de ahora, sabemos que el STM es una linealización de un sistema dinámico que nos permite considerarlo como un sistema lineal en un corto período de tiempo. Entonces, necesitamos linealizar la dinámica del sistema en torno a un estado dado, referencia aquí . Esta referencia se basa en el tiempo y se actualiza a través del STM. En otras palabras, calculamos el STM inicial, calculamos el estado la próxima vez y luego volvemos a calcular el STM alrededor de ese nuevo estado.

El siguiente es un extracto de una conferencia del Dr. McMahon. Lo que está marcado con una estrella corresponde al estado de referencia.

Podemos ver claramente aquí que simplemente estamos calculando la serie de Taylor de la $F$ función en el primer orden! Así que matemáticamente esto es simple. Sin embargo, en la práctica, esto corresponde a la derivada de la aceleración, por lo que es un poco molesto de calcular (pero Mathematica o Sage Math (ahora CoCalc) pueden ayudar mucho con sus derivadas simbólicas, esto podría ayudar ). De todos modos, este parcial generalmente se conoce como el $A$ matriz (al menos en mi experiencia).

Relación entre la matriz A y el STM, de "Análisis del entorno Lagrangiano Sol-Tierra para el New Worlds Observer (NWO)", Deccia 2017 ( enlace )

Creo que un buen ejemplo es ver cómo se puede hacer esto en el código (estos son de mi biblioteca de astrodinámica que está en Golang, lo siento... Creo/espero que todavía sea relativamente legible). En primer lugar, el cálculo de la matriz A con una serie de posibles perturbaciones en función de la configuración de la misión. En segundo lugar, una serie de casos de prueba . Entre otras cosas, la prueba verifica que la norma de la diferencia entre el estado anterior y el nuevo estado (calculado a través del STM) está dentro de $0.1$ (esto es algo arbitrario, pero el estado tiene posiciones y velocidades de una nave espacial LEO, por lo que esta es una pequeña diferencia). En tercer lugar, es posible que desee consultar el código fuente de GMAT (que he puesto a disposición en Github para su comodidad; consulte su repositorio sourceforge para obtener las últimas actualizaciones).

Halo Orbits y el STM

Según su pregunta, parece que ya conoce las órbitas de Halo, por lo que no me sumergiré en estas (de todos modos, no soy un experto en ellas, por lo que podría decir cosas incorrectas). En resumen, Halo orbita en órbitas cuasi-periódicas alrededor de los puntos de libración (son periódicas en el CRTPB). Los puntos de libración son puntos de equilibrio entre dos cuerpos masivos. En efecto, una órbita será periódica durante un tiempo dado. $T$ (y por lo tanto ser una órbita de Halo) si y sólo si en la mitad de su período, el movimiento (es decir, la velocidad) de la nave espacial es cero en todas las direcciones excepto en una. Este folleto del Dr. Davis (de CCAR en CU Boulder) sobre cómo encontrar órbitas de Halo a partir de una suposición inicial detalla cómo programar esto. Añadiré las siguientes aclaraciones:

Todos los cálculos se realizan después de una normalización entre ambos cuerpos.
Esto resuelve las órbitas de Halo solo en el problema circular restringido de tres cuerpos. En otras configuraciones de problemas, es posible que este método no se aplique como tal, o en absoluto.
$T/2$ corresponde al tiempo de la mitad del período
El STM está integrado entre el tiempo cero y el tiempo $T/2$ : esto es todo el período de discretización. (Si proviene de un fondo de OD estadístico, este tiempo es mucho más grande de lo que usaría).
El método de disparo único permite encontrar órbitas que tienen al menos un período. Las órbitas de Halo son inestables por naturaleza, por lo que es probable que la propagación de la órbita de Halo "final" haga que diverja después de más de una órbita (ver figura a continuación).

Respondiendo a tu pregunta (con suerte)

¿Por qué quieres usar el STM para encontrar órbitas de Halo en lugar de forzarlo todo?

La fuerza bruta rara vez es una buena idea. Es lento porque busca todas las soluciones posibles. Depende enteramente de su dicretización del espacio de solución. Imagine que establece el tamaño del paso en 0,5 en la posición del marco normalizado, pero la solución es en realidad un incremento de 0,2, entonces su método nunca convergerá.
El STM permite realizar varias iteraciones que se acercan cada vez más a la órbita de Halo. Debería esperar que el algoritmo converja en menos de 5-6 iteraciones (eso no es nada comparado con una fuerza bruta).
Te refieres a un descenso más empinado. Creo que esto implicaría un método de descenso de gradiente para encontrar soluciones globales a los problemas de optimización. El descenso de gradiente podría aplicarse al STM, pero no puede funcionar con la dinámica completa (el sistema no es lineal). Además, un método de descenso de gradiente es aplicable a problemas convexos, pero su problema no es necesariamente convexo (para ser honesto, no creo que lo sea): es posible que no encuentre una solución. Así que tendrías que encontrar un problema dual que sea convexo y resolver el problema dual. Pero convertirlo al problema dual sería muy complicado dado que tienes un sistema no lineal. Finalmente, y más importante que todas las matemáticas anteriores, ¿cuál sería la función de costo que estás minimizando? ¿Dónde está el problema óptimo?

¿Código?

Descargo de responsabilidad: no he validado este código de Matlab. Puede tener errores, tener casos límite, descomponerse en casos específicos, etc., etc. Pero, puede ser útil tener una idea de cómo implementar esto: código no validado . (Creo que he incluido todos los archivos necesarios para ejecutar esto, pero si no lo he hecho, házmelo saber en los comentarios y lo agregaré; no tengo ningún problema en compartir mi código, todo lo contrario)

Enfoqué esta respuesta principalmente en el STM porque me pareció que tienes una comprensión sólida de las órbitas de Halo. Es posible que haya olvidado una serie de "enlaces" importantes entre los párrafos de mi respuesta que eliminan el "¡Ajá!" momento que buscas. Si es así, házmelo saber aquí y lo mejoraré. Además, le pedí a un amigo que está investigando las órbitas de Halo que lea mi respuesta, por lo que podría hacer algunos cambios en ese momento.
¡Ajá! Tenía la sensación de que algo bueno se avecinaba. OK, este merece café recién hecho, lo cavará mañana por la mañana. ¡¡Gracias!!
Su primera oración ya ha aclarado algo de mi confusión, puedo decir que esto va a ser realmente bueno, ¡gracias! Hoy me he " enredado " con otra respuesta, pero le daré un vuelo de prueba este fin de semana.
Cuando menciono el descenso más pronunciado (o gradiente), en realidad solo estoy pensando en comenzar desde algo cercano y converger en una solución mucho más cercana, o "salir" de un halo a uno adyacente para encontrar una familia. No es la mejor manera de hacerlo, pero sería mejor que la fuerza bruta. Además, está usando el medio período $T/2$ ¿solo un ahorro de tiempo de CPU? Si uno se expandiera a una variedad más amplia de órbitas de puntos de libración cada vez más extrañas, ¿cambiaría simplemente a un período completo, o hay algo especial en $T/2$ ?
@uhoh, perdón por la demora en mi respuesta. AFAIK, las órbitas de Halo no están necesariamente "cercanas" entre sí. Incluso con el procedimiento descrito en el folleto de Davis, debe tener una estimación razonable de una órbita de Halo válida. Así que tratar de adivinarlos probablemente tomaría un tiempo. Sospecho, pero no estoy seguro, que algunas publicaciones tienen una lista completa de órbitas de Halo válidas. (1/n)
@uhoh, usar el medio período de T/2 es para la física de la órbita. A la mitad de la órbita, llegarás a un punto donde la velocidad de la nave espacial es solo en una dirección con respecto al marco giratorio de ambos cuerpos. Piense en ello como cuando una nave espacial está en el apogeo en un marco ECEF donde la referencia está exactamente en el nadir: los componentes de velocidad de la velocidad cambiarán de signo.
@uhoh, Entonces sí, hay algo especial en $T/2$ . Este es el momento en el que espera "ver" la periodicidad y, por lo tanto, corrige la suposición inicial de la órbita del halo (a través del STM) para que tenga solo un componente de velocidad no nulo en ese momento. (3/n; n=3).
quiero decir $T$ en vez de $T/2$ , y la razón por la que pregunto es porque hay todo un zoológico de órbitas CR3BP (Fig. 40) y probablemente algunas no tengan la simetría adecuada para usar $T/2$ (por ejemplo, renacuajos o cualquier cosa asociada con los puntos de libración triangular Fig 8 ) para que en esos casos tendría que usar $T$ . Así que tengo curiosidad si uno podría usar $T$ en todas partes, o si hay algún beneficio matemático más allá del tiempo de CPU que $T/2$ proporciona sobre $T$ en los casos en que sepa o quiera requerir la simetría con anticipación.
Cifras de órbitas periódicas elementales asociadas con los puntos de libración en el problema circular restringido de 3 cuerpos , de EJ DOEDEL et al, Int. J. Bifurcación Caos 17, 2625 (2007). doi.org/10.1142/S0218127407018671 (enlace directo worldscientific.com/doi/abs/10.1142/S0218127407018671 ). O lo que yo llamo el papel de espagueti orbital .
Y lo que quiero decir con "estrechamente espaciados" sería, por ejemplo, dos órbitas adyacentes en una de esas gráficas de espagueti, lo que sería un incremento a lo largo de una sección continua de una de las líneas que se muestran en las gráficas de bifurcación. Entonces, en lugar de las seis trayectorias que se muestran en mi pregunta, podría haber seiscientas, más juntas. Si tiene uno, puede encontrar una solución adyacente compensando arbitrariamente y luego volviendo a converger. No irá a la solución anterior, encontrará una nueva solución adyacente poco espaciada. ¡¡El quid de mi pregunta!! "... para generar una familia densamente espaciada de órbitas de halo".
Bounty hambriento de nuevo todavía? He agregado uno a ¿Cómo encontrar órbitas de halo casi rectilíneas (NHRO)? buscando respuestas para complementar la respuesta básica que he dejado hasta ahora. ¿Quizás una explicación simple de la linealización, o cómo elegir un punto de partida?
@uhoh, le pediré a un compañero de trabajo que genere algunos gráficos. Actualmente estamos haciendo un trabajo financiado por NRHO ;-)

usuario20022 · Answer 2

¡Vamos a intentarlo! Para mantenerlo simple, consideraré una ecuación de movimiento unidimensional

\begin{matrix} (1) & metro \ddot{X (t)} = a (t) X (t) + b (t) \dot{X} (t) \end{matrix}

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

Las aplicaciones a la órbita del halo son en realidad más simples porque los coeficientes $a(t)$ y $b(t)$ no dependería del tiempo.

La teoría de las ecuaciones diferenciales lineales nos dice dos resultados importantes:

Condiciones iniciales $x(0)=x_0,\ \dot{x}(0)=\dot{x}_0$ arreglar completamente la solución;
Cualquier combinación lineal de dos soluciones también es una solución.

El primer resultado implica que debe existir una función que mapea $(x_0,\dot{x}_0)$ sobre $x(t)$ . El segundo resultado garantiza que esta función es lineal, es decir

X (t) = α (t) X_{0} + β (t) {\dot{X}}_{0}

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

Pero entonces la velocidad tiene la misma forma.

\dot{X} (t) = \dot{α} (t) X_{0} + \dot{β} (t) {\dot{X}}_{0}

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

y por lo tanto podemos poner todo junto

\begin{matrix} (2) & (\begin{matrix} X (t) \\ \dot{X} (t) \end{matrix}) = \underset{Φ (t, t_{0})}{\underset{⏟}{(\begin{matrix} Φ_{11} (t, t_{0}) & Φ_{12} (t, t_{0}) \\ Φ_{21} (t, t_{0}) & Φ_{22} (t, t_{0}) \end{matrix})}} (\begin{matrix} X_{0} \\ {\dot{X}}_{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}$

Y $\Phi(t,t_0)$ se llama la matriz de transición del tiempo $t_0$ al tiempo $t$ .

De esta ecuación, ya que $x(t)$ satisface la ecuación diferencial (1) de la que partimos, podemos esperar razonablemente $\Phi(t,t_0)$ para satisfacer a uno también. Para encontrarlo, solo necesitamos diferenciar (2)

\begin{matrix} (3a) & (\begin{matrix} \dot{X} (t) \\ \ddot{X} (t) \end{matrix}) = \dot{Φ} (t, t_{0}) (\begin{matrix} X_{0} \\ {\dot{X}}_{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}$

dónde $\dot{\Phi}(t,t_0)$ denota la diferenciación con respecto a $t$ , manteniendo $t_0$ constante. Pero luego el lado izquierdo dice

(\begin{matrix} \dot{X} (t) \\ \ddot{X} (t) \end{matrix}) = \underset{A (t)}{\underset{⏟}{(\begin{matrix} 0 & 1 \\ \frac{1}{metro} a (t) & \frac{1}{metro} b (t) \end{matrix})}} (\begin{matrix} X (t) \\ \dot{X} (t) \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}$ Entonces usamos (2) para reemplazar

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

$\begin{pmatrix} x(t) \\ \dot{x}(t) \end{pmatrix}$ al lado derecho.

\begin{matrix} (3b) & (\begin{matrix} \dot{X} (t) \\ \ddot{X} (t) \end{matrix}) = A (t) Φ (t, t_{0}) (\begin{matrix} X_{0} \\ {\dot{X}}_{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}$

Igualando el lado derecho de (3a) y (3b), obtenemos

\dot{Φ} (t, t_{0}) (\begin{matrix} X_{0} \\ {\dot{X}}_{0} \end{matrix}) = A (t) Φ (t, t_{0}) (\begin{matrix} X_{0} \\ {\dot{X}}_{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}$

Esta igualdad debe ser cierta para cualquier $x_0$ y cualquier $\dot{x}_0$ . Así, las matrices que actúan sobre $\begin{pmatrix} x_0 \\ \dot{x}_0 \end{pmatrix}$ en ambos lados de la ecuación será igual, y obtenemos la ecuación diferencial que buscamos,

\begin{matrix} (4) & \dot{Φ} (t, t_{0}) = A (t) Φ (t, t_{0}) . \end{matrix}

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

Después de escribir todo eso, siento que tengo que explicar el último truco del artículo de Howell. Entonces tenemos $x(t)$ y queremos entender qué podría hacer que varíe un poquito. $x(t)$ depende de $t$ , tan variable $t$ por $\delta t$ induce una variación, proporcional a la derivada: $\dot{x}(t)\delta t$ . Pero $x(t)$ también depende de $x_0$ y $\dot{x}_0$ y esa dependencia viene dada por (2). La segunda fila de la matriz para ser exactos, y la variación es $\Phi_{21}(t,t_0)\delta x_0 + \Phi_{22}(t,t_0)\delta \dot{x}_0$ . Entonces, si consideramos solo pequeñas variaciones, podemos simplemente sumar esas dos contribuciones y obtener:

d \dot{X} (t) = Φ_{21} (t, t_{0}) d X_{0} + Φ_{22} (t, t_{0}) d {\dot{X}}_{0} + \dot{X} (t) d t

$\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$

En el problema de su interés, $t$ es el medio periodo $T/2$ , y la variación $\delta \dot{x}(T/2)$ proviene de una pequeña variación de $T/2$ , para las mismas condiciones iniciales, o de una pequeña variación de las condiciones iniciales, para el mismo semiperíodo.

¡Espero que traiga algo de iluminación y te deseo lo mejor para tu proyecto!

oh caray ahora he empezado a buscar. Me puede decir qué $\dot{\mathbf{\Phi}}(t, t_0)$ en la cita en bloque en mi pregunta significa? ¿Es tan simple como la derivada parcial con respecto a $t$ y sería $\frac{\partial}{\partial t}$ aplicado en su lugar a cada elemento de ${\mathbf{\Phi}(t, t_0})$ ?
@uhoh, desafortunadamente, no es tan simple como las derivadas del tiempo. Las ecuaciones de movimiento no son funciones explícitas de tiempo, por lo que debe integrarlas numéricamente en su código Python. Debe acumular el STM a lo largo del tiempo, a partir de la identidad. La clave es la matriz F, siendo las derivadas parciales de las ecuaciones de movimiento.
@DuffBeerBaron entonces, ¿qué significa el punto en la parte superior? $\mathbf{\Phi}$ significa si no $\partial / \partial t$ ¿elemento sabio? Estoy empezando a recordar ahora que fue ese punto lo que me detuvo en seco.
En general, el STM es "simplemente" la linerización de la dinámica en un corto período de tiempo. Por lo tanto, estás tomando los parciales de la dinámica de tus sistemas con respecto a cada variable que estás usando (generalmente 3 para posición y 3 más para velocidad). La matriz A es la segunda derivada de su STM. Tapley, Schultz, Born lo explican bastante bien en su libro Orbital Determination.
En defensa de Luc, el problema con los STM es que son difíciles de comprender al principio, pero una vez que comienzas a usarlos, se vuelven muy fáciles de entender y, por lo tanto, se vuelven más difíciles de explicar. =/
@uhoh sí, me disculpo por esa frase. Los (ex) académicos son insufribles, ¿no? Edité mi respuesta para tratar de explicar mejor cada paso. El último truco restante utilizado por Connor Howell es la ecuación que se encuentra en la parte inferior de la página 56. Después de explicar esto, debe estar listo para su proyecto.
Debajo de la ecuación 1, hay $x(t) = \alpha(t)x_0 + \beta(t)\dot{x}_0$ . Si $\alpha$ y $\beta$ son independientes del tiempo como dices, es solo $x(t) = \alpha x_0 + \beta \dot{x}_0$ ¿Es eso correcto? Esta respuesta de "Cómo pensar mejor en..." se está convirtiendo en una derivación paso a paso, donde se explican los pasos matemáticos. Creo que esto funciona, así que no tienes que estar convencido, no tiene que ser probado.
" No tengo B&B... pero por lo que recuerdo es que está escrito al estilo de un libro de texto, en el sentido de que es claro, metódico y procede de principio a fin, pero no tenía muchos asistentes intuitivos, esos atajos que ayudan a algunas personas algunas veces a "captarlo" todo de una vez. ¿Qué es una matriz de transición de estado? ¿De qué transición del estado de una órbita es la matriz? Una vez que tengo ese sentido, es mucho más fácil para mí pasar por la derivación o la aplicación. La respuesta podría ser unas pocas oraciones " .
Dicho esto (y repetido descaradamente y gratuitamente), la derivación y la explicación paso a paso aquí es extremadamente útil. A veces, romper el coco requiere un asalto por todos lados. ¡Gracias!
@uhoh Sí, no hay problema, en realidad tuve a Howell como profesor en mi licenciatura. Ya que tengo tu atención, ¿puedes explicarme el último "truco"? No entiendo por qué hay un término 𝑥˙(𝑡)𝛿𝑡. Pensé que STM asignó el estado en t1 al estado en t2, entonces, ¿no es ese último término hacer una integración directa de euler sobre el STM?
@Nickolai ¡Suerte! :-) Echaré un vistazo, pero no creo que pueda responder a su pregunta y creo que hay algunos otros que pueden hacerlo. Le recomiendo que simplemente publique eso como una nueva pregunta y deje que los expertos disfruten haciéndolo.
Creo que lo descubrí. Parece un término de Euler directo, pero solo superficialmente. 𝛿𝑡 realmente debería ser 𝛿(Punto). El estado delta final puede variar no solo en función de la posición inicial, sino también del período. No he hecho los cálculos yo mismo, pero podría funcionar si agrega Período al vector de estado y calcula el STM con ese vector de estado de 7 elementos. Aquí hay un artículo de uno de los estudiantes de Howell que también lo explica (muchos estudiantes han explicado esto en diferentes tesis de maestría) engineering.purdue.edu/people/kathleen.howell.1/Publications/…

DuffCervezaBaron · Answer 3

Trataré de responder a sus dos preguntas simplemente primero. Si estas respuestas son demasiado simples o no dan en el blanco, házmelo saber y editaré la respuesta.

1) ¿Qué son el vector de propagación de estado y la matriz de transición de estado (STM)?

El vector de propagación de estado es simplemente la posición y la velocidad en un momento dado.

El STM es una matriz que captura la sensibilidad de la propagación al estado inicial. Entonces, responde a la pregunta "Si cambio mi coordenada x inicial en 5 metros, ¿cuánto cambiarán mi posición y velocidad finales?"

2) ¿Cómo puedo usar el STM para mejorar la convergencia en las nuevas órbitas de Halo?

Puede usar el STM para lograr una convergencia más rápida en las nuevas órbitas de Halo mapeando el cambio que necesita en el eje Y que cruza de regreso al estado inicial. (Por ejemplo, si llega al cruce con una velocidad Z de +2, puede usar el STM para calcular un estado inicial diferente que tendrá una velocidad Z reducida en aproximadamente 2. (sujeto a errores de linealización) Dr. Davis de CU Boulder ( CCAR) proporciona el siguiente folleto en el curso de posgrado de Diseño de misiones interplanetarias que imparte:

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

Además, aquí está el resumen de un proyecto sobre las órbitas de Halo que incluye una serie de figuras útiles: http://ccar.colorado.edu/asen5050/projects/projects_2012/dowling/introduction.html

Ese es un enlace ingenioso, ¡gracias, le echaré un vistazo!
Las frases ingeniosas son particularmente útiles, ¡gracias! Ese también es un enlace ingenioso; Le echaré un vistazo.
@uhoh, acabo de proponer otra fuente CCAR del profesor que enseña IMD. En caso de que no se apruebe este cambio, aquí está el enlace: ccar.colorado.edu/imd/2015/documents/SingleShootingHandout.pdf . El enlace anterior es de un proyecto de estudiante.
@uhoh, eso corresponde a la diferencia entre el estado que espera y el que conduce su propagador en el tiempo T/2 (es decir, en medio período). Desde un estado inicial en su marco giratorio normalizado, propaga la órbita hasta que alcanza el punto de cruce esperado (por ejemplo, y=0). Detenga la propagación (usando odeeventsen Matlab) y calcule la diferencia. Luego haga los cálculos para encontrar la corrección necesaria en el estado inicial.
@uhoh, es la resta de vectores entre ambos estados. El estado en la mitad del período y el estado deseado (donde un componente de la posición es cero y un componente de la velocidad es cero).
@uhoh, claro, puedo hacer eso en unas pocas horas y eventualmente también escribir un ejemplo similar al que Davis hizo que la clase hiciera como tarea hace unos meses (tomé su clase).
@uhoh, si alguien con historial de edición pudiera reparar eso, sería genial. Estoy en mi teléfono actualmente y a punto de irme a dormir, así que lo arreglaré en unas horas. ¡Lo siento por eso!
El enlace para el folleto de tiro individual está roto. Obtuve el archivo y lo puse en STORJ link.us1.storjshare.io/s/jvdnmmytv2ivraw4a47zpqfkzrfq/halo/… Podría romperse una vez que alcance los límites en el nivel gratuito en STORJ, pero según mis cálculos, eso tomará aproximadamente 230k descargas por lo que debería estar bien.

¿Cómo pensar mejor en la Matriz de Transición de Estado y cómo usarla para encontrar órbitas periódicas de Halo?

UH oh

usuario20022

UH oh

Respuestas (3)

CrisR

La Matriz de Transición del Estado (STM)

Sistemas dinámicos

Propósito del STM

Cálculo del STM

Halo Orbits y el STM

Respondiendo a tu pregunta (con suerte)

¿Código?

CrisR

UH oh

UH oh

UH oh

CrisR

CrisR

CrisR

UH oh

UH oh

UH oh

UH oh

CrisR

usuario20022

UH oh

DuffCervezaBaron

UH oh

CrisR

CrisR

usuario20022

UH oh

UH oh

UH oh

Nickolai

UH oh

Nickolai

DuffCervezaBaron

UH oh

UH oh

CrisR

CrisR

CrisR

CrisR

CrisR

Nickolai