Ascent Optimization as an Optimal Policy Problem

[K^2]

The recent discussion on optimal ascent problem has reminded me that I saw some interesting ways of dealing with optimization problems in finance that I wanted to try out. I've put together some basic Mathematica code as a proof of concept, and actually got nice results out of it that I wanted to share.

The idea behind optimal policy approach is that you start out with policy function, which is the desired outputs based on current state. In case of rocket's ascent in KSP, the outputs are going to be the thrust and attitude of the rocket, while the state will be the altitude and horizontal and vertical velocities. Once we have a policy function, we can compute state evolution, from which we can compute an integral that we wish to optimize. Again, in KSP we would compute trajectory the rocket follows under given policy, and use it to estimate delta-V required to reach orbit. Once we have that, we can vary the policy function to try and optimize the result.

There are two major obstacles in doing this in KSP. First, we have three input variables and two output variables. That is a huge space of available policy functions. Second is the fact that we need to establish boundary conditions, which is tricky as stated. Fortunately, both of these are resolved if we re-formulate the problem. Instead of looking at velocity as the state, we'll consider it to be the desired output. After all, for a given ascent profile, there will be a unique horizontal and vertical velocity for each altitude, and thrust can be computed from the way these vary with altitude.

Final obstacle is the fact that we are dealing with a continuous range of variables. The simple solution is to use parametrized polynomial approximations and then the search space for optimization consists only of the parameters.

Lets start with the 1D problem of vertical ascent, which we know analytic solution to, and state it explicitly.

Given altitude h, we wish to find optimal ascent velocity v(h) that is most fuel efficient. In other words, we wish to minimize integral over T(t)dt, where T is thrust and t is time. We shall now compute the necessary thrust in units of rocket mass.

T(t) = g + k(h) v² + dv/dt

Here, k(h) is drag coefficient, k(h) = k₀e^-h/H for some scale height H. Now we shall take dv/dt using chain rule.

dv/dt = (∂v/∂h)(dh/dt)

I shall denote ∂v/∂h as v'(h), and, of course, dh/dt = v(h).

T(h) = g + k₀e^-h/Hv²(h) + v'(h)v(h)

Of course, we are interested in integral over T(t)dt, but noting that dt = dh/v(h), we can re-write it as an integral (T(h)/v(h))dh, which we can integrate over having only the policy function v(h) at hand.

Finally, we need to specify the parametrized form for v(h). For simplicity, I've used Legendre polynomials, L_n(x).

v(x) = v₀ + (v_max - v₀)(a L₃(x) + b L₅(x) + c L₇(x) + (1 - a - b - c)L₁(x))

This means that velocity will vary from v₀ to v_max for x = h/h_max from 0 to 1. Since this is to serve as a proof of concept, we'll assume that v₀ and v_max are the terminal velocities at h = 0 and h = h_max respectively. If the approach is valid, minimizing integral of T(t)dt with respect to parameters a, b, c should produce terminal velocity curve as the policy function. And, of course, it does.

The result of optimization is in gray, while the terminal velocity plot is in green.

Now, for a more useful problem. The policy is replaced with a pair of functions v_x(h) and v_y(h). The thrust function T(h) is far more complex, but is derived in the similar way. The integral we optimize is still (T(h)/v_y(h))dh. But the best part is our ability to set initial and final velocities. At h = 0, both v_x and v_y are zero. At h = 80km, v_y is zero, while v_x is orbital. Please note, I did ignore Coriolis force and rotation of the planet. Other parameters are set to Kerbin from sea level. The parametrization for v_x(h) is very similar to v(h) above. In fact, it's the same, with v₀ = 0. Parametrization for v_y is way more complex, because I needed an asymmetric function with zeros at both ends.

I'm not going to go through all of the results and potential problems, mostly because I want to build a finished version that accounts for all of KSP physics first, but I do have a rough result for the ascent profile. This is the angle of rocket's velocity with respect to horizon from sea level to the 80km orbit.

As expected, it is 90° off the pad, and 0° in orbit. The gravity turn does, in fact, begin immediately, but it is more gradual than some people have suggested. Of course, the TWR for this ascent peaks at roughly 3:1. I have not considered caps on thrust yet.

The great thing about this approach is that it is light weight enough to be used on the fly. Once I have all the gremlins worked out and have it working with thrust cap and variable I_SP, it'd be nice to see if this can be written into a plugin that does automatic ascents in KSP.

Edited August 11, 2014 by K^2

[Z-Man]

Very nice. I see a small problem at the end of orbital insertion, though, probably not big enough to cause worry and easy enough to fix if it turns out to be relevant. Since you want v_y = 0 at the end, you naturally have dh/dt = 0. If you have any x-acceleration at the end (and you should have), you have

dv_x/dt > 0

and in the parameter target function

dv_x/dh = dv_x/dt / (dh/dt)

Which would be infinite and the polynomial parametrization can't deliver that.

The problem should manifest itself in the T(h) function: if the above is right, it should go to zero for h = h_max, which may not be terribly optimal.

Easy solution: Add to the v_y target function

d * ((h_max-h)^(1/3) - h_max^(1/3))

d being a new parameter. I hope the exponent 1/3 is right, it should be what allows for a finite acceleration at the end.

[indroth]

T(h) = g + k₀e^-h/H + v'(h)v(h)

Shouldn't there be a v² in there? Or am I missing something?

[Ralathon]

Looking forward to the finished method. What kind of steering losses are you looking at with that derived curve? I imagine people who use the Ferram Aerospace Research addon wouldn't want to see angles more than a few degrees from prograde due to aerodynamic stress.

[Laie]

What kind of steering losses are you looking at with that derived curve? I imagine people who use the Ferram Aerospace Research addon wouldn't want to see angles more than a few degrees from prograde due to aerodynamic stress.

Looks reasonably gentle to me. 75° by 10km, 55° by 20km -- just try it.

[K^2]

Looking forward to the finished method. What kind of steering losses are you looking at with that derived curve? I imagine people who use the Ferram Aerospace Research addon wouldn't want to see angles more than a few degrees from prograde due to aerodynamic stress.

No idea. Though, there are a couple of places where thrust vector differs significantly from velocity vector, so I think this wouldn't work well in FAR. I'd have to take a look at FAR model and take it into account to get a somewhat different profile. As it is, this is just for vanilla drag.

The problem should manifest itself in the T(h) function: if the above is right, it should go to zero for h = h_max, which may not be terribly optimal.

Hm. It has no problem getting "infinite" thrust off the pad, to quickly build up to terminal velocity. Shouldn't it work the same way for h_max by your logic? The other thing is that rocket is directed to burn mostly down at the top. So it might make sense for throttle to ease off. It seems that situation here is similar to a landing burn, where you start burning towards the planet/moon initially in the beginning. (Constant altitude descent method.)

Shouldn't there be a v2 in there? Or am I missing something?

Indeed, there should be. I'll edit it in.

[Z-Man]

It has no problem getting "infinite" thrust off the pad, to quickly build up to terminal velocity.

You built that in explicitly with the v₀ parameter for v_y. If you force that to zero, If you take that out (and I think you should, now that you brought attention to it) you only get finite thrust at the start. The same is true at the end: You can only get the behaviour out of the process that you explicitly allow for in the parametrization.

But yeah, if you only optimize for total dv spent without TWR limitations, easing the throttle down at the end may be a good thing, Oberth and all. And even with TWR limitations, Kerbin has such an insane atmosphere thickness to radius ratio that just leaving the atmosphere gives you more than enough time to push to orbit efficiently even with lower thrust.

Maybe h is not the best argument to parametrize after? Maybe the total energy per mass ( g h + (v_x² + v_y²)/2 ) or some modification would be better? If the rocket is not increasing that one steadily, it's doing something wrong

[CjStaal]

This is brilliant. I would hope you could formulate one that would work with FAR as well though, since that's more realistic, there might be a solution already online for earth based trajectories that you can derive off of for kerbal trajectories.

[K^2]

You built that in explicitly with the v₀ parameter for v_y.

Not for the 2D case. In the 2D case, vertical velocity is precisely zero both at sea level and at target altitude. Which is precisely what you want from an ascent that takes you to circular orbit. Thrust is still extremely large at takeoff. And that's precisely because initial velocity is zero. In vertical ascent case, it's well known that thrust needs to be infinite to bring up ascent velocity to terminal instantly, thereafter settling on TWR of 2:1. The results I'm getting agree with that. Rocket blasts straight up, achieves terminal velocity, then starts to pitch over to initiate gravity turn.

P.S. Lets forget all the terms that obviously aren't going to contribute to an infinite result. T_y(h) = v'(h)v(h). Yes, v(h) goes to zero as h goes to zero. But v'(h) can, and does, diverge. Precisely at zero, the result is actually undefined. But there is a limit. And lim T_y(h) as h goes to 0 can be zero, finite, or infinite. In this particular case, it does happen to diverge to infinity. At h very small, but above zero, I get a numerical result that's finite, but very large.

Edited August 11, 2014 by K^2

[Z-Man]

Not for the 2D case. In the 2D case, vertical velocity is precisely zero both at sea level and at target altitude.

And you're also absolutely correct about getting a divergent thrust at liftoff even then from the term that gives v_x approximately proportional to h at the start. I hadn't noticed. That is definitely going to bite you in the rear later when you want to add TWR limits, though. You have no term there that would produce constant thrust (>1, that is) at the start. Choose your parametrization wisely, is all I'm saying.

Ah, right. I misread what you wrote in your original post then when you said you set v₀ to 0 specifically for v_x. And didn't you also say TWR peaks at 3?

[K^2]

TWR starts off infinite, drops to a value of approximately 2 during early ascent, but then climbs to a bit past 3 as the rocket settles into gravity turn. Once well into gravity turn, TWR starts to drop again. It's somewhere between 1 and 2 through most of the ascent.

I think I'll have to go full Lagrange on the max throttle. I can't think of any good way to enforce the limit other than doing formal constraint. So the integral will be [T(h)/(I_SP(h)v_y(h)) + λ(h)(T(h) - T_max)]dh. And I'll just have to solve for λ(h) at every integration point.

[mhoram]

@K^2

Very interesting approach. Lookig forward so see your progress.

This is brilliant. I would hope you could formulate one that would work with FAR as well though, since that's more realistic, there might be a solution already online for earth based trajectories that you can derive off of for kerbal trajectories.

A major factor of this theoretical ascent optimization is drag.

Stock drag is very simple for rockets of all buildforms and easily put into a formula.

Far's drag calculation way more complex and depends much more on the form of the rocket and unil someone will come up with a reasonable way to theoretically calculate FAR-drag in a simulation, there will not be such a thing.

Maybe one could use the FAR-sourcecode itself for such an Optimizer. ;-)