Wie stellt man sich die Zustandsübergangsmatrix am besten vor und wie verwendet man sie, um periodische Halo-Umlaufbahnen zu finden?

Question

Wie stellt man sich die Zustandsübergangsmatrix am besten vor und wie verwendet man sie, um periodische Halo-Umlaufbahnen zu finden?

äh

Ich werde zuerst meine mathematische Frage zu den Zustandsausbreitungs- und Zustandsübergangsmatrizen stellen und Ihnen dann ein einfaches Problem zeigen, für das ich diese Konzepte verwenden möchte, um eine dicht beabstandete Familie von Halo-Orbits zu erzeugen.

Ich werde auch mit der Aussage vorangehen, dass ich nach einem Aha suche ! Antwort eingeben. Ich hoffe nicht auf eine Erklärung, solange diese ausgezeichnete, intuitive Erklärung von Quaternionen . Ich brauche nicht alles ausgearbeitet, nur eine Erklärung, wie man die Zustandsübergangsmatrix in diesem Zusammenhang verstehen, erhalten und verwenden würde.

Das Folgende ist ziemlich normal, ich zitiere aus einem Papier, das ich gerade zur Hand habe, Juan Senent, Cesar Ocampo und Antonio Capella; Variable-spezifische Impulsübertragungen mit niedrigem Schub und Führung zu instabilen periodischen Umlaufbahnen. Journal of Guidance, Control, and Dynamics, 28. (2) März–April 2005:

Für das dynamische System

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

ausgewertet von $t_0=0$ für manchen $t=t_f$ , das Endzustandsdifferential bei $t_f$ wird von gegeben

$d x_{f} = Φ (t_{f}, t_{0}) δ 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$

wobei die Zustandsübergangsmatrix erfüllt

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

und

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

und $\mathbf{F}$ ist der Jacobi des Vektorfelds, das als Zustandsausbreitungsmatrix verwendet wird,

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

Ich habe mit dem klassischen Artikel von Kathleen Connor Howell Three-Dimensional, Periodic 'Halo' Orbits Celestial Mechanics 32 (1984) 53-71 begonnen. Es beschreibt eine Technik, um Lösungen für Halo-Umlaufbahnen im Circular Restricted 3-body Problem (CR3BP) zu finden, die sich eng an eine Technik anlehnt, die von Breakwell, JV und Brown, JV: 1979, The "Halo" Family of 3-Dimensional Periodic Orbits , beschrieben wurde im Erde-Mond-eingeschränkten 3-Körper-Problem Celest. Mech. 20 , 389.

Howell 1984 beschreibt im Detail ein schrittweises Verfahren, um Mitglieder einer Familie von Halo-Orbits um die kolinearen Lagrange-Librationspunkte zu finden, die Symmetrie um die xz-Ebene haben, indem die Tatsache ausgenutzt wird, dass für diese Gruppe von Orbits drei der sechs Komponenten des Zustandsvektors sollte an dem Punkt, an dem die Umlaufbahn die Ebene schneidet, gegen Null konvergieren.

Das Papier listet sechs Beispiele von Halo-Umlaufbahnen auf, und mit den dort angegebenen Zahlen kann ich die Zustandsvektoren integrieren und verifizieren, dass die drei Zustandsvektorkomponenten $y, v_x, v_z$ gehen tatsächlich am Mittelpunkt der Bahnen durch Null und machen ein schönes Diagramm.

Was ich tun möchte, ist intuitiv zu verstehen, was ein Zustandsausbreitungsvektor und ein Zustandsübergangsvektor sind und wie man diese verwendet, um schneller auf einem neuen Mitglied der Halo-Orbit-Familie zu konvergieren, als wenn ich gerade angefangen hätte, Orbits in einem Cluster zu schießen um einen Startpunkt herum und benutzte etwas Einfaches wie den steilsten Abstieg, um die nächste Umlaufbahn zu finden $y, v_x, v_z$ alle gleich Null.

\ddot{x} = x + 2 \dot{j} - \frac{(1 - μ) (x + μ)}{r_{1}^{3}} - \frac{μ (x - 1 + μ)}{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{j} = - 2 \dot{x} + j (1 - \frac{1 - μ}{r_{1}^{3}} - \frac{μ}{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 - μ}{r_{1}^{3}} + \frac{μ}{r_{2}^{3}})

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

wo

r_{1} = \sqrt{(x + μ)^{2} + j^{2} + z^{2}}

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

r_{2} = \sqrt{(x - 1 + μ)^{2} + j^{2} + z^{2}}

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

HINWEIS! Ich glaube, dass die Bezeichnungen für die Positionen von L ${}_1$ und ich ${}_2$ im GIF und Skript sind vertauscht (falsche Bezeichnungen/Namen). Ich werde das Bild bald aktualisieren.

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

Benutzer20022

Ich weiß, dass mich das Schreiben dieses Kommentars für ein „Aha!“ disqualifiziert. Momentan bin ich aber leider kein Hellseher ;-) Ist dein Problem die Erklärung

δ 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$ ? Oder ist es die Anwendung, die Howell oder Breakwell und Brown daraus machen?

äh

@LucJ.Bourhis ja, das ist eine schwierige Frage, die ich kenne. Ich habe B&B in dieser Sekunde nicht zur Hand, aber soweit ich mich erinnere, ist es im Lehrbuchstil geschrieben, insofern es klar und methodisch ist und von Anfang bis Ende vorgeht, aber nicht viele intuitive Helfer hatte, diese Abkürzungen, die einigen Leuten manchmal helfen, alles auf einmal zu "verstehen". Was ist eine Zustandsübergangsmatrix? Von welchem Übergang des Zustands einer Umlaufbahn ist sie die Matrix? Wenn ich dieses Gefühl habe, fällt es mir viel leichter, die Ableitung oder Anwendung durchzugehen. Die Antwort könnte aus ein paar Sätzen bestehen.

Antworten (3)

Wie stellt man sich die Zustandsübergangsmatrix am besten vor und wie verwendet man sie, um periodische Halo-Umlaufbahnen zu finden?

Ich weiß, dass mich das Schreiben dieses Kommentars für ein „Aha!“ disqualifiziert. Momentan bin ich aber leider kein Hellseher ;-) Ist dein Problem die Erklärung $\delta x_f = \Phi(t_f, t_0)\delta x_0 + \dot{x}_f \delta t$ ? Oder ist es die Anwendung, die Howell oder Breakwell und Brown daraus machen?
@LucJ.Bourhis ja, das ist eine schwierige Frage, die ich kenne. Ich habe B&B in dieser Sekunde nicht zur Hand, aber soweit ich mich erinnere, ist es im Lehrbuchstil geschrieben, insofern es klar und methodisch ist und von Anfang bis Ende vorgeht, aber nicht viele intuitive Helfer hatte, diese Abkürzungen, die einigen Leuten manchmal helfen, alles auf einmal zu "verstehen". Was ist eine Zustandsübergangsmatrix? Von welchem Übergang des Zustands einer Umlaufbahn ist sie die Matrix? Wenn ich dieses Gefühl habe, fällt es mir viel leichter, die Ableitung oder Anwendung durchzugehen. Die Antwort könnte aus ein paar Sätzen bestehen.

ChrisR · Answer 1

Die Zustandsübergangsmatrix (STM)

Das STM ist ein Linearisierungsverfahren eines dynamischen Systems. Es kann für jedes nichtlineare dynamische System verwendet werden und wird verwendet, um die Dynamik eines Systems über kurze Zeiträume zu approximieren . In der Astrodynamik wird es vor allem für die statistische Bahnbestimmung (stat OD) und das Circular Restricted Third Body Problem (CRTBP) verwendet.

Die Berechnung des STM für die statistische OD wird ausführlich in "Statistical Orbit Determination" von Tapley, Schultz, Born, Elsevier 2004 erläutert. Insbesondere Abschnitt 1.2.5 und 4.2.1. Von hier an wird diese Referenz als "(1)" bezeichnet.

Systemdynamik

Lassen $\boldsymbol{X}$ der Zustand Ihres Systems in einem kartesischen Rahmen sein. Im Folgenden, $\bf{r}$ und $\bf{v}$ jeweils der Position und der Geschwindigkeit des Raumfahrzeugs entsprechen; $\dot\gamma$ entspricht der zeitlichen Ableitung der $\gamma$ Variable. Die Wahl von Position und Geschwindigkeit ist oft das, was Sie für Einsteigerprobleme verwenden werden. Wenn Sie eine ernsthaftere statistische OD durchführen, sollten Sie auch den Gravitationsparameter, die Position Ihrer Bodenstationen usw. hinzufügen. Es ist jedoch wichtig zu beachten, dass eine Änderung Ihres Zustandsvektors auch das STM und die A-Matrix ändert (vgl. unter).

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

Wir können dann die zeitliche Ableitung des Zustands ausdrücken $\boldsymbol{X}$ folgendermaßen:

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

In dieser Formulierung ist die $F$ Funktion entspricht der vollen Dynamik des Systems: Diese Funktion wird über einen Zeitraum integriert, wenn Sie die reale Dynamik berechnen, dh sie ist eine Darstellung der Bewegungsgleichungen. Geht man vom Zweikörperproblem aus, $\boldsymbol{\dot{v}}$ ist die Beschleunigung nur auf den Primärkörper zurückzuführen, dh $-\frac{\mu}{r^3}\boldsymbol{X}$ . Bei der Modellierung komplexerer Dynamiken ist die $F$ Funktion wird auch diese enthalten.

Zweck des STM

Wie oben gesagt, ist das STM eine Linearisierung Ihrer Dynamik. Wir beginnen also mit der Diskretisierung der Zeit und der Annahme , dass sich das System während dieser Zeit linear verhält. Dies ist eine sehr nützliche Näherung. Tatsächlich ermöglicht es, die Simulation zu vereinfachen: Anstatt Ihre Dynamik (dh die $F$ Funktion) über eine gegebene Integrationszeit, müssen Sie einfach den Zustand multiplizieren $X_{i-1}$ mit dem STM $\Phi$ um zu bekommen $X_i$ . Darüber hinaus hat das STM gemäß (1) die folgenden Eigenschaften (Abschnitt und Seitenzahl in der ersten Zeile als Referenz angegeben):

Berechnung des STM

Bis jetzt wissen wir, dass das STM eine Linearisierung eines dynamischen Systems ist, was uns erlaubt, es über einen kurzen Zeitraum als ein lineares System zu betrachten. Wir müssen also die Dynamik des Systems um einen gegebenen Zustand herum linearisieren, hier Referenz . Diese Referenz basiert auf der Zeit und wird über das STM aktualisiert. Mit anderen Worten, wir berechnen das anfängliche STM, berechnen den Zustand beim nächsten Mal und berechnen dann das STM um diesen neuen Zustand herum neu.

Das Folgende ist ein Auszug aus einem Vortrag von Dr. McMahon. Was mit einem Stern markiert ist, entspricht dem Referenzzustand.

Wir können hier deutlich sehen, dass wir einfach die Taylor-Reihe von berechnen $F$ Funktion bei der ersten Bestellung! Also mathematisch ist das einfach. In der Praxis entspricht dies jedoch der Ableitung der Beschleunigung, daher ist die Berechnung etwas lästig (aber Mathematica oder Sage Math (jetzt CoCalc) können mit ihren symbolischen Ableitungen eine Menge helfen, dies könnte helfen ). Jedenfalls wird dieser Teil allgemein als der bezeichnet $A$ Matrix (zumindest meiner Erfahrung nach).

Beziehung zwischen einer Matrix und dem STM, aus "Analysis of the Sun-Earth Lagrangian environment for the New Worlds Observer (NWO)", Deccia 2017 ( Link )

Ich denke, ein gutes Beispiel ist, wie dies im Code gemacht werden kann (diese stammen aus meiner Astrodynamik-Bibliothek, die sich in Golang befindet, tut mir leid ... ich denke / hoffe, dass sie noch relativ lesbar ist). Erstens die Berechnung der A-Matrix mit einer Anzahl möglicher Störungen basierend auf der Missionskonfiguration. Zweitens eine Reihe von Testfällen . Der Test prüft unter anderem, ob die Norm der Differenz zwischen dem vorherigen Zustand und dem neuen Zustand (berechnet über das STM) innerhalb liegt $0.1$ (Dies ist etwas willkürlich, aber der Staat hat Positionen und Geschwindigkeiten eines LEO-Raumfahrzeugs, also ist dies ein winziger Unterschied). Drittens möchten Sie vielleicht die Codequelle von GMAT überprüfen (die ich der Einfachheit halber auf Github zur Verfügung gestellt habe – überprüfen Sie deren Sourceforge-Repository auf die neuesten Updates).

Halo Orbits und das STM

Aus Ihrer Frage geht hervor, dass Sie Halo-Umlaufbahnen bereits kennen, daher werde ich nicht auf diese eingehen (ich bin sowieso kein Experte für sie, daher könnte ich falsche Dinge sagen). Kurz gesagt, Halo umkreist quasi-periodische Umlaufbahnen um Librationspunkte (sie sind im CRTPB periodisch). Librationspunkte sind Gleichgewichtspunkte zwischen zwei massiven Körpern. Tatsächlich wird eine Umlaufbahn für eine gegebene Zeit periodisch sein $T$ (und daher eine Halo-Umlaufbahn sein), wenn und nur wenn bei der Hälfte seiner Periode die Bewegung (dh die Geschwindigkeit) des Raumfahrzeugs in allen Richtungen außer einer Null ist. Dieses Handout von Dr. Davis (von CCAR an der CU Boulder) zum Finden von Halo-Umlaufbahnen anhand einer anfänglichen Vermutung beschreibt, wie dies programmiert werden kann. Ich füge folgende Klarstellungen hinzu:

Alle Berechnungen werden nach einer Normalisierung zwischen beiden Körpern durchgeführt
Dies löst die Halo-Bahnen nur im kreisförmig eingeschränkten Drei-Körper-Problem. In anderen Problemsituationen kann diese Methode als solche nicht oder überhaupt nicht angewendet werden.
$T/2$ entspricht der Halbperiodenzeit
Der STM wird zwischen Zeit Null und Zeit integriert $T/2$ : Dies ist die gesamte Diskretisierungsperiode. (Wenn Sie von einem statistischen OD-Hintergrund kommen, ist diese Zeit viel viel größer als das, was Sie verwenden würden).
Die Single-Shooting-Methode ermöglicht es, Bahnen zu finden, die mindestens eine Periode haben. Halo-Umlaufbahnen sind von Natur aus instabil, daher ist es wahrscheinlich, dass die Fortpflanzung der "endgültigen" Halo-Umlaufbahn dazu führt, dass sie nach mehr als einer Umlaufbahn divergiert (vgl. Abbildung unten).

Beantwortung Ihrer Frage (hoffentlich)

Warum möchten Sie das STM verwenden, um Halo-Umlaufbahnen zu finden, anstatt alles brutal zu erzwingen?

Brute Force ist selten eine gute Idee. Es ist langsam, weil es nach allen möglichen Lösungen sucht. Es hängt ganz von Ihrer Dikretisierung des Lösungsraums ab. Stellen Sie sich vor, Sie setzen die Schrittweite auf 0,5 an der Position des normalisierten Frames, aber die Lösung liegt tatsächlich bei einem Inkrement von 0,2, dann wird Ihre Methode niemals konvergieren.
Das STM ermöglicht es, mehrere Iterationen durchzuführen, die der Halo-Umlaufbahn immer näher kommen. Sie sollten damit rechnen, dass der Algorithmus in weniger als 5-6 Iterationen konvergiert (das ist nichts im Vergleich zu einer Brute Force).
Sie beziehen sich auf einen steilsten Abstieg. Ich glaube, dies würde eine Gradientenabstiegsmethode beinhalten, um globale Lösungen für Optimierungsprobleme zu finden. Gradientenabstieg könnte auf das STM angewendet werden, aber es kann nicht mit der vollen Dynamik arbeiten (das System ist nicht linear). Darüber hinaus ist eine Gradientenabstiegsmethode auf konvexe Probleme anwendbar, aber Ihr Problem ist nicht unbedingt konvex (ich glaube nicht, dass es das ist, um ehrlich zu sein): Sie finden möglicherweise keine Lösung. Sie müssten also ein duales Problem finden, das konvex ist, und das duale Problem lösen. Die Umwandlung in das duale Problem wäre jedoch sehr kompliziert, da Sie ein nichtlineares System haben. Zu guter Letzt, und das ist noch wichtiger als der ganze mathematische Kram oben, was wäre die Kostenfunktion, die Sie minimieren? Wo ist das optimale Problem?

Code?

Haftungsausschluss: Ich habe diesen Matlab-Code nicht validiert. Es kann fehlerhaft sein, Grenzfälle haben, in bestimmten Fällen zusammenbrechen usw. usw. Aber es kann hilfreich sein, sich eine Vorstellung davon zu machen, wie dies implementiert werden kann: nicht validierter Code . (Ich glaube, ich habe alle Dateien hinzugefügt, die zum Ausführen erforderlich sind, aber wenn nicht, lassen Sie es mich in den Kommentaren wissen, und ich werde sie hinzufügen - ich habe kein Problem damit, meinen Code zu teilen, ganz im Gegenteil)

Ich habe diese Antwort hauptsächlich auf das STM konzentriert, weil es mir so vorkam, als hätten Sie ein solides Verständnis für Halo-Umlaufbahnen. Möglicherweise habe ich in meiner Antwort eine Reihe wichtiger "Verbindungen" zwischen den Absätzen vergessen, die das "Aha!" Moment, den Sie suchen. Wenn ja, lassen Sie es mich bitte hier wissen und ich werde es verbessern. Außerdem habe ich einen Freund, der Halo-Umlaufbahnen erforscht, gebeten, meine Antwort Korrektur zu lesen, damit ich dann möglicherweise einige Änderungen vornehmen kann.
Aha! Ich hatte das Gefühl, dass etwas Gutes kommen würde. OK, dieser verdient frischen Kaffee, wird morgen früh reingehauen. Vielen Dank!!
Dein erster Satz hat schon etwas von meiner Verwirrung beseitigt, ich kann sagen, dass das wirklich gut wird, danke! Heute habe ich mich mit einer anderen Antwort " verheddert ", aber ich werde dieses Wochenende einen Testflug geben.
Wo ich den steilsten (oder steilsten) Abstieg erwähne, denke ich wirklich nur daran, von etwas Nahem auszugehen und auf eine viel nähere Lösung zu konvergieren oder von einem Heiligenschein zu einem benachbarten "herauszugehen", um eine Familie zu finden. Es ist nicht der beste Weg, es zu tun, aber es wäre besser als rohe Gewalt. Außerdem wird die halbe Periode verwendet $T/2$ nur eine CPU-Zeitersparnis? Wenn man auf eine größere Vielfalt von immer seltsam geformten Librationspunkt-Umlaufbahnen ausdehnt, würde man einfach zu einer vollen Periode wechseln, oder gibt es da etwas Besonderes $T/2$ ?
@uhoh, Entschuldigung für die Verspätung meiner Antwort. AFAIK, Halo-Umlaufbahnen sind nicht unbedingt "nah" beieinander. Selbst mit dem in Davis' Handout beschriebenen Verfahren müssen Sie eine gültige Halo-Umlaufbahn vernünftig schätzen. Der Versuch, sie zu erraten, würde wahrscheinlich eine Weile dauern. Ich vermute, bin mir aber nicht sicher, dass einige Veröffentlichungen eine ganze Liste gültiger Halo-Umlaufbahnen haben. (1/n)
@uhoh, Die Verwendung der halben Periode von T / 2 dient der Physik der Umlaufbahn. Auf halbem Weg durch die Umlaufbahn erreichen Sie einen Punkt, an dem die Geschwindigkeit des Raumfahrzeugs in Bezug auf den rotierenden Rahmen beider Körper nur in einer Richtung ist. Stellen Sie sich das so vor, als ob sich ein Raumschiff in einem ECEF-Rahmen am Höhepunkt befindet, wo die Referenz genau am Nadir liegt: Die Geschwindigkeitskomponenten der Geschwindigkeit ändern das Vorzeichen.
@uhoh, also ja, da ist etwas Besonderes $T/2$ . Dies ist die Zeit, in der Sie erwarten, die Periodizität zu "sehen", und deshalb korrigieren Sie die anfängliche Schätzung der Halo-Umlaufbahn (über das STM), damit sie zu diesem Zeitpunkt nur eine Nicht-Null-Geschwindigkeitskomponente hat. (3/n; n=3).
Ich meine $T$ Anstatt von $T/2$ , und der Grund, warum ich frage, ist, dass es einen ganzen Zoo von CR3BP-Umlaufbahnen gibt (Abb. 40) und wahrscheinlich einige nicht die geeignete Symmetrie haben, um sie zu verwenden $T/2$ (z. B. Kaulquappen oder irgendetwas, das mit den dreieckigen Librationspunkten in Abb. 8 verbunden ist), so dass Sie in diesen Fällen verwenden müssten $T$ . Da bin ich mal gespannt, ob man das einfach verwenden könnte $T$ überall, oder ob es einen mathematischen Vorteil jenseits der CPU-Zeit gibt $T/2$ bietet über $T$ in Fällen, in denen Sie die Symmetrie im Voraus kennen oder benötigen möchten.
Abbildungen von elementaren periodischen Umlaufbahnen, die mit den Librationspunkten im kreisförmigen beschränkten 3-Körper-Problem verbunden sind, von EJ DOEDEL et al, Int. J. Bifurcation Chaos 17, 2625 (2007). doi.org/10.1142/S0218127407018671 (direkter Link worldscientific.com/doi/abs/10.1142/S0218127407018671 ). Oder was ich das orbitale Spaghettipapier nenne .
Und was ich mit "eng beieinander" meine, wären zum Beispiel zwei benachbarte Umlaufbahnen in einem dieser Spaghetti-Plots, die ein Inkrement entlang eines kontinuierlichen Abschnitts einer der in den Bifurkationsplots gezeigten Linien wären. Anstelle der in meiner Frage gezeigten sechs Trajektorien könnten es also sechshundert sein, die enger beieinander liegen. Wenn Sie eine haben, können Sie eine benachbarte Lösung finden, indem Sie willkürlich ausgleichen und dann wieder konvergieren. Sie werden nicht zur alten Lösung gehen, sondern eine neue, eng benachbarte Lösung finden. Der Kern meiner Frage!! "... um eine dicht beabstandete Familie von Halo-Umlaufbahnen zu erzeugen."
Schon wieder Kopfgeldhunger? Ich habe einen zu How to find Near Rectilinear Halo Orbits (NHROs) hinzugefügt? Ich suche nach Antworten, um die einfache Antwort zu ergänzen, die ich bisher hinterlassen habe. Vielleicht eine einfache Erklärung der Linearisierung oder wie man einen Startpunkt wählt?
@uhoh Ich werde einen Kollegen bitten, ein paar Diagramme zu erstellen. Wir machen derzeit einige finanzierte NRHO-Arbeiten ;-)

Benutzer20022 · Answer 2

Lass es uns versuchen! Der Einfachheit halber werde ich eine eindimensionale Bewegungsgleichung betrachten

\begin{matrix} (1) & m \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}$

Die Anwendungen in der Halo-Umlaufbahn sind aufgrund der Koeffizienten tatsächlich einfacher $a(t)$ und $b(t)$ würde nicht von der Zeit abhängen.

Die Theorie der linearen Differentialgleichungen zeigt uns zwei wichtige Ergebnisse:

Anfangsbedingungen $x(0)=x_0,\ \dot{x}(0)=\dot{x}_0$ die Lösung vollständig reparieren;
Jede Linearkombination zweier Lösungen ist auch eine Lösung.

Das erste Ergebnis impliziert, dass es eine abbildende Funktion geben muss $(x_0,\dot{x}_0)$ auf zu $x(t)$ . Das zweite Ergebnis garantiert, dass diese Funktion linear ist, dh

x (t) = a (t) x_{0} + β (t) {\dot{x}}_{0}

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

Aber dann hat die Geschwindigkeit die gleiche Form

\dot{x} (t) = \dot{a} (t) x_{0} + \dot{β} (t) {\dot{x}}_{0}

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

und wir können daher alles zusammenstellen

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

Und $\Phi(t,t_0)$ heißt die Übergangsmatrix von der Zeit $t_0$ zur Zeit $t$ .

Aus dieser Gleichung, da $x(t)$ die Differentialgleichung (1) erfüllt, von der wir ausgegangen sind, können wir vernünftigerweise erwarten $\Phi(t,t_0)$ auch einen zufrieden zu stellen. Um es zu finden, müssen wir nur differenzieren (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}$

wo $\dot{\Phi}(t,t_0)$ bezeichnet die Differenzierung nach $t$ , halten $t_0$ Konstante. Aber dann steht auf der linken Seite

(\begin{matrix} \dot{x} (t) \\ \ddot{x} (t) \end{matrix}) = \underset{EIN (t)}{\underset{⏟}{(\begin{matrix} 0 & 1 \\ \frac{1}{m} a (t) & \frac{1}{m} 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}$ Dann verwenden wir (2) zum Ersetzen

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

$\begin{pmatrix} x(t) \\ \dot{x}(t) \end{pmatrix}$ auf der rechten Seite.

\begin{matrix} (3b) & (\begin{matrix} \dot{x} (t) \\ \ddot{x} (t) \end{matrix}) = EIN (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}$

Durch Gleichsetzen der rechten Seite von (3a) und (3b) erhalten wir

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

Diese Gleichheit muss für alle gelten $x_0$ und alle $\dot{x}_0$ . Damit wirken die Matrizen auf $\begin{pmatrix} x_0 \\ \dot{x}_0 \end{pmatrix}$ auf beiden Seiten der Gleichung gleich sein, und wir erhalten die gesuchte Differentialgleichung,

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

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

Nachdem ich das alles geschrieben habe, habe ich das Gefühl, den letzten Trick in der Howell-Studie erklären zu müssen. Also haben wir $x(t)$ und wir wollen verstehen, was dazu führen könnte, dass es ein wenig variiert. $x(t)$ kommt drauf an $t$ , so unterschiedlich $t$ durch $\delta t$ induziert eine Variation proportional zur Ableitung: $\dot{x}(t)\delta t$ . Aber $x(t)$ hängt auch davon ab $x_0$ und $\dot{x}_0$ und diese Abhängigkeit ist durch (2) gegeben. Die zweite Zeile der Matrix, um genau zu sein, und die Variation ist $\Phi_{21}(t,t_0)\delta x_0 + \Phi_{22}(t,t_0)\delta \dot{x}_0$ . Wenn wir dann nur kleine Variationen betrachten, können wir diese beiden Beiträge einfach summieren und erhalten:

δ \dot{x} (t) = Φ_{21} (t, t_{0}) δ x_{0} + Φ_{22} (t, t_{0}) δ {\dot{x}}_{0} + \dot{x} (t) δ 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$

In dem für Sie interessanten Problem, $t$ ist die halbe Periode $T/2$ , und die Variation $\delta \dot{x}(T/2)$ stammt entweder aus einer kleinen Variation von $T/2$ , für die gleichen Anfangsbedingungen oder von einer kleinen Variation der Anfangsbedingungen für die gleiche Halbperiode.

Ich hoffe, es bringt etwas Erleuchtung und ich wünsche Ihnen das Beste für Ihr Projekt!

oh je, jetzt habe ich angefangen zu suchen. Können Sie mir sagen, was $\dot{\mathbf{\Phi}}(t, t_0)$ im Blockzitat in meiner Frage bedeutet? Ist es so einfach wie die partielle Ableitung bzgl $t$ und wäre $\frac{\partial}{\partial t}$ an Ort und Stelle auf jedes Element angewendet ${\mathbf{\Phi}(t, t_0})$ ?
@uhoh, leider ist es nicht so einfach wie Zeitableitungen. Die Bewegungsgleichungen sind keine expliziten Funktionen der Zeit, weshalb Sie sie numerisch in Ihren Python-Code integrieren müssen. Sie müssen den STM im Laufe der Zeit akkumulieren, beginnend mit der Identität. Der Schlüssel ist die F-Matrix, die die partiellen Ableitungen der Bewegungsgleichungen sind.
@DuffBeerBaron worauf steht dann der Punkt drauf $\mathbf{\Phi}$ meine wenn nicht $\partial / \partial t$ elementweise? Jetzt fange ich an, mich daran zu erinnern, dass es dieser Punkt war , der mich aufgehalten hat.
Insgesamt ist das STM „nur“ die Linearisierung der Dynamik über einen kurzen Zeitraum. Daher nehmen Sie die Teiltöne der Dynamik Ihrer Systeme in Bezug auf jede Variable, die Sie verwenden (normalerweise 3 für die Position und 3 weitere für die Geschwindigkeit). Die A-Matrix ist die zweite Ableitung Ihres STM. Tapley, Schultz, Born erklären das ziemlich gut in ihrem Buch Orbital Determination.
Um Luc zu verteidigen, das Problem mit STMs ist, dass sie anfangs schwer zu verstehen sind, aber sobald Sie anfangen, sie zu verwenden, werden sie sehr einfach zu verstehen und daher schwieriger zu erklären. =/
@uhoh ja, ich entschuldige mich für diese Formulierung. (Ehemalige) Akademiker sind unerträglich, oder? Ich habe meine Antwort bearbeitet, um zu versuchen, jeden Schritt besser zu erklären. Der letzte verbleibende Trick von Connor Howell ist die Gleichung unten auf Seite 56. Nachdem dies erklärt wurde, sollten Sie für Ihr Projekt bereit sein.
Unter Gleichung 1 gibt es $x(t) = \alpha(t)x_0 + \beta(t)\dot{x}_0$ . Wenn $\alpha$ und $\beta$ sind zeitunabhängig, wie Sie sagen, es ist nur $x(t) = \alpha x_0 + \beta \dot{x}_0$ Ist das korrekt? Diese "Wie man am besten an die ... denkt"-Antwort wird zu einer schrittweisen Herleitung, in der die mathematischen Schritte erklärt werden. Ich glaube, das funktioniert, muss also nicht überzeugt werden, es muss nicht bewiesen werden.
" Ich habe kein B&B... aber soweit ich mich erinnere, ist es im Lehrbuchstil geschrieben, insofern es klar und methodisch ist und von Anfang bis Ende vorgeht, aber nicht viele intuitive Helfer hatte, Diese Abkürzungen, die einigen Leuten manchmal helfen, alles auf einmal zu „verstehen". Was ist eine Zustandsübergangsmatrix? Von welchem Übergang eines Orbit-Zustands ist es die Matrix? Sobald ich das Gefühl habe, ist es das viel einfacher für mich, die Ableitung oder Anwendung durchzugehen. Die Antwort könnte aus ein paar Sätzen bestehen. “
Abgesehen davon (und schamlos und grundlos wiederholt) ist die Herleitung und Schritt-für-Schritt-Erklärung hier äußerst hilfreich. Manchmal erfordert das Knacken der Kokosnuss einen Angriff von allen Seiten. Vielen Dank!
@uhoh Ja, kein Problem, ich hatte Howell tatsächlich als Professor in meinem Grundstudium. Da ich Ihre Aufmerksamkeit habe, können Sie mir den letzten "Trick" erklären? Ich verstehe nicht, warum es einen 𝑥˙(𝑡)𝛿𝑡-Begriff gibt. Ich dachte, STM hat den Zustand bei t1 dem Zustand bei t2 zugeordnet, also führt dieser letzte Begriff nicht die Euler-Vorwärtsintegration auf dem STM durch?
@Nickolai Glück gehabt! :-) Ich werde es mir ansehen, aber ich glaube nicht, dass ich Ihre Frage beantworten kann, und ich denke, es gibt ein paar andere, die das können. Ich empfehle Ihnen, das einfach als neue Frage zu posten und die Experten daran teilhaben zu lassen.
Ich glaube, ich habe es herausgefunden. Es sieht aus wie ein Vorwärts-Euler-Term, aber nur oberflächlich. 𝛿𝑡 sollte wirklich 𝛿(Punkt) sein. Der endgültige Delta-Zustand kann nicht nur in Abhängigkeit von der Anfangsposition, sondern auch von der Periode variieren. Ich habe nicht selbst nachgerechnet, aber es könnte funktionieren, wenn Sie dem Zustandsvektor Period hinzufügen und das STM mit diesem 7-Elemente-Zustandsvektor berechnen. Hier ist ein Papier von einer von Howells Studenten, das es auch erklärt (sie hat viele Studenten dies in verschiedenen Masterarbeiten erklären lassen): engineering.purdue.edu/people/kathleen.howell.1/Publications/…

DuffBeerBaron · Answer 3

Ich werde versuchen, Ihre beiden Fragen zunächst einfach zu beantworten. Wenn diese Antworten zu einfach sind oder das Ziel verfehlen, lassen Sie es mich wissen, und ich werde die Antwort bearbeiten.

1) Was sind der Zustandsausbreitungsvektor und die Zustandsübergangsmatrix (STM)?

Der Zustandsausbreitungsvektor ist einfach die Position und Geschwindigkeit zu einem bestimmten Zeitpunkt.

Das STM ist eine Matrix, die die Empfindlichkeit der Ausbreitung gegenüber dem Anfangszustand erfasst. Es beantwortet also die Frage "Wenn ich meine Start-x-Koordinate um 5 Meter ändere, um wie viel ändern sich meine Endposition und Geschwindigkeit?"

2) Wie kann ich das STM verwenden, um die Konvergenz auf neuen Halo-Orbits zu verbessern?

Sie können das STM verwenden, um eine schnellere Konvergenz auf neuen Halo-Umlaufbahnen zu erreichen, indem Sie die Änderung, die Sie beim Kreuzen der Y-Achse benötigen, wieder auf den Ausgangszustand abbilden. (Wenn Sie z. B. mit einer Z-Geschwindigkeit von +2 an der Kreuzung ankommen, können Sie mit dem STM einen anderen Anfangszustand berechnen, bei dem die Z-Geschwindigkeit um etwa 2 verringert wird. (vorbehaltlich Linearisierungsfehlern) Dr. Davis von CU Boulder ( CCAR) bietet das folgende Handout in dem Graduiertenkurs Interplanetary Mission Design, den sie unterrichtet:

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

Darüber hinaus ist hier die Zusammenfassung eines Projekts zu Halo-Umlaufbahnen, die eine Reihe nützlicher Abbildungen enthält: http://ccar.colorado.edu/asen5050/projects/projects_2012/dowling/introduction.html

Die Einzeiler sind besonders hilfreich, danke! Das ist auch ein raffinierter Link; Ich werde es mir ansehen.
@uhoh, ich habe gerade eine weitere CCAR-Quelle von dem Professor vorgeschlagen, der IMD unterrichtet. Falls diese Änderung nicht genehmigt wird, hier ist der Link: ccar.colorado.edu/imd/2015/documents/SingleShootingHandout.pdf . Der vorherige Link stammt aus einem Studentenprojekt.
@uhoh, das entspricht der Differenz zwischen dem Zustand, den Sie erwarten, und dem, zu dem Ihr Propagator zum Zeitpunkt T / 2 (dh zu einer halben Periode) führt. Ausgehend von einem Anfangszustand in Ihrem normalisierten Rotationsrahmen propagieren Sie die Umlaufbahn, bis Sie den erwarteten Kreuzungspunkt (z. B. y=0) erreichen. Stoppen Sie die Ausbreitung ( odeeventsin Matlab) und berechnen Sie den Unterschied. Führen Sie dann die Mathematik durch, um die Korrektur zu finden, die für den Anfangszustand erforderlich ist.
@uhoh, es ist die Vektorsubtraktion zwischen beiden Zuständen. Der Zustand bei halber Periode und der gewünschte Zustand (wobei eine Komponente der Position Null ist und eine Komponente der Geschwindigkeit Null ist).
@uhoh, sicher, ich kann das in ein paar Stunden machen und schließlich auch ein Beispiel schreiben, das dem ähnelt, das Davis vor ein paar Monaten für die Klasse als Hausaufgabe gemacht hat (ich habe ihre Klasse belegt).
@uhoh, wenn jemand mit Bearbeitungsverlauf das reparieren könnte, wäre es großartig. Ich bin gerade auf meinem Handy und gehe gleich schlafen, also werde ich es in ein paar Stunden reparieren. Das tut mir leid!
Der Link für das Einzelschießen-Handout ist defekt. Ich habe die Datei erhalten und auf STORJ link.us1.storjshare.io/s/jvdnmmytv2ivraw4a47zpqfkzrfq/halo/… hochgeladen. Sie könnte brechen, sobald ich die Grenzen des kostenlosen Kontingents auf STORJ erreiche, aber nach meinen Berechnungen wird das ungefähr dauern 230.000 Downloads, also sollte es in Ordnung sein.

Wie stellt man sich die Zustandsübergangsmatrix am besten vor und wie verwendet man sie, um periodische Halo-Umlaufbahnen zu finden?

äh

Benutzer20022

äh

Antworten (3)

ChrisR

Die Zustandsübergangsmatrix (STM)

Systemdynamik

Zweck des STM

Berechnung des STM

Halo Orbits und das STM

Beantwortung Ihrer Frage (hoffentlich)

Code?

ChrisR

äh

äh

äh

ChrisR

ChrisR

ChrisR

äh

äh

äh

äh

ChrisR

Benutzer20022

äh

DuffBeerBaron

äh

ChrisR

ChrisR

Benutzer20022

äh

äh

äh

Nickolai

äh

Nickolai

DuffBeerBaron

äh

äh

ChrisR

ChrisR

ChrisR

ChrisR

ChrisR

Nickolai