Erklärung, wie moderne numerische Lösungen für die Umlaufbahngenerierung, Umlaufbahnplanung und optimale Steuerungsprobleme erhalten werden können.

[Ergänzung März 2017] Ich habe es als Python-Bibliothek als Erweiterung des hier geschriebenen Inhalts veröffentlicht. Bitte beachten Sie, dass sich die unten implementierte Theorie von der Legendre-Gauss-Lobatto-Pseudospektrum-Methode unterscheidet und der auf der aktuellen Seite geschriebene Inhalt sich von der Legendre-Gauss-Pseudo-Spektral-Methode unterscheidet. https://github.com/istellartech/OpenGoddard https://istellartech.github.io/OpenGoddard/

Verwendung von OpenGoddard 1 Verwendung von OpenGoddard 2 Verwendung von OpenGoddard 3 Verwendung von OpenGoddard 4

[Ergänzung hier]

Persönlich beabsichtige ich, die Grundlage für die Planung der Raketenumlaufbahn und die Erzeugung der Umlaufbahn zu sein.

Die Erzeugung der Umlaufbahn besteht beispielsweise darin, den Winkel zum Fliegen des Flugobjekts zum Zielort durch Minimieren des Treibstoffs zu ermitteln. Mathematisch kann es als optimales Steuerungsproblem definiert werden. In Raketenlehrbüchern wird eine Methode angewendet, die in der indirekten Methode als Variationsmethode bezeichnet wird. Dies dient dazu, das Problem der optimalen Steuerung in ein Hamilton-Randwertproblem umzuwandeln und es zu lösen.

Obwohl die Variationsmethode garantiert optimal ist, ist sie nicht vielseitig, da die Ausdruckstransformation von Problem zu Problem unterschiedlich ist und Variablen definiert, die nicht intuitiv sind. Eine Lösung für kluge Leute, die mit mathematischen Formeln spielen können. Es scheint, dass die Variantenmethode in der heutigen Zeit immer noch verwendet wird, aber jetzt, da reichlich Rechenressourcen verwendet werden können, ist eine Methode, die als direkte Methode bezeichnet wird, auch praktisch.

Die direkte Methode diskretisiert das ursprüngliche zeitlineare nichtlineare optimale Steuerungsproblem, wandelt es in ein nichtlineares Programmierproblem um und löst es. Unter den direkten Methoden gibt es zwei Typen: die Aufnahmemethode, bei der nur Steuervariablen unterschieden werden, und die direkte Kollokationsmethode, bei der Steuervariablen und Zustandsvariablen unterschieden werden. Darüber hinaus verwendet das Direktkollokationsverfahren auch ein Verfahren, das als Pseudospektralverfahren bezeichnet wird und die Dispersion nicht monoton teilt, so dass die Lösung selbst bei starken und nichtlinearen Problemen relativ stabil konvergiert. Die mathematische Grundlage der Pseudospektralmethode ist das Gaußsche Produkt, das ohne mathematische Ausbildung schwer zu verstehen ist (ich habe aufgegeben).

Hier wird das nichtlineare zeitoptimale Steuerungsproblem für die Umlaufbahnerzeugung unter Verwendung der Pseudospektralmethode, die eine Art direkte Kollokationsmethode ist, insbesondere der Legendre-Gauss-Pseudospektralmethode, und des Scipy-Optimierungsmoduls in ein nichtlineares Optimierungsproblem umgewandelt Ich habe es numerisch mit dem SLSQP-Solver gelöst.

Übrigens denke ich, dass der Grund, warum das Problem der Luft- und Raumfahrtoptimierung normalerweise das "Problem der eingeschränkten nichtlinearen Minimierung" ist, schwer zu verstehen ist. Ohne Einschränkungen, lineare oder quadratische Formate rollen unzählige nützliche Löser herum, aber mit "eingeschränkter Nichtlinearität" wird es schnell schwierig. Als Löser scheint die sequentielle quadratische Programmierung (SQP) in der heutigen Zeit effektiv zu sein.

Anstatt einfach einen Allzwecklöser für nichtlineare Optimierungsprobleme zu verwenden, scheint es einen handwerklichen Teil zu geben, der an große Probleme angepasst werden kann und durch verschiedene Einfallsreichtümer schneller konvergiert. Es gibt einige Beispiele für die Verwendung von Solvern wie IPOPT und SNOPT. Hier wird jedoch der Solver von scipy.optimize verwendet, da der Quellcode einfach ist und Priorität eingeräumt wird.

Während ich dies schrieb, habe ich versucht, die Aufnahmemethoden zu vergleichen, aber es konvergierte überhaupt nicht, nur Tsurai. Die Aufnahmemethode ist schwierig anzuwenden. Die Idee ist einfach, aber es ist besser, sie nicht zu übernehmen.

Beispiel

Lösen Sie als Beispiel ein einfaches.

Eine Rakete, die über eine zweidimensionale Ebene fliegt und weder Schwerkraft noch Luftwaffe empfängt, erzeugt eine Umlaufbahn mit dem minimalen Treibstoff (minimale Flugzeit), um die Zielhöhe und die Quergeschwindigkeit zu erhalten. Es gibt jedoch keine Zunahme oder Abnahme der Masse und die Beschleunigung ist konstant. Es ist eine kühne Vereinfachung. Zu diesem Zeitpunkt werden der Anfangswert und die Abbruchbedingung bestimmt, und die Abbruchzeit wird als Bewertungsfunktion (die schnellste Flugbahn) festgelegt. Das optimale Steuerungsproblem sind die folgenden Bewertungsfunktionen und Randbedingungen. (Eigentlich ist es sehr wichtig, ob es durch mathematische Formeln als Bewegungsgleichungen und optimale Steuerungsprobleme beschrieben werden kann)

\begin{aligned}
{\rm minimize} & : t_f \\
{\rm subject \ to} & : u(0) = v(0) = x(0) = y(0) = 0 \\
& : u(t_f) = 1,\  v(t_f) = 0, \ y(t_f) = 1 \\
& : -\frac{\pi}{2} \leq \beta \leq \frac{\pi}{2}
\end{aligned}

Die Bewegungsgleichung (Zustandsgleichung) ist wie folgt

\begin{aligned}
\dot{u} & = a \cos \beta \\
\dot{v} & = a \sin \beta \\
\dot{x} & = u \\
\dot{y} & = v
\end{aligned}

Der Quellcode ist unten. 5 Punkte

Es gibt eine Funktion zum Finden der diskreten Punkte (Legandre-Gauß-Knoten, Gauß-Punkt), Gewichte und Differentialmatrix.
Die zu optimierenden Variablen sind in einer einzigen Spalte angeordnet, wobei die Zustandsvariablen und Steuervariablen als diskrete Werte vorliegen.
Enthält Einschränkungen für die Übereinstimmung der numerischen Integration der Direct Collocation-Methode
Im numerischen Integrationsteil gibt es eine Randbedingung, ob das Ergebnis der trapezförmigen Integration vom Anfang bis zum Ende, so dass das System geschlossen ist, (letzter Wert) - (erster Wert) ist.
Normalisieren Sie in der Optimierungsberechnung die Zeitrichtung auf (-1,1)? Kehren Sie später zur ursprünglichen Zeitachse zurück

# -*- coding: utf-8 -*-
#Optimales Raketensteuerungsproblem, extrem vereinfachtes Problem
#Konvertieren Sie das Problem der optimalen Steuerung in ein nichtlineares Programmierproblem und suchen Sie nach einer Lösung mithilfe der sequentiellen Minimum-Square-Methode (SLSQP) von Scipy.
#Unter den direkten Methoden ist die Methode zum Ändern des Problems die direkte Kollokationsmethode und weiter
# Legendre-Es wird die Pseudospektralmethode von Gauß verwendet.
#
# Copyright (c) 2016 Takahiro Inagawa
# Released under the MIT license
import sys
reload(sys)
import platform
#Ändern Sie den Standardzeichencode.
sys.setdefaultencoding('utf-8')

import numpy as np
from scipy.integrate import odeint
from scipy.optimize import minimize
import matplotlib as mpl
import matplotlib.pyplot as plt
import matplotlib.font_manager
from matplotlib.font_manager import FontProperties
from matplotlib.backends.backend_pdf import PdfPages
import time

if 'Windows' == platform.system():
	font_path = r'C:\WINDOWS\Fonts\MSGothic.ttf'

if 'Darwin' == platform.system(): # for Mac
	font_path = '/Library/Fonts/Osaka.ttf'

font_prop = FontProperties(fname=font_path)
mpl.rcParams['font.family'] = font_prop.get_name()
plt.close('all')
plt.style.use('ggplot')
mpl.rcParams['axes.grid'] = True
mpl.rcParams['savefig.dpi'] = 300 #Standard-dpi=Wechsel von 100

# plt.ion()

def make_node_derivative_matrix(n):
	# Legendre-Aus der Gauß-Pseudospektralmethode als Kollokationspunkt
	# Legendre-Gauß-Knoten und Gewicht, Gauß-Differentialmatrix(Spektrale Differentialmatrix)Berechnung
	# tau: [-1,1]Legendre-Gauß-Knoten(collocation point)
	# w:Gewichte der Gaußschen numerischen Integrationsformel(n reelle Zahlen)
	# D:Gauß-Differentialmatrix
	beta = np.array([0.5 / np.sqrt(1-(2*i)**(-2)) for i in range(1,n)])
	T = np.diag(beta,1) + np.diag(beta,-1) #Jacobian
	D_, V = np.linalg.eig(T) #Einzigartige Wertzerlegung
	tau = np.sort(D_) # Legendre-Gauß-Knoten (Gauß-Punkt)
	i = np.argsort(D_)
	w = 2 * (V[0,i]**2) #Produktgewicht
	tau = np.hstack((-1,tau,1)) # -1,1(0,N+1)Hinzufügen
	D = np.zeros([n,n+1]) #Gauß-Differentialmatrix(Die Gaußsche Differentialmatrix gibt das Differential dieses Punktes an.)
	for k in range(1,n+1):
		for l in range(0,n+1):
			if k==l:
				D[k-1,l] = 0
				for m in range(0,n+1):
					if m != k:
						D[k-1,l] += 1.0/(tau[k]-tau[m])
			else:
				D[k-1,l] = 1.0/(tau[l]-tau[k])
				for m in range(0,n+1):
					if m != k and m != l:
						D[k-1,l] *= (tau[k]-tau[m])/(tau[l]-tau[m])
	return tau, w, D

def dynamics(p, div): #Bewegungsgleichung
	u = p[1:div[0]-1]
	v = p[div[0]+1:div[1]-1]
	x = p[div[1]+1:div[2]-1]
	y = p[div[2]+1:div[3]-1]
	beta = p[div[3]+1:div[4]-1]
	dx = np.zeros(0)
	a = 1.0 # acceleration
	dx0 = a * np.cos(beta)
	dx1 = a * np.sin(beta)
	dx2 = u
	dx3 = v
	dx = np.hstack((dx0, dx1, dx2, dx3))
	return dx

def evaluation(p, div): #Bewertungsfunktion
	return p[-1]

def equality(p, D, w, t0, div, n): #Gleichungsbeschränkungsbedingung
	u = p[0:div[0]-1]
	v = p[div[0]:div[1]-1]
	x = p[div[1]:div[2]-1]
	y = p[div[2]:div[3]-1]
	beta = p[div[3]:div[4]-1]
	dx = dynamics(p, div)
	result = np.zeros(0)
	derivative = np.hstack((D.dot(u), D.dot(v), D.dot(x), D.dot(y)))
	result = np.hstack((result, derivative - (p[-1]-t0) / 2.0 * dx))
	result = np.hstack((result, p[div[0]-1] - p[0] - np.sum(D.dot(u)*w)))
	result = np.hstack((result, p[div[1]-1] - p[div[0]] - np.sum(D.dot(v)*w)))
	result = np.hstack((result, p[div[2]-1] - p[div[1]] - np.sum(D.dot(x)*w)))
	result = np.hstack((result, p[div[3]-1] - p[div[2]] - np.sum(D.dot(y)*w)))
	result = np.hstack((result, u[0] - 0.0))
	# result = np.hstack((result, p[div[0]-1] - 1))
	result = np.hstack((result, v[0] - 0.0))
	result = np.hstack((result, p[div[1]-1] - 0))
	result = np.hstack((result, x[0] - 0.0))
	result = np.hstack((result, y[0] - 0.0))
	result = np.hstack((result, p[div[3]-1] - 1))
	return result

def inequality(p, D, w, t0, div, n): #Ungleichheitsbeschränkung
	u = p[0:div[0]-1]
	v = p[div[0]:div[1]-1]
	x = p[div[1]:div[2]-1]
	y = p[div[2]:div[3]-1]
	beta = p[div[3]:div[4]-1]
	dx = dynamics(p, div)
	result = np.zeros(0)
	result = np.hstack((result, beta + np.pi/2))
	result = np.hstack((result, np.pi/2 - beta))
	result = np.hstack((result, p[div[0]-1] - 1))
	return result

if __name__ == '__main__':
	t0 = 0 #Startzeit
	n = 20 #Diskriminiertes Netz (Legendre)-Anzahl der Gauß-Knoten)
	num_p = 5 #Gesamttypen von Zustands- und Steuervariablen
	iterator = 10 #Standardmäßig Anzahl der Optimierungsiterationen:1
	tau, w, D = make_node_derivative_matrix(n)
	div = [i*(n+2) for i in range(1,num_p+1)] #Index beim Übergang von Variablen
	num_variables = div[-1] + 1 #Gesamtzahl der Zustandsvariablen, Steuervariablen und statischen Variablen

	#Vektoren der oberen und unteren Grenze
	# vec_lower_limit = np.zeros(num_variables)
	# vec_upper_limit = np.ones(num_variables) * 100
	# vec_upper_limit[-1] = 0.5

	p = np.ones(num_variables) * 0.5 #Anfangswerte der zu optimierenden Variablen
	cons = ({'type': 'eq',
			 'fun': lambda p, D, w, t0, div, n: equality(p, D, w, t0, div, n),
			 'args': (D, w, t0, div, n,)},
			{'type': 'ineq',
			 'fun': lambda p, D, w, t0, div, n: inequality(p, D, w, t0, div, n),
			 'args': (D, w, t0, div, n,)})
	for i in range(iterator):
		opt = minimize(evaluation, p, args=(div,), constraints=cons,
		 			   method='SLSQP')
		p = opt.x
		print(u"iteration: %d" % (i+1))

	#Zerlegen Sie die resultierenden Variablen
	u = p[0:div[0]]
	v = p[div[0]:div[1]]
	x = p[div[1]:div[2]]
	y = p[div[2]:div[3]]
	beta = p[div[3]:div[4]]
	tf = p[-1]
	time = (tf-t0)/2.0 * tau + (tf+t0)/2.0

	#Ergebnis PLOT
	plt.figure()
	plt.subplot(3,1,1)
	plt.plot(time,u,label=u"vel x")
	plt.plot(time,v,label=u"vel y")
	plt.legend(loc="best")
	plt.ylabel(u"velocity [m/s]")
	plt.subplot(3,1,2)
	plt.plot(time,x,label=u"x")
	plt.plot(time,y,label=u"y")
	plt.legend(loc="best")
	plt.ylabel(u"position [m]")
	plt.subplot(3,1,3)
	plt.plot(time,beta,label=u"beta")
	plt.legend(loc="best")
	plt.ylabel(u"angle [rad]")
	plt.xlabel(u"time [s]")
	plt.show()

最適化の基礎.png

Verweise

Benson, David A., et al. "Direct trajectory optimization and costate estimation via an orthogonal collocation method." Journal of Guidance, Control, and Dynamics 29.6 (2006): 1435-1440.

https://goo.gl/BdLsvH　（PDF） ↑ Es wird empfohlen, das Papier des Originalmaterials der Methode zu lesen.

Koichiro Kato "Optimaler technischer Kontrollansatz zur Nichtlinearität" 1988 Da es die Variantenmethode verwendet, ist es nicht direkt verwandt, aber es scheint als ein Buch berühmt zu sein, das sich mit Problemen der optimalen Steuerung befasst.

Takeshi Tsuchiya und Shinji Suzuki. "Studie zur Problemlösung bei optimaler Steuerung unter Verwendung mathematischer Programmierung (Teil 2) Vorschlag der blockdiagonalen Cyan-Methode." Journal of Japan Aerospace Society 46.536 (1998): 497-503. https://www.jstage.jst.go.jp/article/jjsass1969/46/536/46_536_497/_article/-char/ja/

Nobuhiro Yokoyama, Shinji Suzuki und Takeshi Tsuchiya. "Studie zur Optimierung der aufsteigenden Flugbahn der Raumfläche zur ISS." JOURNAL DER JAPANISCHEN GESELLSCHAFT FÜR AERONAUTISCHE UND RAUMWISSENSCHAFTEN 49.570 (2001): 208-213. https://www.jstage.jst.go.jp/article/jjsass/49/570/49_570_208/_article/-char/ja/

↑ Es gibt mehrere Artikel von Suzuki und Tsuchiya Labs der Universität Tokio, die jedoch sehr hilfreich sind.

[PYTHON] Grundlagen der Umlaufbahngenerierung (Lösung eines nichtlinearen Problems der optimalen Steuerung mit Scipys Optimize)

Beispiel

Verweise