Wie soll ich meine Rakete drosseln, um die höchste Höhe zu erreichen? [geschlossen]

Dies ist ein Optimalsteuerungsproblem , daher werde ich die Regeln der Optimalsteuerung verwenden. Zuerst stellen wir die Zustandsraumgleichungen dar. Außerdem nehmen wir die Gesamtmasse als Zustand und die verbrannte Kraftstoffmenge als Eingangskontrolle. Also haben wir:

$$\begin{cases} \tag{1} \dot{x}_1=x_2 \\ \dot{x}_2=\frac{\eta \theta-k(x_1)x_2^2}{x_3}-g \\ \dot{x}_3=-\theta \end{cases}$$ mit diesen Randbedingungen: $$\begin{cases} \tag{2} x_1(0)=0\\ x_2(0)=0\\ x_3(0)=m_d+m_f(0)\\ x_3(t_f)=m_d \end{cases}$$ Wir wollen Funktion finden $\theta(t)$die die folgende Kostenfunktion minimieren: $$J(\theta)=h(x(t_f),t_f)+\int_0^{t_f}g(x,\theta,t)dt=\int_{0}^{t_f}(-x_2)dt\tag{3}$$ In den obigen Gleichungen die Endzeit $t_f$ist unbekannt und muss ermittelt werden.

Um die Lösung zu finden, definieren wir den Hamiltonoperator dieses Systems:

$$\begin{eqnarray*}\tag{4} H(t,\theta,x,p)&=&g+p^Ta\\ &=&-x_2+p_1x_2+p_2(\frac{\eta \theta-k(x_1)x_2^2}{x_3}-g)-p_3\theta \end{eqnarray*}$$ Wo $a$ist der Vektor der Zustandsdifferenzialgleichungen und $(p_1,p_2,p_3)$sind costate Variablen.

Unter Verwendung der Variationsrechnung werden die notwendigen Bedingungen für Optimalität gefunden:

$$\begin{cases}\tag{5} \dot{x}=\frac{\partial H}{\partial p}\\ \dot{p}=-\frac{\partial H}{\partial x}\\ 0=\frac{\partial H}{\partial \theta} \end{cases}$$ und wir erhalten diese Randbedingung: $$\begin{eqnarray*}\tag{6} &&\Biggl(\frac{\partial h}{\partial x}(x(t_f),t_f)-p(t_f) \Biggl)\delta x_f\\ &&+ \Biggl(H(x(t_f),p(t_f),\theta(t_f),t_f)+\frac{\partial h}{\partial t}(x(t_f),t_f) \Biggl)\delta t_f \\ && =0 \end{eqnarray*}$$

Verwenden von Gleichungen$(5)$wir können Zustände und Differenzialgleichungen ableiten:

$$\begin{cases} \tag{7} \dot{x}_1=x_2\\ \dot{x}_2=\frac{\eta \theta-k(x_1)x_2^2}{x_3}-g \\ \dot{x}_3=-\theta\\ \dot{p}_1=\frac{p_2x_2^2}{x_3}\frac{\partial k}{\partial x_1}\\ \dot{p}_2=1-p_1+2\frac{p_2k(x_1)x_2}{x_3}\\ \dot{p}_3=-\frac{p_2}{x_3^2}(\eta \theta-k(x_1)x_2^2) \end{cases}$$ und von $0=\frac{\partial H}{\partial \theta}$In $(5)$wir bekommen: $$\frac{p_2}{p_3}=\frac{x_3}{\eta} \tag{8}$$ Auch aus Gleichung $(6)$Diese Endbedingungen werden gefunden: $$\begin{cases}\tag{9} p_1(t_f)=0\\ p_2(t_f)=0 \end{cases}$$ Und: $$\begin{equation}\tag{10} H(t_f)=0 \end{equation}$$

Jetzt haben wir sieben unbekannte Funktionen$\Bigl(x_1(t),x_2(t),x_3(t),p_1(t),p_2(t),p_3(t),\theta(t)\Bigl)$mit sechs Differentialgleichungen in$(7)$plus eine Nebenbedingungsgleichung in$(8)$. Zur Lösung dieser Gleichung benötigen wir hinreichende Randbedingungen. Verwendung von vier Randbedingungen in Gleichung$(2)$und zwei Endbedingungen in Gleichung$(9)$alle benötigten Randbedingungen sind bekannt.

Wir haben keine Gleichung verwendet$(10)$noch. Diese Gleichung wird verwendet, um die Endzeit zu finden$t_f$.

Schließlich zur Lösung des Randwertproblems in$(7)$und wegen der Nichtlinearität müssen wir numerische Methoden wie die Schießmethode verwenden .

Ich fand diese Frage interessant und wollte heute Nachmittag keine richtige Arbeit machen, also habe ich ein einfaches Modell gemacht, um herauszufinden, was passieren würde. Mein Matlab-Code ist am Ende der Frage.

Bisher habe ich drei Fälle getestet und überlegt, den Anfangsschub zu ändern und für jeden Fall eine lineare Erhöhung des Schubs hinzuzufügen. Der Schub wird als Bruchteil des gesamten pro Sekunde verbrannten Treibstoffs angegeben.

Offensichtlich sind die Zahlen, die ich verwendet habe, nicht unglaublich realistisch oder vernünftig, aber sie geben eine erste Vorstellung davon, was vor sich geht. Die Farbe in den Diagrammen stellt die Höhe dar, die mit diesem Schub erreicht wird.

Bisher habe ich nur konstante Luftwiderstandsbeiwerte betrachtet. Zuerst habe ich versucht k = 0 (kein Widerstand)

Wie erwartet ergibt dies die größte Höhe, wenn Sie Kraftstoff so schnell wie möglich verbrennen. Dies ist hauptsächlich nur ein Test, um zu überprüfen, ob ich nichts falsch gemacht habe.

Wenn wir k=1 wählen, erhalten wir

Interessanterweise sehen wir hier den gegenteiligen Effekt. Sie wollen den niedrigsten und längsten Schub, der Sie vom Boden abhebt, da der Luftwiderstand die dominierende Form des Luftwiderstands ist.

Schließlich ein Zwischenfall von k = 0,1

Hier herrscht ein Gleichgewicht zwischen Luftwiderstand und Schwerkraft und ein mittlerer Schubwert wird bevorzugt.

Interessanterweise (zumindest für mich) war es nicht vorteilhaft, nach dem Start Gas zu geben. Ich vermute, dass dies auftreten kann, wenn die Variable k berücksichtigt wird.

Hier ist mein Matlab-Code. Sollte auch in Oktave laufen.

function[x]=rocket_test(thrustparam)

g=9.81;

md=1000;
mf=1000;  
k=0.1;
nu=1000;

dt=0.001;

v=0;
x=0;
t=0;

c=thrustparam(1)*mf;
p_1=thrustparam(2)*mf;

count=1;

while v>=0                  %want to find max height keep going until you start descending.

    throttle=throttle_set(t, c, p_1);

    if mf-throttle*dt<0    %if you would burn more fuel than left burn all the fuel  left;
        throttle=mf/dt;
    end

m=md+mf;
Fd=k*v^2+m*g;
Ft=nu*throttle;

F=Ft-Fd;

a=F/m;

x=x+v*dt+0.5*a*dt^2;
v=v+a*dt;
mf=mf-throttle*dt;
t=t+dt;

%store(count, :)=[t, x, v, a, throttle, Fd, Ft, F, m];
%count=count+1;

end
end

function [throttle]=throttle_set(t, c, p_1)

throttle=c+p_1*t;

if throttle<0
    throttle=0;
end

end