Soft Landing on the Moon with Mathematica

Marian Mureşan

We consider the problem of landing a spacecraft on the Moon, assuming that aerodynamic and gravitational forces of bodies other than the Moon are negligible, and lateral motion can be ignored. Accordingly, the descent trajectory is vertical, and the thrust vector is tangent to the trajectory.

Because the spacecraft is near the Moon, we assume that the lunar acceleration of gravity has the constant value g=1.63, that the relative velocity of the exhaust gases with respect to the spacecraft is constant, and that the mass rate m'(t) is constrained by -μ⩽m'(t)⩽0, where μ is constant and gives the maximum rate of change of the mass due to burning the fuel.

Mathematical Approach to a Soft Landing

The problem is to make a soft landing on the surface of the Moon with the minimum amount of fuel.

Here is a sketch of the system immediately preceding the landing.

Graphics[{
{Thick,Line[{{-8,0},{8,0}}]},
Line[{{-2,5},{2,5},{2,7},{-2,7},{-2,5}}],Line[{{-1,7},{1,7},{1,8},{-1,8},{-1,7}}],
Line[{{-2,5},{-3,4}}],
Line[{{2,5},{3,4}}],
{Red,Line[{{0,5},{#,4}}]&/@Range[-.4,.4,.2]},
Style[Text[Row[{"total mass = ",Style["m",Italic],"(",Style["t",Italic],")"}],{-5,6}],12],
Arrow[{{0,6},{0,11}}],
Style[Text[Row[{"thrust = ","-",Style["k",Italic]," ",Style["m",Italic],"'(",Style["t",Italic],")"}],{3,9.5}],12],
{Dashed,Arrow[{{5,0},{5,6}}]},
Style[Text[Row[{"height = ",Style["h",Italic],"(",Style["t",Italic],")"}],{7,3}],12],
Arrow[{{0,6},{0,2}}],
Style[Text[Row[{Style["m",Italic]," ",Style["g",Italic]}],{1,3}],12],
PointSize[.015],Point[{0,6}]
}]

Following [, pp. 247–248] and [], we introduce the following notation and assumptions:

t is time

m(t) is the mass of the spacecraft, which varies as fuel is burned

m'(t) is the rate of change of mass, constrained by -μ≤m'(t)⩽0

g=1.63, the gravitational constant near the Moon

k is a constant, the relative velocity of the exhaust gases with respect to the spacecraft

T(t)=-k m'(t), the thrust

h(t) is the the height, with h(t)≥0

v(t)=h'(t), the velocity

u(t)=m'(t), the control function

Recalling our assumptions, aerodynamic forces and gravitational forces of bodies other than the Moon are negligible and lateral motion is ignored. Thus the descent trajectory is vertical and the thrust vector is perpendicular to the ground.

We also suppose that (Subscript[m, 0])=m(0)=M+F, where M is the mass of the spacecraft without fuel and F is the initial mass of fuel; m(t)>M, since as we expect that the spacecraft will return to Earth, it needs some fuel for takeoff.

Equations of Motion

By Newton’s second law ([, p. 128] and []),

m(t)h''(t)=-g m(t)+T(t)=-g m(t)-k m'(t),

which can be written as a system of equations

h'(t)=v(t),

v'(t)=-g-k ((u(t))/(m(t))),

m'(t)=u(t),

where k is a constant. The third equation states that the loss of mass per second (the fuel burned by the jet per second) is proportional to the thrust of the jet.

The Optimal Control Problem

Our goal is to minimize the fuel consumption, so the cost functional is

Λ(h,v,m,u)=-(∫_0^b)m'(t)dt=m(0)-m(b),

where b is the first time for which

({{h(b)=v(b)=0.}})

Thus the horizon is [0,b], where b remains to be determined.

In vector form, if x(t)=((Subscript[x, 1])(t),(Subscript[x, 2])(t),(Subscript[x, 3])(t))=(h(t),v(t),m(t)), then

x'(t)=f[x(t),u(t)],

and the problem can be written

(Subscript[x, 1])'=(Subscript[x, 2,])
(Subscript[x, 2])'=-g-k (Subscript[x, 3]^(-1)) u,
(Subscript[x, 3])'=u (dynamics of the motion),
(Subscript[x, 1])(0)=h(0),
(Subscript[x, 2])(0)=(Subscript[v, 0]),
(Subscript[x, 3])(0)=M+F (initial constraints),
(Subscript[x, 1])(b)=0,(Subscript[x, 2])(b)=0,(Subscript[x, 3])(b)≥M (final conditions),
u:[0,b]⟶[-μ,0] (u integrable,bounded control), 
t∈[0,b],b≥0 (finite horizon),

and finally,

Λ(x,u)=-(∫_0^b)((d(Subscript[x, 3])(t))/(dt))dt=(Subscript[x, 3])(0)-(Subscript[x, 3])(b) (the cost functional to minimize).

From () we have that

v'(t)=-g-k (d/(dt))ln(m(t)),

and by integration over the interval [0,t], we get that

h'(t)=-g t-k ln(((m(t))/(m(0))))+h'(0).

It follows that h'(b)=0 if and only if

k ln(((m(b))/(m(0))))=h'(0)-g b.

Solving for m(b),

m(b)=m(0)(e^((h'(0)-g b)/k)).

Now we substitute this into () to get

Λ(u)=m(0)(1-(e^((h'(0)-g b)/k))).

This result was published in [].

Theorem 1

Consider equation () with conditions () and (). Then for some given h(0), h'(0), g, and k, the amount of fuel required to stop the spacecraft, that is, to force h(b)=v(b)=0, is a strictly increasing monotonic function of the terminal time b.

Corollary 1

Minimizing the terminal time b for equation () with () and () is equivalent to minimizing the fuel consumption.

From here it follows that instead of () we can consider the following cost functional

Λ(x,u)=(∫_0^b)dt=b to be minimized,

and thus equation () with () becomes () with (). This is a Mayer optimal control problem (see Chapter 4 of []).

Necessary Conditions for the Mayer Problem

To avoid a lengthy discussion, we state a short version of theorem 4.2.i in []. Let the Mayer problem be expressed as

Λ(x,u)=g(a,x(a),b,x(b)) (cost functional),

((dx )/(dt))=f(t,x(t),u(t)), t∈[a,b] (a.e. (almost everywhere)), (differential constraint),
e[x]=g(a,x(a),b,x(b))∈B⊂(ℝ^(1+n+1+n)) (boundary conditions),
(t,x(t))∈T,t∈[a,b] (time-state constraint),
u(t)∈U(t),t∈[a,b] (control constraint).

A pair (x(t),u(t)), a≤t≤b, is said to be admissible or feasible provided x:[a,b]->(ℝ^n) is absolutely continuous [], u:[a,b]->(ℝ^m) is measurable, and x and u satisfy (). Let Ω be the class of admissible pairs (x,u). The goal is to find the minimum of the cost functional Λ(x,u) over Ω, that is, to find an element ((x^⋆),(u^⋆))∈Ω such that -∞<Λ((x^⋆),(u^⋆))≤Λ(x,u) for all (x,u)∈Ω. We introduce the variables p=((Subscript[p, 1]),…,(Subscript[p, n])), called multipliers, and an auxiliary function H(t,x,u,p), called the Hamiltonian, defined on T×U×(ℝ^n) by

H(t,x,u,p)=(∑_(i=1)^n)(Subscript[p, i])(Subscript[f, i])(t,x,u).

We define

M(t,x,p)=(inf_(u∈U(t)))H(t,x,u,p).

More assumptions are necessary:

There exists an element ((x^⋆),(u^⋆))∈Ω such that -∞<Λ[(x^⋆),(u^⋆)]≤Λ[x,u] for all (x,u)∈Ω.

T is closed in (ℝ^(1+n)).

The set S={(t,x,u)|(t,x)∈T, u∈U(t)} is closed in (ℝ^(1+m+n)).

f∈(C^1)(S,ℝ).

Notation:

(Subscript[f, it])=((∂(Subscript[f, i]))/(∂t)), (Subscript[f, i,(Subscript[x, j])])=((∂(Subscript[f, i]))/(∂(Subscript[x, j]))), (Subscript[H, (Subscript[x, j])])=((∂H)/(∂(Subscript[x, j])))=(∑_(i=1)^n)(Subscript[p, i]) (Subscript[f, i,(Subscript[x, j])]), (Subscript[H, t])=((∂H)/(∂t))=(∑_(i=1)^n)(Subscript[p, i]) (Subscript[f, i,t]), (Subscript[H, (Subscript[p, j])])=((∂H)/(∂(Subscript[p, j])))=(Subscript[f, j]).

The graph {(t,(x^⋆)(t))|a≤t≤b} of the optimal trajectory (x^⋆) belongs to the interior of T.

U does not depend on time and is a closed set.

The end point e((x^⋆))=g(a,(x^⋆)(a),b,(x^⋆)(b)) of the optimal trajectory (x^⋆) is a point of B, where B possesses a tangent variety B' (of some dimension k, 0≤k≤2n+2), whose vectors are denoted by

h={(Subscript[τ, 1]),(Subscript[ξ, 1]),(Subscript[τ, 2]),(Subscript[ξ, 2])},(Subscript[ξ, 1])={(Subscript[ξ, 1]^1),…,(Subscript[ξ, 1]^n)},(Subscript[ξ, 2])={(Subscript[ξ, 2]^1),…,(Subscript[ξ, 2]^n)},

or by

h={d(Subscript[t, 1]),d(Subscript[x, 1]),d(Subscript[t, 2]),d(Subscript[x, 2])},d(Subscript[x, 1])={d(Subscript[ξ, 1]^1),…,d(Subscript[ξ, 1]^n)},d(Subscript[x, 2])={d(Subscript[ξ, 2]^1),…,d(Subscript[ξ, 2]^n)}.

Theorem 2

Assume these eight hypotheses and let ((x^⋆)(u^⋆)) be an optimal pair of the Mayer problem () and (). Then the optimal pair ((x^⋆),(u^⋆)) necessarily has the following properties:

((d(Subscript[p, i]))/(dt))=(Subscript[H, (Subscript[x, i])])[t,(x^⋆)(t),(SuperStar[x])'(t)],i=1,…,n,t∈[a,b] (a.e.).

((dM)/(dt ))=(d/(dt))M[t,x(t),p(t)]=(Subscript[H, t])[t,x(t),u(t),p(t)],t∈[a,b] (a.e.).

The Hamiltonian and the equations for the multipliers to () and () are

H=(Subscript[p, 1])(Subscript[x, 2])+(Subscript[p, 2])(-g-k(Subscript[x, 3]^(-1))u)+(Subscript[p, 3])u=((Subscript[p, 1])(Subscript[x, 2])-(Subscript[p, 2])g)+(-(Subscript[p, 2])k(Subscript[x, 3]^(-1))+(Subscript[p, 3]))u,

(Subscript[p, 1])'=0, (Subscript[p, 2])'=-(Subscript[p, 1]), (Subscript[p, 3])'=(Subscript[p, 2])k (Subscript[x, 3]^(-2)),

so that (Subscript[p, 1])=(Subscript[c, 1]), (Subscript[p, 2])=-(Subscript[c, 1])t+(Subscript[c, 2]), 0≤t≤b, where (Subscript[c, 1]) and (Subscript[c, 2]) are constants.

For -(Subscript[p, 2])k (Subscript[x, 3]^(-1))+(Subscript[p, 3])>0, the minimum of H is attained with u=0, and then

H=(Subscript[p, 1])(Subscript[x, 2])-(Subscript[p, 2])g.

This corresponds to free fall for the spacecraft.

For -(Subscript[p, 2])k(Subscript[x, 3]^(-1))+(Subscript[p, 3])<0, the minimum of H is attained with u=-μ, and then

H=(Subscript[p, 1])(Subscript[x, 2])+(Subscript[p, 2])(-g-k(Subscript[x, 3]^(-1))u)+(Subscript[p, 3])u=((Subscript[p, 1])(Subscript[x, 2])-(Subscript[p, 2])g)-(-(Subscript[p, 2])k(Subscript[x, 3]^(-1))+(Subscript[p, 3]))μ.

Thus we find that the control function u takes only extreme values 0 and μ.

If on an interval [(Subscript[τ, 1]),(Subscript[τ, 2])] we have that u=0, and hence

(Subscript[x, 1])'=(Subscript[x, 2]),(Subscript[x, 2])'=-g,(Subscript[x, 3])'=0,(Subscript[p, 3])'=0,

then for t∈[(Subscript[τ, 1]), (Subscript[τ, 2])], we have

(Subscript[x, 3])(t)=(Subscript[x, 3])((Subscript[τ, 1])),

(Subscript[x, 2])(t)=(Subscript[x, 2])((Subscript[τ, 1]))-g(t-(Subscript[τ, 1])),

(Subscript[x, 1])(t)=(Subscript[x, 1])((Subscript[τ, 1]))-(1/2)g(((t-(Subscript[τ, 1])))^2)+(Subscript[x, 2])((Subscript[τ, 1]))(t-(Subscript[τ, 1])).

In this case, ((Subscript[x, 1]),(Subscript[x, 2])) describes an arc of a parabola of equation

(Subscript[x, 1])=-0.5 (g^(-1))(Subscript[x, 2]^2)+c(Subscript[x, 2])+d, with  (Subscript[x, 2])<0.

If on an interval [(Subscript[τ, 1]),(Subscript[τ, 2])] we have u=-μ, and hence

(Subscript[x, 1])'=(Subscript[x, 2]),(Subscript[x, 2])'=-g+k μ (Subscript[x, 3]^(-1)), (Subscript[x, 3])'=-μ,

then for t∈[(Subscript[τ, 1]), (Subscript[τ, 2])], we find

(Subscript[x, 3])(t)=(Subscript[x, 3])((Subscript[τ, 1]))-μ(t-(Subscript[τ, 1])),

(Subscript[x, 2])(t)=(Subscript[x, 2])((Subscript[τ, 1]))-g(t-(Subscript[τ, 1]))-k ln[1-(Subscript[x, 3]^(-1))((Subscript[τ, 1]))μ(t-(Subscript[τ, 1]))],

(Subscript[x, 1])(t)=(Subscript[x, 1])((Subscript[τ, 1]))+(Subscript[x, 2])((Subscript[τ, 1]))(t-(Subscript[τ, 1]))+k(t-(Subscript[τ, 1]))-0.5((g(t-(Subscript[τ, 1])))^2)+k (μ^(-1))[(Subscript[x, 3])((Subscript[τ, 1]))+μ(t-(Subscript[τ, 1]))]ln[1-(Subscript[x, 3]^(-1))((Subscript[τ, 1]))μ(t-(Subscript[τ, 1]))].

Theorem 3

If all the assumptions in the section Mathematical Approach to a Soft Landing are satisfied, τ is the time of ignition of the engine, b is the time of landing, that is (Subscript[x, 1])(b)=(Subscript[x, 2])(b)=0, and if the system of equations

(Subscript[x, 2])((Subscript[τ, 1]))=g(b-(Subscript[τ, 1]))-k ln(1-(Subscript[x, 3]^(-1))((Subscript[τ, 1]))μ(b-(Subscript[τ, 1]))),

(Subscript[x, 1])((Subscript[τ, 1]))=-((Subscript[x, 2])((Subscript[τ, 1]))+k)(b-(Subscript[τ, 1]))+0.5((g(b-(Subscript[τ, 1])))^2)+k (μ^(-1))[(Subscript[x, 3])((Subscript[τ, 1]))+μ(b-(Subscript[τ, 1]))]ln(1-(Subscript[x, 3]^(-1))((Subscript[τ, 1]))μ(b-(Subscript[τ, 1]))),

Program for Soft Landing on the Moon

MoonLanding is a Mathematica program for a soft landing on the Moon. Here h0 is the initial height, v0 is the initial velocity, mass is the mass of the lander without fuel, fuel is the initial fuel, g is acceleration due to gravity, k is the relative velocity of the exhaust gases, and μ is the rate of change of the mass by burning.

The correctness of the results drastically depends on the initial values of the variables z and g that we use in solving the nonlinear system of equations in the program.

MoonLanding[h0_,v0_,mass_,fuel_,g_,k_,μ_]:=Module[
{m,z,b,sol,t1,t2},
m=mass+fuel;
sol=FindRoot[{v0-g z-g(b-z)-k Log[1-(μ/m)(b-z)]==0,
-0.5 g (z^2)+v0 z+h0+(k+v0-g z)(b-z)-0.5 g (((b-z))^2)+(m k/μ)(1-(μ/m)(b-z))Log[1-(μ/m)(b-z)]==0},{{z,130},{b,200}}];
t1=z/.sol;
t2=b/.sol;
x1[t_]:=-0.5 g (t^2)+v0 t+h0;
x2[t_]:=v0-g t;
x32[t_]:=m-μ(t-t1);
x22[t_]:=x2[t1]-g(t-t1)-k Log[1-(μ/m)(t-t1)];
x12[t_]:=x1[t1]+(k+x2[t1])(t-t1)-0.5 g (((t-t1))^2)+(m k/μ)(1-(μ/m)(t-t1))Log[1-(μ/m)(t-t1)];

Show[{
ParametricPlot[{v0-g t,x1[t]},{t,0,t1}],
ParametricPlot[{x22[t],x12[t]},{t,t1,t2}],
Graphics[{
{
Arrow[{{-440,0},{800,0}}],
Arrow[{{0,-4000},{0,h0+h0/7}}]
},
{ 
Text[Style["P",Italic,12], {v0+30,h0}],
Text[Style["A",Italic,12], {x2[t1]-20,x1[t1]}] ,Text[Style["O",Italic,12], {20,1500}],
Text[(Subscript[Style["h",Italic,12], 0]), {40,h0}] ,Text[(Subscript[Style["v",Italic,12], 0]), {v0,-2000}],
Text[Style["(Subscript[h, A])",Italic,12], {40,x1[t1]}] ,
Text[Style["(Subscript[v, A])",Italic,12], {x2[t1],-2000}],
Text[Style["h",Italic,12], {-35,h0+h0/10}] ,
Text[Style["v",Italic,12], {775,2500}],

Text[Style[Row[{Style["P" ,Italic]," = (",(Subscript[Style["v" ,Italic], 0]),", ",(Subscript[Style["h" ,Italic], 0]),") "," = (",v0,", ",h0,")" }],12],{450,h0-h0/10}],
Text[Style[Row[{(Subscript[Style["m",Italic], 0])," = mass + fuel = ",m}],12],{450,h0-2 h0/10}],
Text[Style[Row[{Style["k",Italic]," = ",k}],12],{450,h0-3 h0/10}],
Text[Style[Row[{"μ = ",μ}],12],{450,h0-4 h0/10}],
Text[Style[Row[{Style["A" ,Italic]," = (",Style["(Subscript[v, A])" ,Italic],", ",Style["(Subscript[h, A])" ,Italic],") = (",Chop@x2[t1],", ",Chop@x1[t1],")"}],12], {450,h0-5 h0/10}]
},
{
PointSize[.01],
Point[{{v0,h0},{x2[t1],x1[t1]},{0,0},{0,h0},{0,x1[t1]},{v0,0},{x2[t1],0}}]
}
}]
},PlotRange->{{Automatic,All},{0,All}},ImageMargins->{{10,50},{10,10}},Axes->False,AxesOrigin->{0,0},AspectRatio->1/1.3]
]

MoonLanding[50000,-150,2010.51,719.63,1.63,636,16.5]

This Manipulate lets you vary the parameters in real time.

Manipulate[
Quiet@MoonLanding[h0,v0,mass,fuel,g,k,μ],
{{h0,50000,"intial height (Subscript[h, 0])"},5000,60000,1,Appearance->"Labeled"},
{{v0,-150,"initial velocity (Subscript[v, 0])"},0,-200,-1,Appearance->"Labeled"},
{{mass,2010.51,"mass without fuel"},1000,3000,1,Appearance->"Labeled"},
{{fuel,719.63,"initial fuel"},200,1500,1,Appearance->"Labeled"},
{{g,1.63,"acceleration due to gravity g"},1,5,.01,Appearance->"Labeled"},
{{k,636,"relative velocity of the exhaust gases k"},500,1000,1,Appearance->"Labeled"},
{{μ,16.5,"rate of change of the mass μ"},12,20,.01,Appearance->"Labeled"},
SaveDefinitions->True
]

Acknowledgments

The author expresses his gratitude to Horia F. Pop from Babes-Bolyai University, Faculty of Mathematics and Computer Science in Cluj-Napoca, Romania, for valuable discussions.

References

D. E. Kirk, Optimal Control Theory, Englewood Cliffs, NJ: Prentice-Hall, Inc., 1970.

J. Meditch, “On the Problem of Optimal Thrust Programming for a Lunar Soft Landing,” IEEE Transactions on Automatic Control, 9(4), 1964 pp. 477–484. doi: 10.1109/TAC.1964.1105758.

G. Leitmann (ed.), Optimization Techniques: With Applications to Aerospace Systems, Mathematics in Science and Engineering, Vol. 5, New York: Academic Press, 1962.

L. Cesari, Optimization—Theory and Applications, Problems with Ordinary Differential Equations, Applications of Mathematics, Vol. 17, New York: Springer, 1983.

M. Mureşan, A Concrete Approach to Classical Analysis, New York: Springer, 2009.