Over-engineering a Mun landing in Kerbal Space Program

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.

Methodology

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.

Apply 90% throttle at 10 degrees pitch, until the apoapsis reaches 10 km above the surface.
Coast to the apoapsis.
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.

https://gitlab.com/joost.meulenbeld/ksp_calculator

https://gitlab.com/joost.meulenbeld/ksp_interface