Explication de la façon d'obtenir des solutions numériques modernes pour la génération d'orbite, la planification d'orbite et les problèmes de contrôle optimal.

[Ajout mars 2017] Je l'ai publié en tant que bibliothèque python en tant qu'extension du contenu écrit ici. Veuillez noter que la théorie implémentée ci-dessous est différente de la méthode pseudo-spectrale Legendre-Gauss-Lobatto, et le contenu écrit sur la page courante est différent de la méthode Legendre-Gauss-Pseudo-spectrale. https://github.com/istellartech/OpenGoddard https://istellartech.github.io/OpenGoddard/

Comment utiliser OpenGoddard 1 Comment utiliser OpenGoddard 2 Comment utiliser OpenGoddard 3 Ouvrir Goddard Comment utiliser 4

[Ajout ici]

Personnellement, j'ai l'intention d'être la base de la planification de l'orbite des fusées et de la génération des orbites.

La génération d'orbite consiste, par exemple, à trouver l'angle de vol vers l'emplacement cible en minimisant le carburant du projectile. Mathématiquement, il peut être défini comme un problème de contrôle optimal. Dans les manuels de fusée, une méthode appelée méthode de variation dans la méthode indirecte est adoptée. Il s'agit de convertir le problème de contrôle optimal en un problème de valeur aux limites hamiltonien et de le résoudre.

Bien que la méthode de variation soit garantie optimale, elle n'est pas polyvalente car la transformation de l'expression diffère d'un problème à un autre et définit des variables qui ne sont pas intuitives. Une solution pour les personnes intelligentes qui peuvent jouer avec des formules mathématiques. Il semble que la méthode des variantes soit encore utilisée dans les temps modernes, mais maintenant que des ressources informatiques abondantes peuvent être utilisées, une méthode appelée méthode directe est également pratique.

La méthode directe discrétise le problème original de contrôle optimal non linéaire en temps continu, le convertit en un problème de programmation non linéaire et le résout. Parmi les méthodes directes, il existe deux types: la méthode de prise de vue, qui ne discrimine que les variables de contrôle, et la méthode de collocation directe, qui discrimine les variables de contrôle et les variables d'état. En outre, la méthode de collocation directe adopte également une méthode appelée méthode pseudo-spectrale qui ne divise pas de manière monotone la dispersion, de sorte que la solution converge de manière relativement stable même pour des problèmes forts et non linéaires. La base mathématique de la méthode pseudo-spectrale est le produit gaussien, difficile à comprendre sans éducation mathématique (j'ai abandonné).

Ici, le problème de contrôle optimal en temps continu non linéaire pour la génération d'orbite est converti en un problème d'optimisation non linéaire à l'aide de la méthode pseudo-spectrale, qui est une sorte de méthode de collocation directe, en particulier la méthode pseudo-spectrale Legendre-Gauss, et le module d'optimisation Scipy Je l'ai résolu numériquement en utilisant le solveur SLSQP à l'intérieur.

En passant, je pense que la raison pour laquelle le problème d'optimisation aérospatiale est généralement le "problème de minimisation non linéaire contrainte" est difficile à comprendre. En l'absence de contraintes, de formats linéaires ou quadratiques, il existe des tonnes de solveurs utiles, mais avec une "non-linéarité contrainte", cela devient rapidement difficile. En tant que solveur, il semble que la programmation quadratique séquentielle (SQP) soit efficace à l'époque moderne.

Au lieu d'utiliser simplement un solveur polyvalent pour les problèmes d'optimisation non linéaires, il semble qu'il existe une partie artisanale qui peut être adaptée à des problèmes à grande échelle et converge plus rapidement en faisant diverses ingéniosité. Il existe quelques exemples d'utilisation de solveurs tels que IPOPT et SNOPT, mais ici, le solveur de scipy.optimize est utilisé car le code source est simple et la priorité est donnée.

En écrivant ceci, j'ai essayé de comparer les méthodes de prise de vue, mais cela n'a pas du tout convergé, seulement Tsurai. La méthode de prise de vue est difficile à utiliser. L'idée est simple, mais il vaut mieux ne pas l'adopter.

Exemple

Par exemple, résolvez-en une simple.

Une fusée qui survole un plan bidimensionnel et ne reçoit ni gravité ni force aérienne crée une orbite avec le minimum de carburant (temps de vol minimum) afin d'obtenir l'altitude cible et la vitesse latérale. Cependant, il n'y a pas d'augmentation ou de diminution de masse et l'accélération est constante. C'est une simplification audacieuse. A ce moment, la valeur initiale et la condition de fin sont déterminées, et le temps de fin est défini comme fonction d'évaluation (la trajectoire la plus rapide). Le problème de contrôle optimal est la fonction d'évaluation et les conditions de contrainte suivantes. (En fait, il est assez important de savoir s'il peut être décrit par des formules mathématiques comme des équations de mouvement et des problèmes de contrôle optimal)

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

L'équation de mouvement (équation d'état) est la suivante

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

Le code source est ci-dessous. 5 points

--Il existe une fonction pour trouver les points discrets (nœud de Legandre-Gauss, point de Gauss), les poids et la matrice différentielle. --Les variables à optimiser sont organisées dans une seule colonne avec les variables d'état et les variables de contrôle sous forme de valeurs discrètes.

Contient des contraintes pour l'intégration numérique correspondante de la méthode Direct Collocation --Dans la partie d'intégration numérique, il y a une condition de contrainte si le résultat de l'intégration trapézoïdale du début à la fin pour que le système soit fermé soit (dernière valeur) - (première valeur). --Dans le calcul d'optimisation, normaliser la direction du temps à (-1,1)? Et puis revenez à l'axe du temps d'origine plus tard

# -*- coding: utf-8 -*-
#Problème de contrôle optimal des fusées, problème extrêmement simplifié
#Convertissez le problème de contrôle optimal en problème de programmation non linéaire et recherchez une solution à l'aide de la méthode des carrés minimum séquentiels de Scipy (SLSQP).
#Parmi les méthodes directes, la méthode pour changer le problème est la méthode Direct Collocation, et en outre
# Legendre-La méthode pseudo-spectrale de Gauss est utilisée.
#
# Copyright (c) 2016 Takahiro Inagawa
# Released under the MIT license
import sys
reload(sys)
import platform
#Modifiez le code de caractère par défaut.
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 #dpi par défaut=Changement de 100

# plt.ion()

def make_node_derivative_matrix(n):
	# Legendre-De la méthode pseudo-spectrale de Gauss, comme point de collocation
	# Legendre-Nœud et poids de Gauss, matrice différentielle de Gauss(Matrice différentielle spectrale)Calculer
	# tau: [-1,1]Legendre-Nœud de Gauss(collocation point)
	# w:Poids de formule d'intégration numérique gaussienne(n nombres réels)
	# D:Matrice différentielle de Gauss
	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) #Jacobien
	D_, V = np.linalg.eig(T) #Décomposition de valeur unique
	tau = np.sort(D_) # Legendre-Nœud de Gauss (point de Gauss)
	i = np.argsort(D_)
	w = 2 * (V[0,i]**2) #Poids du produit
	tau = np.hstack((-1,tau,1)) # -1,1(0,N+1)Ajouter
	D = np.zeros([n,n+1]) #Matrice différentielle de Gauss(La matrice différentielle gaussienne donne le différentiel de ce point.)
	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): #Équation du mouvement
	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): #Fonction d'évaluation
	return p[-1]

def equality(p, D, w, t0, div, n): #Condition de contrainte d'équation
	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): #Contrainte d'inégalité
	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 #Heure de début
	n = 20 #Maille discriminée (Legendre)-Nombre de nœuds de Gauss)
	num_p = 5 #Total des types de variables d'état et de contrôle
	iterator = 10 #Nombre d'itérations d'optimisation par défaut:1
	tau, w, D = make_node_derivative_matrix(n)
	div = [i*(n+2) for i in range(1,num_p+1)] #Index à la transition des variables
	num_variables = div[-1] + 1 #Nombre total de variables d'état, de variables de contrôle et de variables statiques

	#Vecteurs de limites supérieure et inférieure
	# 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 #Valeurs initiales des variables à optimiser
	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))

	#Décomposer les variables résultantes
	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

	#RÉSULTAT 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

Les références

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） ↑ Il est recommandé de lire le papier du matériel original de la méthode.

Koichiro Kato «Approche optimale de contrôle technique de la non-linéarité» 1988 Puisqu'il utilise la méthode de variation, il n'est pas directement lié, mais il semble être célèbre comme un livre traitant des problèmes de contrôle optimal.

Takeshi Tsuchiya et Shinji Suzuki. «Étude sur la résolution optimale des problèmes de contrôle à l'aide de la programmation mathématique (partie 2) Proposition de la méthode du cyan en diagonale par blocs.» 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 et Takeshi Tsuchiya. "Étude sur l'optimisation de la trajectoire ascendante d'un avion spatial vers l'ISS." JOURNAL DE LA SOCIÉTÉ JAPONAISE POUR LES SCIENCES AÉRONAUTIQUES ET SPATIALES 49.570 (2001): 208-213. https://www.jstage.jst.go.jp/article/jjsass/49/570/49_570_208/_article/-char/ja/

↑ Il existe plusieurs articles de Suzuki et Tsuchiya Labs de l'Université de Tokyo, mais ils sont très utiles.

[PYTHON] Bases de la génération d'orbite (résolution d'un problème de contrôle optimal non linéaire en utilisant l'optimisation de Scipy)

Exemple

Les références