¿Cómo debo acelerar mi cohete para alcanzar la mayor altitud? [cerrado]

Este es un problema de control óptimo , así que usaré las reglas de control óptimo. Primero, representamos las ecuaciones del espacio de estado. También tomamos la masa total como estado y la cantidad de combustible quemado como control de entrada. Entonces tenemos:

$$\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}$$ con estas condiciones de contorno: $$\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}$$ Queremos encontrar la función. $\theta(t)$que minimizan la siguiente función de costo: $$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}$$ En las ecuaciones anteriores el tiempo terminal $t_f$es desconocido y debe ser determinado.

Para encontrar la solución, definimos el hamiltoniano de este sistema:

$$\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*}$$ dónde $a$es el vector de ecuaciones diferenciales de estado y $(p_1,p_2,p_3)$son variables de costa.

Utilizando el cálculo de variaciones, se encontrarán las condiciones necesarias para la optimalidad:

$$\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}$$ y obtenemos esta condición de contorno: $$\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*}$$

Usando ecuaciones$(5)$podemos derivar ecuaciones diferenciales de estados y costas:

$$\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}$$ y de $0=\frac{\partial H}{\partial \theta}$en $(5)$obtenemos: $$\frac{p_2}{p_3}=\frac{x_3}{\eta} \tag{8}$$ También de la ecuación $(6)$estas condiciones finales serán: $$\begin{cases}\tag{9} p_1(t_f)=0\\ p_2(t_f)=0 \end{cases}$$ y: $$\begin{equation}\tag{10} H(t_f)=0 \end{equation}$$

Ahora, tenemos siete funciones desconocidas.$\Bigl(x_1(t),x_2(t),x_3(t),p_1(t),p_2(t),p_3(t),\theta(t)\Bigl)$con seis ecuaciones diferenciales en$(7)$más una ecuación de restricción en$(8)$. Para resolver estas ecuaciones necesitamos suficientes condiciones de contorno. Usando cuatro condiciones de contorno en la ecuación$(2)$y dos condiciones finales en la ecuación$(9)$se conocerán todas las condiciones de contorno necesarias.

No usamos la ecuación$(10)$todavía. Esta ecuación se utilizará para encontrar el tiempo terminal$t_f$.

Finalmente, para resolver el problema de valores en la frontera en$(7)$y debido a la no linealidad, debemos usar métodos numéricos como el método de disparo .

Pensé que esta pregunta era interesante y no quería hacer ningún trabajo adecuado esta tarde, así que hice un modelo simple para averiguar qué sucedería. Mi código de matlab está al final de la pregunta.

Hasta ahora, probé tres casos y consideré cambiar el empuje inicial y agregar un aumento lineal en el empuje para cada caso. El empuje se da como una fracción del combustible total quemado por segundo.

Obviamente, los números que he usado no son increíblemente realistas ni sensatos, pero dan una idea inicial de lo que está pasando. El color de los gráficos representa la altura obtenida con ese empuje.

Hasta ahora solo he considerado coeficientes de arrastre constantes. Primero probé k = 0 (sin arrastre)

Como era de esperar, esto da la mayor altura cuando quema combustible lo más rápido que puede. Esto es principalmente solo una prueba para comprobar que no he hecho nada demasiado malo.

Si elegimos k=1 obtenemos

Curiosamente aquí vemos el efecto contrario. Desea el empuje más bajo y más largo que lo hará despegar, ya que la resistencia del aire es la forma dominante de arrastre.

Finalmente un caso intermedio de k=0.1

Aquí hay un equilibrio entre la resistencia del aire y la gravedad y se favorece un valor intermedio de empuje.

Curiosamente (al menos para mí), no fue ventajoso acelerar después del lanzamiento. Sospecho que esto puede ocurrir si se considera la variable k.

Aquí está mi código de matlab. Debería funcionar en octava también.

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