Как я должен дросселировать свою ракету, чтобы достичь максимальной высоты? [закрыто]

Это задача оптимального управления , поэтому я буду использовать правила оптимального управления. Во-первых, мы представляем уравнения пространства состояний. Также мы берем общую массу как состояние и количество сожженного топлива как входное управление. Итак, у нас есть:

$$\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}$$ с этими граничными условиями: $$\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}$$ Мы хотим найти функцию $\theta(t)$которые минимизируют следующую функцию стоимости: $$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}$$ В приведенных выше уравнениях конечное время $t_f$неизвестно и должно быть определено.

Для нахождения решения определим гамильтониан этой системы:

$$\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*}$$ где $a$– вектор дифференциальных уравнений состояний и $(p_1,p_2,p_3)$являются стоимостными переменными.

Используя вариационное исчисление, найдем необходимые условия оптимальности:

$$\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}$$ и мы получаем это граничное условие: $$\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*}$$

Использование уравнений$(5)$мы можем вывести дифференциальные уравнения состояний и костатов:

$$\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}$$ и из $0=\frac{\partial H}{\partial \theta}$в $(5)$мы получаем: $$\frac{p_2}{p_3}=\frac{x_3}{\eta} \tag{8}$$ Также из уравнения $(6)$эти окончательные условия будут найдены: $$\begin{cases}\tag{9} p_1(t_f)=0\\ p_2(t_f)=0 \end{cases}$$ и: $$\begin{equation}\tag{10} H(t_f)=0 \end{equation}$$

Теперь у нас есть семь неизвестных функций$\Bigl(x_1(t),x_2(t),x_3(t),p_1(t),p_2(t),p_3(t),\theta(t)\Bigl)$с шестью дифференциальными уравнениями в$(7)$плюс одно уравнение связи в$(8)$. Для решения этих уравнений нужны достаточные граничные условия. Использование четырех граничных условий в уравнении$(2)$и два конечных условия в уравнении$(9)$все необходимые граничные условия будут известны.

Мы не использовали уравнение$(10)$еще. Это уравнение будет использоваться для нахождения терминального времени$t_f$.

Наконец, для решения краевой задачи в$(7)$и из-за нелинейности мы должны использовать численные методы, такие как метод стрельбы .

Я подумал, что этот вопрос интересен, и я не хотел делать какую-либо работу сегодня днем, поэтому я сделал простую модель, чтобы узнать, что произойдет. Мой код Matlab находится в конце вопроса.

До сих пор я тестировал три случая и рассматривал возможность изменения начальной тяги и добавления линейного увеличения тяги для каждого случая. Тяга дается как доля от общего количества сжигаемого топлива в секунду.

Очевидно, что цифры, которые я использовал, не являются невероятно реалистичными или разумными, но они дают первоначальное представление о том, что происходит. Цвет на графиках представляет высоту, полученную с помощью этой тяги.

До сих пор я рассматривал только постоянные коэффициенты лобового сопротивления. Сначала я попробовал k=0 (без перетаскивания)

Как и ожидалось, это дает наибольшую высоту, когда вы сжигаете топливо так быстро, как только можете. Это в основном просто тест, чтобы убедиться, что я не сделал ничего плохого.

Если мы выберем k=1, мы получим

Интересно, что здесь мы видим обратный эффект. Вам нужна самая низкая и самая длинная тяга, которая оторвет вас от земли, поскольку сопротивление воздуха является доминирующей формой сопротивления.

Наконец, промежуточный случай k=0,1

Здесь существует баланс между сопротивлением воздуха и гравитацией, и предпочтение отдается промежуточному значению тяги.

Интересно (по крайней мере, для меня), что было невыгодно сбрасывать газ после запуска. Я подозреваю, что это может произойти, если рассматривается переменная k.

Вот мой код матлаба. Должен работать и в октаве.

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