Need help with math for suicide burn

[Ypsilon]

I'm trying to write a python script (using the kRPC mod) to do a suicide burn. However, so far I always reach a velocity of 0 a few hundred meters above the ground and then proceed to plummet to my death. My assumption is that the problem is a mathematical one, as the code is fairly simple. So I'd be very grateful if someone could check my math. Alternatively, if anyone has a different solution for what I'm trying to accomplish, I'd be happy to hear that to.

For simplicity's sake I assume that my spacecraft is always pointing straight up and falling straight down. I also assume that there is no atmosphere and thus no drag. So here is what I have come up with:

(Handwritten version of math)
t: Time since start of the burn in seconds (s)

Fe: Thrust of the spacecraft's engines in Newton (N)

M0: Wet Mass of the spacecraft at the beginning of the burn in Kilogram (Kg)

g: Surface gravity of the body I'm landing on in meters per second squared (m/s^2)

W: Weight of the spacecraft in Newton (N)

W = M0 * g

F: Net thrust of the engine when decelerating the spacecraft in Newton (N)

F = Fe - W

K: Fuel consumption of the engine in Kilograms per second (Kg/s)

a(t): Acceleration (or deceleration, depending how you look at it) of the spacecraft after t seconds, taking into account the decreased mass due to fuel being burned. In meters per second squared (m/s^2)

a(t) = F / (M0 - K * t)

Tb: Duration of engine burn in seconds (s)

Dv(Tb): Speed change of spacecraft for a burn of Tb seconds in meters per second (m/s)

Dv(Tb) = Integral from 0 to TB [a(t)] dt = Integral from 0 to Tb [F / (M0 - K * t)] dt = (F/K) * ln(1 + Tb * (K/M0))

Up until this point everything is based on this reddit comment. However, I am not sure if the integration for Dv(Tb) is correct. If I do it myself I get a different result, but the dimensional analysis for my result doesn't work out, so I've been sticking to this one.

v0: velocity of spacecraft at the beginning of the suicide burn in meters per second (m/s)

To calculate Tb, let

Dv(Tb) = (F/K) * ln(1 + Tb * (K/M0)) = v0

and solve for Tb:

Tb = (e^((v0 * K) / F) - 1) * (M0 / K)

Da(Tb): Distance the spacecraft falls during burn of Tb seconds in meter (m)

Da(Tb) = Integral from 0 to Tb [-v0 + Dv(Tb)] dt = Integral from 0 to Tb [-v0 + (F/K) * ln(1 + t * (K / M0))] dt

Because I couldn't find an analytical solution for Da(Tb) I decided to solve it numerically every tick.

 

I now run the above calculations every tick (multiple time per second) and then check if the resulting Da(Tb) is close to (within 10m) of the current altitude of my spacecraft. If it is I know I have to start the suicide burn. However, as mentioned above, it always starts a few hundred meters to early.

It shouldn't be relevant, but this is the spacecraft I'm using to test my code.

In case the above is to confusing to follow, here are my handwritten notes. They should contain the same information as above but more neatly formatted. Beware my handwriting though Also, is there a way to insert latex formulas in my post? I think that would have helped a lot.

And this is a copy of my code so far, if anyone is interested:

Spoiler

import time
import krpc
import math

import scipy.integrate as integrate
import numpy as np

engineThrust = 216394 #kN
fuelConsumption = 16.048 #Kg/s

def shouldBurn(vessel):
    rf = vessel.orbit.body.reference_frame
    velocity = -vessel.flight(rf).vertical_speed
    altitude = vessel.flight(rf).surface_altitude
    mass = vessel.mass
    gravity = vessel.orbit.body.surface_gravity
    weight = mass * gravity
    thrust = engineThrust - weight

    burnTime = ((math.exp((velocity * fuelConsumption) / thrust) - 1) * (mass / fuelConsumption))
    burnLength = integrate.quad(lambda t: -velocity + (thrust / fuelConsumption) * np.log(1 + t * (fuelConsumption / mass)), 0, burnTime)

    print("Burn time: " + str(burnTime) + " Burn start altitude: " + str(burnLength[0]))

    return altitude - burnLength[0] < 10

con = krpc.connect(name='Suicide-burn')
vessel = con.space_center.active_vessel

while(not shouldBurn(vessel)):
    time. sleep (0.001)

aStart = vessel.flight(vessel.orbit.body.reference_frame).surface_altitude
tStart = time.perf_counter()
print('burning!')
vessel.control.throttle = 1

while(vessel.flight(vessel.orbit.body.reference_frame).vertical_speed < -0.5):
    time.sleep(0.01)

vessel.control.throttle = 0
tEnd = time.perf_counter()
aEnd = vessel.flight(vessel.orbit.body.reference_frame).surface_altitude

print("Burned for " + str(tEnd - tStart) + "s over a distance of " + str(aEnd - aStart) + "m")

 

So, again, if anyone is willing to check my math, I'd be very grateful.

 

Edited April 30, 2020 by Ypsilon

[steuben]

Feels like you might be having a problem with floating point math. The errors are just compounding until things go poof.

What altitude/speed are you testing from? have you tried lower speeds/altitude? Can you land the bird manually?

[Laie]

(SCNR)

• You are calculating with surface gravity, yet local gravity changes quite noticeably with altitude.
   
• While I'm here, let me add another bit of advice: you cannot do a perfect, pre-calculated suicide burn. Your script only loops so fast, and KSP itself only updates the physics only every so often. So you can only hit the button at a few distinct points in time: one of these is bound to be a little too early, and the next one is already too late.
  
  It stands to reason that you should err on the side of caution, but then you have to take measurements along the way and adjust accordingly.
   
• Oh, and finally, I notice that you're working with a fixed burnTime = 2.7171 -- that might matter, too.

[Ypsilon]

35 minutes ago, steuben said:

Feels like you might be having a problem with floating point math. The errors are just compounding until things go poof.

What altitude/speed are you testing from? have you tried lower speeds/altitude? Can you land the bird manually?

I don't think floating point imprecision is the problem. I'm only doing ~10 operations and in the past I had to do a lot more before that became a problem.
I've been testing either at Kerbin or the Mun. At Kerbin from ~15000m. This gives me a velocity of ~250m/s when the burn starts at ~900m. It slows me to 0 m/s which I reach at ~150m above ground. At the Mun I've done tests from ~35000m. This gives me a top speed of ~300m/s. The burn starts at ~835m and finishes at ~135m.

I can totally land manually. The suicide burn leaves me about ~150m above the ground at 0m/s, so the hardest part from then on is to not overshoot the actual landing burn, as the engines are quite powerful for such a small craft.

[Ypsilon]

5 minutes ago, Laie said:

(SCNR)

• You are calculating with surface gravity, yet local gravity changes quite noticeably with altitude.
   
• While I'm here, let me add another bit of advice: you cannot do a perfect, pre-calculated suicide burn. Your script only loops so fast, and KSP itself only updates the physics only every so often. So you can only hit the button at a few distinct points in time: one of these is bound to be a little too early, and the next one is already too late.
  
  It stands to reason that you should err on the side of caution, but then you have to take measurements along the way and adjust accordingly.
   
• Oh, and finally, I notice that you're working with a fixed burnTime = 2.7171 -- that might matter, too.

I've been wondering about the surface gravity, but my thinking was that it only needs to be reasonably accurate for the last calculation I do, as thoose are the only values I'm actually using to control the spacecraft. As that is usually below 1000m I was hoping to get away with using the surface gravity. Do you think that's to optimistic?

I've been anticipating that my calculation won't be 100% accurate, but I think there's more at play here than that, as I am consistently and repeatably breaking to early. If the tick rate of my script was too slow, I'd expect to have more inconsistent results. As you say, sometimes burning too early, sometimes too late. But that's not the case.
 

You're right of course about the 2.7171. I'd put that in there for debugging purposes and forgot to take it out before posting. But the problem happens as described regardless.

[Laie]

2 hours ago, Ypsilon said:

As that is usually below 1000m I was hoping to get away with using the surface gravity. Do you think that's to optimistic?

Can't say. Compare local and surface gravity, multiply by burn time, and see how far that takes you. With a low TWR and a long burn time, even half a percent will add up.

I honestly cannot check your math. I cannot even follow most of it, sorry: My last integration was some 30 years ago. But. To me it looks as if you're making it needlessly complicated. For example,

4 hours ago, Ypsilon said:

Dv(Tb): Speed change of spacecraft for a burn of Tb seconds in meters per second (m/s)

Dv(Tb) = Integral from 0 to TB [a(t)] dt = Integral from 0 to Tb [F / (M0 - K * t)] dt = (F/K) * ln(1 + Tb * (K/M0))

If you already know the burn time and fuel flow, you know mass before/after. Why not plug that into the rocket equation? That's a tried-and-true solution you can copy verbatim out of a textbook and be done with it.

FWIW (digs in old scripts), here's the core of my low-math solution:

Spoiler

def impact_calc(alt, v):
    # I don't know how to directly calculate time-to-impact and impact velocity.
    # so I'm tracking back to a virtual point why my free fall has started,
    # do a simple free fall from v=0,
    # then adjust for how long and far I've fallen since then.
    alt = max(0, alt)
    g = vessel.orbit.body.surface_gravity
    tt0 = v / g                       # time to zero velocity
    alt0 = alt + (0.5 * g * tt0**2)   # altitude we started from
    tti0 = math.sqrt(2 * alt0 / g)    # time-to-impact, free fall from starting alt
    tti = tti0 + tt0                  # time remaining
    vi = math.sqrt(2 * g * alt0)      # impact velocity
    dist = (0.5 * g * tti0**2)
    return (tti, vi)

def burntime(delta_v):
    F = vessel.available_thrust
    Isp = vessel.specific_impulse * 9.82
    m0 = vessel.mass
    m1 = m0 / math.exp(delta_v/Isp)
    flow_rate = F / Isp
    burn_time = (m0 - m1) / flow_rate
    return (burn_time)

def decel_dist(v):
    g = vessel.orbit.body.surface_gravity # not needed, already in impact calc
    a = (vessel.available_thrust / vessel.mass / 0.981)
    t = v / a
    dist = (0.5 * a * (t**2))
    return dist

1. I've calculated a free fall, how long until I hit the ground, and how fast.
2. "How fast" was assumed to be the required dV for the suicide burn, burn time derived from that.
3. Finally, another free fall (though in reverse) for distance traveled during burn time under my vessel's acceleration. That was the altitude where I lit the engines.

Of course, as soon as I light the engines, that changes the time until impact and thus everything else. If applied to long burns, that approach will be horribly lossy. I only needed it as a final flourish, though, and it was definitely good enough for that purpose.  Here's how it played out.

 

2 hours ago, Ypsilon said:

If the tick rate of my script was too slow, I'd expect to have more inconsistent results.

That was more philosophical than aimed at your problem at hand. The precision of any action cannot be better than however far you're moving between updates -- if you want to stop at altitude X +-0.1m, but start the maneuver while moving 12m between updates, you will have to refine your solution along the way.

 

 

Edited April 30, 2020 by Laie

[CBase]

On 4/30/2020 at 9:15 PM, Laie said:

That was more philosophical than aimed at your problem at hand. The precision of any action cannot be better than however far you're moving between updates -- if you want to stop at altitude X +-0.1m, but start the maneuver while moving 12m between updates, you will have to refine your solution along the way.

Actually you should plan as well for differences in surface altitude. Predicting the true landing position for the suicide burn is way harder than the speed timing. But without perfect position you do not have perfect altitude... So as Laie wanted to express: Add some safety margin for things that you do can not know exactly when starting the burn.

On 4/30/2020 at 6:56 PM, Ypsilon said:

I've been testing either at Kerbin or the Mun.

Kerbin has atmosphere, the Mun not. The difference is huge even without parachutes, even from 1000m down. No single equation will get you landed on both.

On 4/30/2020 at 9:15 PM, Laie said:

tti0 = math.sqrt(2 * alt0 / g) # time-to-impact, free fall from starting alt

Unfortunately the equations are not that easy even for the Mun as the problem is 2d: You need to kill horizontal speed and gain vertical speed, both are connected by the vehicle angle which constantly changes.

However you are right, basically it is a rocket equation reverse, so iterative solutions or known solvers might help.

I did recently touch this topic by MechJebs Landing scripts, but only changed final descent, which kicks in once your are falling mostly straight.