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.

  • $\begingroup$ The solution in the book gives discrete values for $\alpha (t)$ $\endgroup$ Mar 29, 2016 at 20:31
  • $\begingroup$ Sorry I forgot to mention that $\alpha(t)$ is either max (3000), or min (0), and does not attain any intermediate values. I think the max-min thrust profile was proved by Topcu et al. a few years ago with regards to powered descents. Thus it is a case of when to switch from min to max and vice-versa, and how "big" $tf$ should be. $\endgroup$ Mar 29, 2016 at 21:00
  • $\begingroup$ Isn't this a calculus of variations problem? $\endgroup$
    – user1722
    Mar 30, 2016 at 2:23
  • 1
    $\begingroup$ Btw, surface gravity on the moon is 1.62, not 9.81. Another thought: I think the switching of the engine has to be constrained by the minimal time between states. You can't have the engine continuously switching on and off. $\endgroup$
    – shrx
    Apr 11, 2016 at 8:44
  • $\begingroup$ Is it the same problem that is described here? $\endgroup$
    – C. E.
    Apr 16, 2016 at 18:20

2 Answers 2


Oh boy, what a question! This is very similar to some stuff I played a few weeks ago (Kerbal, what a game!). What follows solves (I think) the question you are asking.

An approach that seemed to to help was to split the problem into two: before, and after the engine burn. I do this with the knowledge that the most efficient landing will comprise just a single engine burn, ending at ground level, a so-called suicide burn, because of gravity losses.

I use the parameters:

g = -9.81;
v0 = 0;
h0 = 100000;
k = 0.001;
m0 = 2000;
T = 15000;

Before the engine burn, the equations of motion are simple free fall. Therefore if we wait τ seconds, we will be at height h0 + v0 τ + 0.5 g τ^2 and velocity v0 + g τ.

Solving for the equations of motion during the burn phase, we use the wonderful ParametricNDSolveValue, I've just discovered this function and it is magic:

pfun = ParametricNDSolveValue[{
   z''[t] == g + T/m[t] ,
   z'[0] == v0 + g τ,
   z[0] == h0 + v0 τ + 0.5 g τ^2,
   m'[t] == -k T,
   m[0] == m0
   }, z, {t, 0, 120}, {τ}]

enter image description here

As this only account for the motion after the engine is switched on, we prepend the freefall motion:

fullZ[τ_, t_] := If[t > τ, pfun[τ][t - τ], h0 + v0 t + 0.5 g t^2]
fullV[τ_, t_] := If[t > τ, pfun[τ]'[t - τ], v0 + g t]

Thus we can plot the trajectories for different engine switch on times:

Plot[Evaluate[Table[fullZ[τ, t], {τ, 70, 80, 2}]], {t, 0, 200}, PlotRange -> All]

enter image description here

We can see one of them gets very close to zero velocity at zero height..! We can find exactly the parameter τ which gives us the perfect landing, and the time of the landing with FindRoot (also amazing):

root = FindRoot[{fullZ[τ, tl], fullV[τ, tl]}, {{τ, 75}, {tl, 190}}]

{τ -> 75.7512, tl -> 187.999}

And there we have it! The most efficient landing is the one where the engines come on full blast after 75.75s and don't stop until the deck.

 Plot[Evaluate@Table[Tooltip[fullZ[τ, t], τ], {τ, 70, 80, 2}], {t, 0, 200}, PlotStyle -> Lighter[Gray, 0.5]],
 Plot[Evaluate[fullZ[τ, t] /. root], {t, 0, Evaluate[tl /. root]}, PlotStyle -> Orange]

enter image description here

If it looks like the velocity or height are not exactly zero at the landing it's just due to the plotting:

 Table[Evaluate[{t, fullZ[τ, t], fullV[τ, t]} /. root], {t, 183, 193., 0.5}],
 TableHeadings -> { None, {"t", "h(t)", "v(t)"}}

enter image description here

  • $\begingroup$ You almost have the same answer as I did, but you are not considering the mass: "The goal is to land softly while simultaneously using as little fuel as possible", $\endgroup$
    – MathX
    Apr 16, 2016 at 17:36
  • $\begingroup$ The problem wants the optimum solution WHILE having a soft landing. So the soft landing is more of a constraint, but the least fuel consumption is the objective function. Add a column of mass to your table. Hopefully it will not become less than 100 at 188? (The rockets mass after all the fuel has been burnt). $\endgroup$
    – MathX
    Apr 16, 2016 at 17:40
  • $\begingroup$ The mass is considered, you can see it in the equations of motion. You are right that I do not explicitly maximise the landing mass, this is because I make use of the fact that it is known that turning the engines on once and at the latest possible time implicitly minimises fuel spent and therefore maximises landing mass. Which form the looks of it is what you do as well? My ParametricNDSolve does what your NDSolve does, but with the switchon time still as a free parameter, which I then find the optimal values for with FindRoot, rather than your brute force approach. $\endgroup$ Apr 16, 2016 at 18:10
  • 1
    $\begingroup$ Yes, this is very close to my answer (316.675). $\endgroup$
    – MathX
    Apr 16, 2016 at 18:10
  • $\begingroup$ But for a general case, you haven't restricted the mass consumption in your solution, which could become negative.I could have gone further with root finding but it didn't seem to solve the maximization problem so it is like solving an under-determined problem, $\endgroup$
    – MathX
    Apr 16, 2016 at 18:13

This is more of a long comment, hopefully you can find some of these ideas helpful because I don't know how to completely implement this. Therefore, I'd really appreciate comments from more experienced users about whether they think this is feasible or not.

If we define a list of "switching times" and then define a function that takes this list as an argument and solves the system of ODE's, we might be able to minimize the objective function by finding the optimum list of switching times.

I put a couple hours on that part but I gave up because I don't know how to define a list of variable length and then optimize it. However, I am sharing what I found so far:

First, in my code, I got better results by defining v separately. Second, I think it's better if you use your thrust as a discrete variable. moreover, it was easier for me to do it as a binary variable and put the 15000 in the equations directly. You can also use WhenEvent to stop the integration. This will give you a few advantages (i.e. no need to define constraints for h[t] and m[t] as long as their initial value is correct and, no need to use stiffness switching method). I didn't go that far but using the domain of the solution's interpolating functions gives you your free variable time tf (I have done this for another problem using the "Domain" attribute and it worked fine). As a demonstration, I bring the code of how the "switching" is done every 20 seconds:

sol= NDSolve[{v'[t] == -g + 15000 power[t] /m[t], h'[t] == v[t], 
      m'[t] == -k 15000 power[t] , h[0] == 100000, v[0] == 0, 
      m[0] == 2000, power[0] == 0, 
      WhenEvent[h[t] == 0, "StopIntegration"], 
      WhenEvent[m[t] == 100, "StopIntegration"], 
      WhenEvent[Mod[t, 20] == 0, power[t] -> Abs[power[t] - 1]]}, {h[t], 
      v[t], m[t], power[t]}, {t, 0, 100000}, 
     DiscreteVariables -> {power[t] \[Element] {0, 1}}]

Which by no means gives a safe landing velocity but imo was an improvement:


enter image description here

Power (or switching time):

enter image description here

Finally, another method might be to use ParametricNDSolve and solve an equation that gives the landing velocity of zero?

Again, please remember that this post is more of a comment, not an answer.

EDIT 1: In this solution, when m becomes equal to 100, all the fuel has been consumed and the power should be zero from there on. So what I did was wrong. However, since m=100 is the maximum fuel consumption, it will not be the optimum solution anyways.

EDIT 2: On second thought, when I read your question again, it seems like you are talking about only one switching time. If that's, it is much easier:

switchon[ts_] := 
 Module[{ts1 = ts}, 
  sol1 = NDSolve[{v'[t] == -g, h'[t] == v[t], h[0] == 100000, 
     v[0] == 0}, {h[t], v[t]}, {t, 0, ts1}]; 
  sol2 = NDSolve[{v'[t] == -g + 15000 power[t]/m[t], h'[t] == v[t], 
     m'[t] == -k 15000 power[t], 
     h[0] == (h[t] /. sol1 /. t -> ts1)[[1]], 
     v[0] == (v[t] /. sol1 /. t -> ts1)[[1]], m[0] == 2000, 
     power[0] == 1, WhenEvent[h[t] == 0, tf = t; "StopIntegration"], 
     WhenEvent[m[t] == 100, power[t] -> 0]}, {h[t], v[t], m[t], 
     power[t]}, {t, 0, 100000}, 
    DiscreteVariables -> {power[t] \[Element] {0, 1}}]; 
  Return[{ts1 + 
     tf, (m[t] /. sol2 /. t -> tf)[[1]], (v[t] /. sol2 /. 
       t -> tf)[[1]]}]]

This way, switchon gives a list of the form {${tf, m(tf),v(tf)}$} which is the landing time and the mass and velocity at this time, given as a function of our switching time. With a single switching time, for current initial values the landing velocity as a function of switching time is:

enter image description here

Based on this, (using a brute force approach) I got a minimum velocity magnitude of 0.987722 at switching time of 75.7512, with the mass of 316.675. Please keep in mind that with this approach I only tried to find the switching time that gives us a soft landing. In other words, m is relaxed, and we are only solving for v==0. It doesn't look like a single-objective problem if you want to find the maximum of m and zero landing velocity at the same time. Maybe the formulation should be changed? Maybe minimizing the landing v/m ?


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.