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
ausgewertet von für manchen , das Endzustandsdifferential bei wird von gegeben
wobei die Zustandsübergangsmatrix erfüllt
und
und ist der Jacobi des Vektorfelds, das als Zustandsausbreitungsmatrix verwendet wird,
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 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 alle gleich Null.
wo
HINWEIS! Ich glaube, dass die Bezeichnungen für die Positionen von L und ich 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])
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.
Lassen der Zustand Ihres Systems in einem kartesischen Rahmen sein. Im Folgenden, und jeweils der Position und der Geschwindigkeit des Raumfahrzeugs entsprechen; entspricht der zeitlichen Ableitung der 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).
Wir können dann die zeitliche Ableitung des Zustands ausdrücken folgendermaßen:
In dieser Formulierung ist die 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, ist die Beschleunigung nur auf den Primärkörper zurückzuführen, dh . Bei der Modellierung komplexerer Dynamiken ist die Funktion wird auch diese enthalten.
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 Funktion) über eine gegebene Integrationszeit, müssen Sie einfach den Zustand multiplizieren mit dem STM um zu bekommen . Darüber hinaus hat das STM gemäß (1) die folgenden Eigenschaften (Abschnitt und Seitenzahl in der ersten Zeile als Referenz angegeben):
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 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 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 (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).
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 (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:
Warum möchten Sie das STM verwenden, um Halo-Umlaufbahnen zu finden, anstatt alles brutal zu erzwingen?
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)
Lass es uns versuchen! Der Einfachheit halber werde ich eine eindimensionale Bewegungsgleichung betrachten
Die Anwendungen in der Halo-Umlaufbahn sind aufgrund der Koeffizienten tatsächlich einfacher und würde nicht von der Zeit abhängen.
Die Theorie der linearen Differentialgleichungen zeigt uns zwei wichtige Ergebnisse:
Das erste Ergebnis impliziert, dass es eine abbildende Funktion geben muss auf zu . Das zweite Ergebnis garantiert, dass diese Funktion linear ist, dh
Aber dann hat die Geschwindigkeit die gleiche Form
und wir können daher alles zusammenstellen
Und heißt die Übergangsmatrix von der Zeit zur Zeit .
Aus dieser Gleichung, da die Differentialgleichung (1) erfüllt, von der wir ausgegangen sind, können wir vernünftigerweise erwarten auch einen zufrieden zu stellen. Um es zu finden, müssen wir nur differenzieren (2)
wo bezeichnet die Differenzierung nach , halten Konstante. Aber dann steht auf der linken Seite
Durch Gleichsetzen der rechten Seite von (3a) und (3b) erhalten wir
Diese Gleichheit muss für alle gelten und alle . Damit wirken die Matrizen auf auf beiden Seiten der Gleichung gleich sein, und wir erhalten die gesuchte Differentialgleichung,
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 und wir wollen verstehen, was dazu führen könnte, dass es ein wenig variiert. kommt drauf an , so unterschiedlich durch induziert eine Variation proportional zur Ableitung: . Aber hängt auch davon ab und und diese Abhängigkeit ist durch (2) gegeben. Die zweite Zeile der Matrix, um genau zu sein, und die Variation ist . Wenn wir dann nur kleine Variationen betrachten, können wir diese beiden Beiträge einfach summieren und erhalten:
In dem für Sie interessanten Problem, ist die halbe Periode , und die Variation stammt entweder aus einer kleinen Variation von , 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!
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
odeevents
in 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.
Benutzer20022
äh