I was hoping to tackle optimal control using Mathematica in order to learn how I can use Mathematica's built in numerical integration and optimization functions together in order to solve an optimal control problem numerically (that is, without using a Hamiltonian approach). As such, I thought it would be best if I pose a simple scenario and possibly have someone explain to me how I can use Mathematica to solve the optimal control problem.
The following was taken from page 55 (example 4.4.4: Moon Lander) of this text.
I have the following set of differential equations: $$ \begin{cases} \dot{v}(t)=-g+ \frac{\alpha(t)}{m(t)} \\ \dot{h}(t)=v(t)\\ \dot{m}(t)=-k\times\alpha(t)\\ \end{cases} $$ where $h(t)$ is the height, $v(t)$ is the velocity, $m(t)$ is the spacecraft mass (which changes as fuel is burned) and $\alpha(t)$ is the thrust at time $t$. The following constraints are applied: $h(t)\geq0$, $m(t)\geq100$ and $0\leq\alpha(t)\leq15000$ for all $t$. Where the thrust is either set to $0$ (min, engine off) or $15000$ (max, engine on), but attains no values in between. And we have the following initial conditions (I just chose these arbitrarily): $h(0)=100000$, $v(0)=0$, $m(0)=2000$.
The goal is to land softly while simultaneously using as little fuel as possible, where the final time $tf$ is a free variable. That is, we want to do as little "thrusting" as possible and thus need to minimize the cost function (or equivalently maximize our final mass upon landing)
$$J=\min \int_{0}^{tf} \alpha(t)dt$$
And since the thrust $\alpha(t)$ can only be set to $0$ or $15000$, we essentially need to find the "switching time" that will turn the spacecraft's rocket engine on to full thrust such that we'll make a perfect soft landing at time $tf$.
I've used NDSolve[] before, so know how to solve systems of ODEs. Furthermore, I've used NMinimize[] as well in order to find an optimal parameter that does not change with time. However, now that the control $\alpha(t)$ changes with time, I simply have no clue how I can combine NDSolve[] and NMinimize[]/NMaximize[] in order to find the optimal solution for minimizing the cost function $J$ and meeting the constraints. Is such a thing possible and am I even approaching the above problem in the correct manner? As you can see, when it comes to actually writing code to solve the problem I'm unfortunately at a loss!
EDIT: This is all the code I have for this particular problem so far, which is very wrong indeed (I guess it should be treated at pseudo-code, and weak pseudo-code at that. Note that $\alpha(t)$ is $a[t]$ in the code)
g = 9.81;
k = 0.001;
ODESystem = NDSolve[{
h''[t] == -g + a[t]/m[t],
m'[t] == -k a[t],
h[0] == 100000, h'[0] == 0, m[0] == 2000}, {h[t], m[t]}, {t, 0,
tf} , MaxSteps -> 1000000, Method -> "StiffnessSwitching"];
Optimization =
NMaximize[{m[tf], h[t] >= 0, m[t] >= 100, h[tf] == 0, h'[tf] == 0,
0 <= a[t] <= 15000}, {tf}]
As can be seen, it's full of errors because I don't know how to incorporate a changing thrust a[t] and a free final time tf into NDSolve[]. This would presumably be found using NDSolve[] and NMaximize in parallel. Any help would be brilliant.
