Over-engineering a Mun landing in Kerbal Space Program

30 Jan

This is an over-engineered attempt to perform a very efficient Mun landing, which I worked on together with my friend Joost Meulenbeld. Using kRPC, we have written an autopilot that executes the landing automatically.



For a perfectly spherical body, the ideal trajectory would be to place the periapsis just centimeters above the surface. An almost horizontal suicide burn then kills of all horizontal velocity. However, Mun is not a perfectly spherical body and this approach would risk us to fly into a mountain. The varying mass also prevents us from making a simple calculation.

Our method is to first generate a reference approach profile. We built a simulator that implements the orbital mechanics. Instead of starting the simulator in orbit and simulating our way towards the surface, the simulation is started on the landing spot. The fuel mass is increasing in this ‘reversed’ simulator, but everything else is the same. The simulator starts with a constant pitch angle of 10 degrees and 90% throttle. The angle can be lowered further, depending on obstacles on the surface of Mun. 90% throttle is chosen to leave some margin for control, as we want to do a precision landing of course.

The ‘reversed’ landing simulator

The state vector for the simulator is:

(1)   \begin{equation*} \mathbf{x} = \begin{bmatrix}r \\ \dot{r} \\\theta \\\dot{\theta} \\m_{fuel}\end{bmatrix}\end{equation*}

Where r is the radius from the center of Mun and \dot{r} is vertical speed. \theta is the angle relative to the landing location at the landing time and \dot{\theta} is the angular velocity. m_{fuel} is the remaining fuel mass.

The differential state equations are:

(2)   \begin{equation*}  \dot{\mathbf{x}} = \begin{bmatrix} \dot{r} \\ \sqrt{\frac{GM}{r^2}} + r\dot{\theta}^2 + \frac{T_{max}}{m} u_t\sin{\phi}  \\ \dot{\theta} \\ -\frac{2\dot{r}\dot{\theta}}{r} + r\dot{\theta}^2 + \frac{1}{r}\frac{T_{max}} {m}u_t\cos{\phi} \\ u_{t} \cdot \dot{m}_{max} \end{bmatrix}\end{equation*}

\dot{m}_{max} is positive! Mass is increasing in this ‘reversed’ simulator. u_t is throttle, T_{max} is the thrust delivered at 100% throttle and \dot{m}_{max} is the mass flow at 100% throttle. \phi is the pitch angle relative to the horizontal direction. m is the vehicle mass, which is the sum of the vehicle fuel mass m_{fuel} and dry mass m_{dry}. G and M are the gravitational constant and mass of Mun.

The differential equations are converted to python and are solved by the solve_ivp function of scipy, with time steps of 0.1s.

Getting to orbit

Now that we have a ‘reversed’ simulator, the task is to find a set of control inputs that results in a circular orbit at 10 km above the surface.

  1. Apply 90% throttle at 10 degrees pitch, until the apoapsis reaches 10 km above the surface.
  2. Coast to the apoapsis.
  3. Perform a circularization burn at apoapsis.

For step 1, we need to know the apoapsis of the orbit in order to stop burning in time. From the simulator, we only have the state vector which does not contain orbit parameters. The following calculations allow us to calculate orbit parameters just from r, \dot{r} and \dot{\theta}:

(3)   \begin{equation*} V = \sqrt{\dot{r}^2 + \left(r\dot{\theta}\right)^2}\end{equation*}

(4)   \begin{equation*} \gamma = \arctan{\frac{r\dot{\theta}}{\dot{r}}}\end{equation*}

(5)   \begin{equation*} a = \left(\frac{2}{r} - \frac{V^2}{GM} \right)^{-1}\end{equation*}

(6)   \begin{equation*} e = \sqrt{1 - \frac{r^2V^2 \sin{\gamma}^2}{GM \cdot a}}\end{equation*}

(7)   \begin{equation*} r_a = a \left( 1 + e\right)\end{equation*}

(8)   \begin{equation*} r_p = a \left( 1 - e\right)\end{equation*}

V is the length of the velocity vector, \gamma is the flight path angle, a is the semi-major axis of the ellipse and e is the eccentricity of the ellipse. r_a and r_p are the radius at apoapsis and periapsis respectively. \gamma is positive when approaching the apoapsis and negative when approaching the periapsis. Using these parameters and the current radius we can determine the position on the ellipse.

(9)   \begin{equation*} \cos{E} = \frac{1}{e} - \frac{r}{ae}\end{equation*}

(10)   \begin{equation*} E = \text{sgn}\left(\gamma\right) \cdot \arccos\left( \frac{1}{e} - \frac{r}{ae}\right)\end{equation*}

(11)   \begin{equation*} M = E - e\sin{E}\end{equation*}

(12)   \begin{equation*} \nu = \arccos{\frac{\cos{E} - e}{1 - e\cos{E}}}\end{equation*}

(13)   \begin{equation*} n = \sqrt{\frac{GM}{a^3}}\end{equation*}

E is the eccentric anomaly, M is the mean anomaly and \nu is the true anomaly. n is the mean motion in [rad/s]. Using these equations, we can calculate the time to apoapsis and velocity at apoapsis V_a:

(14)   \begin{equation*} \left( t_a - t \right) = \frac{\pi - M}{n}\end{equation*}

(15)   \begin{equation*} V_a = \sqrt{GM\left( \frac{2}{r_a} - \frac{1}{a}\right)}\end{equation*}

The required orbital velocity for a circular orbit at radius r is given by:

(16)   \begin{equation*} V_c = \sqrt{\frac{GM}{r}}\end{equation*}

The duration of the cirularization burn is estimated using a constant mass:

(17)   \begin{equation*} \Delta t_{burn} = \left(V_c - V_a\right) \frac{T_{max}}{m}ut\end{equation*}

For the circularization burn, a throttle setting u_t of 50% was used. The burn starts \Delta t_{burn}/2 seconds before reaching apoapsis.

Iterate to find correct fuel mass

The reversed simulator needs to start with an unknown quantity of fuel. The initial run can start with 0 kg. This will result in a certain amount of fuel in circular orbit. The amount of landed fuel is iterated a couple of times until the the orbit fuel mass matches with the fuel on board the ship.

The flight profile is now complete. Altitude and throttle are plotted here against time:

Compensate for rotation of Mun

The angle \theta where we need to start the de-orbit burn is relative to the landing location at t = 0. However, the landing location is moving with approximately 9m/s because Mun rotates. That’s more than 2km error if we don’t compensate. If the de-orbit burn starts at t_c. The longitude at which we need to start the first burn is:

(18)   \begin{equation*} \lambda_{start} = \lambda_{target} - \theta_{start} + \omega t_c\end{equation*}

where \omega is the rotational velocity of Mun in [rad/s].

Feedback control in the landing burn

The de-orbit burn is purely feed-forward. This needs to be timed precisely. A tiny deviation in the de-orbit burn can result in a large error at the start of the landing burn.

The reference trajectory describes the desired vehicle state for every moment in time. However, there will be small deviations from this profile that need to be corrected for. Since Mun itself is rotating, we must aim for the landing to happen at t=0. Every second difference increases the error by 9 meters. The variables we try to control are longitude and radius (altitude). A relatively simple controller is implemented that decouples these. To control longitude, we use thrust. For altitude control, we use the pitch angle.

Full video of the landing burn

Longitude control

The longitude error is defined by

(19)   \begin{equation*} \lambda_e = \lambda_{ref}\left(t\right) - \lambda\end{equation*}

where \lambda_{ref}\left(t\right) is the longitude in the reference trajectory at a given time and \lambda is the longitude of the ship. This is converted to a horizontal distance error, which is more practical for control.

(20)   \begin{equation*} x_e = \lambda_e \cdot r\end{equation*}

The throttle setpoint is defined by the following PD controller with feed-forward from the reference trajectory:

(21)   \begin{equation*} u_t = u_{t,ref}\left(t\right) - P_x \cdot x_e - D_x \cdot \frac{d}{dt}x_e\end{equation*}

The gains P_u and D_u are found kerbal-style by tuning. \u_t is bound between 0.8 and 1.0 (the reference trajectory has a constant value of 0.9).

Altitude control

The altitude error is:

(22)   \begin{equation*} h_e = r_{ref}\left(t\right) - r\end{equation*}

The pitch angle setpoint is controlled with the following PD-controller with feed-forward from the reference trajectory:

(23)   \begin{equation*} u_{\phi} = u_{\phi,ref}\left(t\right) + P_h \cdot h_e + D_h \cdot \frac{d}{dt} h_e\end{equation*}

These gains are also tuned kerbal-style, and the pitch setpoint is bound between 0 degrees and 20 degrees. The reference trajectory is 10 degrees.

Horizontal ‘suicide burn’

For the last 1000 meters, a different horizontal controller is used for improved landing accuracy. It is somewhat similar to a suicide burn, except that it requires slightly less than 100% throttle. When we are this close to the target, it is safe to ignore orbit dynamics and pretend that Mun is flat.

If the following constant horizontal acceleration setpoint \ddot{x}_{sp} would be applied, the horizontal velocity reaches zero exactly at the target:

(24)   \begin{equation*} \ddot{x}_{sp} = \frac{1}{2}\dot{x}^2 \cdot \frac{1}{x}\end{equation*}

where \dot{x} is the horizontal velocity and x is the horizontal distance to the target. The required horizontal acceleration can be translated to a horizontal thrust setpoint:

(25)   \begin{equation*} u_t = \ddot{x}_{sp} \cdot \frac{m}{T_{max}}\end{equation*}

Technically we would need a little more thrust to compensate for the pitch of the vehicle. Also, the target height is increased by 5 meters so that we can quickly rotate the vehicle upright for a smooth vertical touchdown.

What’s next

This post just shows the part where we get from orbit to the surface of the Mun, but we have already automated the entire flight starting at launch from Kerbin. The next step will be to build a Mun base that mines for fuel. Instead of landing on the surface, we can build vehicles that automatically dock to the Mun base and don’t require landing legs. Then, we can transfer massive amounts of fuel from the Mun base to an orbital station, where we can refuel the biggest rockets and spaceships to send to the edge of the Kerbol System and beyond!

Source code

For those who are interested: the code is open-source, but not documented.



Leave a Reply

Your email address will not be published. Required fields are marked *

Time limit is exhausted. Please reload CAPTCHA.

This site uses Akismet to reduce spam. Learn how your comment data is processed.