I need someone help me do some math for launch optimization

[SaturnV]

Some equation can be found here: http://wiki.kerbalspaceprogram.com/wiki/Atmosphere

The problem: What speed you need to keep at certain altitude to maximize your vehicle efficiency?

Detailed Description: In space, just install MJ and watch how much delta-V you remain, that's it. But on Kerbin(mostly where you launch from, we will discuss other planet with atmosphere later, maybe), the air resistance is a big problem to optimize your launch, as fluid mechanism is one of the most difficult science on earth.

But in KSP, we can simplify the most complex part: the drag coefficient & cross-sectional area. Let's consider a single-stage rocket, both parameters is constant during the whole launch. I merge that with some other constant into one, called constant C. That's not very important in our further discuss.

As you speed up in dense air, you slow down more, that drag is related to square of you velocity so it increase dramatically if you're too fast.

But if you are fast, you can raise up quickly. The air density is decreased exponential which results in a lesser drag force, since Fd is related to altitude. What's more, faster means higher thrust output means faster fuel consume rate means faster mass decreased, that's an additional acceleration boost.

So, is there an optimal point can balance the air drag and altitude gain, with certain vehicle with given mass, air drag factor, thrust strength?

Is there a speed map that you could follow, say "you should not exceed the speed 165m/s at altitude 4500m, with this vehicle"?

[PakledHostage]

Have a look at Closette's Mini-challenge: max altitude with this supplied spacecraft thread from about a year and a half ago. That thread documents the great deal of effort that several members of this forum undertook to study the one dimensional (i.e. vertical climb only) version of this very problem. Many of the key contributors in that thread went on to collaborate on MechJeb's initial development. As far as I know, that thread is the origin of the "climb at terminal velocity" and "design for TWR=2" rules of thumb that everyone on these forums cites as gospel.

More recently, Tarmenius' Launch Efficiency Exercise thread attempted to find the optimal two-dimensional launch trajectory (i.e. vertical and horizontal acceleration into an orbital trajectory) by "crowd sourced" trial and error.

Both threads were specific to two unique spacecraft designs, but the principles discussed are still relevant to more general designs.

Edited August 15, 2013 by PakledHostage

[K^2]

Optimal climb rate is, indeed, terminal velocity. For constant density, the proof is trivial. The thrust, in units of ship mass, is equal to F = g + kv². The amount of time needed to rises to some altitude h is t = h/v. In that time you consume some quantity of fuel proportional to Ft, so your goal is to find maximum of (g + kv²)/v. Differentiating that with respect to v and setting derivative to zero yields k - g/v² = 0. Or g = kv². Note that this is the same exact condition as setting F = 0, meaning you just found terminal velocity, except going upwards instead of down. Substituting into original equation we have F = 2g, meaning thrust is twice the weight.

Now, the problem is actually more complicated because air density is decreasing exponentially, so the value k is varying as k(h) = k0 exp(-h/H), where H is scale height. Math for this problem is extremely hairy and involves calculus of variation. I have solved this problem, and if anyone here wants to see the proof, I can outline it. But the result is exactly the same. Climb must be performed at terminal velocity for given altitude. In other words, v(h) = v0 exp(h/(2H)), where v0 is terminal velocity at sea level. Now, however, thrust must compensate for drag and acceleration. By the time acceleration is significant, however, you should already be well into gravity turn, so it turns out that TWR of 2 is still a good one to shoot for.

Two dimensional problem involving gravity turn is extremely complex, and I have not found a solution. I have verified a number of "rules of thumb" to not be true optimums, but even solving for optimal ascent numerically has proven to be beyond the "simple" methods.

Edit: There is an interesting consequence for all of this. Again, I'm not going to get into all the math, but you can estimate the extra delta-V required to climb out of atmosphere as 4gH/v0. Again, H is scale height and v0 is terminal velocity at whatever height you start out from. This means that the total minimum delta-V required to reach stable circular orbit is roughly dV = 4gH/v0 + Sqrt[MG/(h+R)]. Here R is radius of the planet, M is its mass, G is the gravitational constant, and h is the altitude of your final circular orbit. For Kerbin, that works out to 4,128m/s, which is very close to empirically established value.

Edited August 15, 2013 by K^2

[Diche Bach]

Fascinating topic. The parameters re unfamiliar to me, but the concepts are quite intriguing.

Couple ideas/questions that occur to me as a result of (a) watching Scott Manley's video on the Orion nuclear explosive pulse propulsion spacecraft; ( reading about the apparent solution in recent years to the "Pioneer Anomaly;" © just a general sense of how aerodynamics works.

1. Air resistance or drag are largely a function of the shape of the vehicle right? For example, the Concorde with its extra long and pointy nose. Why aren't rockets that employ these kinds of features used? Would it not make enough difference? Would it produce other problems? Would it not be possible to make a nose cone that was sufficiently rigid?

2. Air density is going to be a function of heat right? What about (in addition to more aerodynamic nose) putting some lasers on the nose of the vehicle that project directly up into the prograde direction of the vehicle, thus (maybe!?) creating a heated 'column' of air immediately above/in front of the craft thus reducing drag and improving efficiency? Seems insane, but that's the luxury of being a naïve social scientist forumite on a gaming forum with all you physical scientists around to tell me how hilarious such ideas are

3. Air gets compressed when something pushes on it, so I would think that, to some extent (depending on the bluntness of the craft) the vehicle itself causes the resistance in front of it to increase as it ascends. However, if that pressure upwards were not applied continuously but rather in some oscillating pattern, I would think that instead of simply 'stacking up' the air into a compressed space it might actually push some of it out of the way, at least momentarily. In large part I reckon this is what a highly streamlined pointed nose cone does. But in addition to this, what about some method for actually pushing the air directly above/ in front of the vehicle out of the way? Ideas here could be: a nose-cone that rapidly pumps up and down (aka the pumping pressure plate on the Orion is what is giving me this crazy idea), or a series of small explosive devices that launch out in front of the vehicle and blast air out of the way.

Thank you for reading my absurd, naïve and insane anthropological theories on aeronautics. I appreciate all of your humor in pointing out just how crazy these ideas are!

[PakledHostage]

1. Air resistance or drag are largely a function of the shape of the vehicle right? For example, the Concorde with its extra long and pointy nose. Why aren't rockets that employ these kinds of features used? Would it not make enough difference? Would it produce other problems? Would it not be possible to make a nose cone that was sufficiently rigid?

The Concorde is optimized for supersonic cruise, while launch vehicles are optimized for maximizing payload delivered into orbit. Given that launch vehicles climb through the thickest part of the atmosphere at sub-sonic speeds, their optimal aerodynamic shapes are different than that of supersonic aircraft.

2. Air density is going to be a function of heat right? What about (in addition to more aerodynamic nose) putting some lasers on the nose of the vehicle that project directly up into the prograde direction of the vehicle, thus (maybe!?) creating a heated 'column' of air immediately above/in front of the craft thus reducing drag and improving efficiency?

I would imagine that the mass and cost of the lasers, together with the energy requirements to run them would outweigh the benefit.

3. Air gets compressed when something pushes on it, so I would think that, to some extent (depending on the bluntness of the craft) the vehicle itself causes the resistance in front of it to increase as it ascends. However, if that pressure upwards were not applied continuously but rather in some oscillating pattern, I would think that instead of simply 'stacking up' the air into a compressed space it might actually push some of it out of the way, at least momentarily. In large part I reckon this is what a highly streamlined pointed nose cone does. But in addition to this, what about some method for actually pushing the air directly above/ in front of the vehicle out of the way? Ideas here could be: a nose-cone that rapidly pumps up and down (aka the pumping pressure plate on the Orion is what is giving me this crazy idea), or a series of small explosive devices that launch out in front of the vehicle and blast air out of the way.

The process of compressing the air and pushing it out of the way happens by default. The speed of sound is the speed at which an infinitesimal pressure wave propagates through the air. Stronger pressure waves (such as blast waves from explosions) propagate at higher speeds. A body moving through the air “notifies the air ahead of it that it is coming†by means of pressure. Supersonic objects develop pressure waves (shock waves) around themselves that are strong enough to propagate at the speed that the object is moving. Air passing through a supersonic shock will increase in pressure, density and temperature. The shock can be either normal or oblique. Oblique shocks are weaker than normal shocks and cause less drag.

Again, the faster the object is moving, the stronger the shock must be. The Concorde’s pointy nose was chosen to ensure that the Concorde forms a (lower drag) oblique shock while flying at supersonic speeds. Missiles also have pointy nose fairings like the Concorde, but the aerodynamic fairings on launch vehicles would be optimized to minimize drag over the full range of speeds that they encounter during the flight phase where they are required.

[SaturnV]

Have a look at Closette's Mini-challenge: max altitude with this supplied spacecraft thread from about a year and a half ago. That thread documents the great deal of effort that several members of this forum undertook to study the one dimensional (i.e. vertical climb only) version of this very problem. Many of the key contributors in that thread went on to collaborate on MechJeb's initial development. As far as I know, that thread is the origin of the "climb at terminal velocity" and "design for TWR=2" rules of thumb that everyone on these forums cites as gospel.

I finished it, but I need more time to catch on. That's a great thread, so do you. Thanks a lot dude.

[tavert]

Not sure if my derivation matches K^2's as I don't think he's posted it anywhere yet, but I went through this a few weeks ago here:

http://forum.kerbalspaceprogram.com/showthread.php/43706-Rocket-Acceleration-g?p=573727&viewfull=1#post573727

[Mrwood]

http://sbir.nasa.gov/SBIR/successes/ss/8-048text.html NASA plasma drag reduction research.

The Concorde is optimized for supersonic cruise, while launch vehicles are optimized for maximizing payload delivered into orbit. Given that launch vehicles climb through the thickest part of the atmosphere at sub-sonic speeds, their optimal aerodynamic shapes are different than that of supersonic aircraft.

I would imagine that the mass and cost of the lasers, together with the energy requirements to run them would outweigh the benefit.

The process of compressing the air and pushing it out of the way happens by default. The speed of sound is the speed at which an infinitesimal pressure wave propagates through the air. Stronger pressure waves (such as blast waves from explosions) propagate at higher speeds. A body moving through the air “notifies the air ahead of it that it is coming†by means of pressure. Supersonic objects develop pressure waves (shock waves) around themselves that are strong enough to propagate at the speed that the object is moving. Air passing through a supersonic shock will increase in pressure, density and temperature. The shock can be either normal or oblique. Oblique shocks are weaker than normal shocks and cause less drag.

Again, the faster the object is moving, the stronger the shock must be. The Concorde’s pointy nose was chosen to ensure that the Concorde forms a (lower drag) oblique shock while flying at supersonic speeds. Missiles also have pointy nose fairings like the Concorde, but the aerodynamic fairings on launch vehicles would be optimized to minimize drag over the full range of speeds that they encounter during the flight phase where they are required.

[K^2]

Not sure if my derivation matches K^2's as I don't think he's posted it anywhere yet, but I went through this a few weeks ago here:

http://forum.kerbalspaceprogram.com/showthread.php/43706-Rocket-Acceleration-g?p=573727&viewfull=1#post573727

You are looking at mass changes. I basically disregarded that because in KSP drag is proportional to the mass, so everything can be done in units of ship's mass.

What I was working on is accounting for dÃÂ/dt terms. Your derivation simply assumes that ÃÂ is constant, but it really isn't, and it would throw a major monkey wrench in your computations, because then you can't just say, "What is optimal v?" You have to say, "Which function, v(t), optimizes total fuel consumption." And that's Variation Calculus problem. Fortunately, some really smart people found some really nice shortcuts there.

[tavert]

Care to share, if you've been deriving something more detailed? I wasn't necessarily considering rho as constant, but rather a variable that you don't have control over if you are trying to determine the (locally, greedy) best v at a given altitude.

I've looked at the optimal control literature on this problem and haven't found many detailed analytical derivations. There's a good paper by Rugescu that goes through the history of developments on the problem, but I wish it was more detailed in the math and most of the references are either in German or difficult to find.

I've been making headway on coding up a numerical solution of the full problem including gravity turn, but actual work keeps interfering.

Edited August 20, 2013 by tavert

[K^2]

Well, the basic idea is that total fuel used (ignoring I_SP changes due to pressure) is proportional to trust. Now, there is choice in what we are optimizing with respect to, but if an integral over F(t) from some t₁ to t₂ can be reduced with minor variations in v(t), then you haven't done the best possible job. Now, given altitude as h(t), velocity is h'(t) and acceleration is h''(t). (Just making note of notation.) Then in units of ship's mass, we can write.

F(t) = h''(t) + g + k(t)h'²(t)

Here k(t) is the drag coefficient given by k(t) = k₀ exp(-h(t)/H). This is proportional to ÃÂ in your computations.

So what we want to optimize is a functional S[h(t)] = ∫F(t)dt. Any function h(t) that satisfies such an optimization will also satisfy the Euler-Lagrange Equation. For a function F(t) which depends on h(t), h'(t), and h''(t), the following is satisfied.

∂F/∂h - (d/dt)∂F/∂h' + (d²/dt²)∂F/∂h'' = 0

In other words, whatever the optimal ascent is going to be, it's going to satisfy the following equation.

(d/dt) (∂k(t)/∂h h'²(t) - 2 k(t)h'(t) = 0

Two things should be noted right away. First, if k(t) was constant, we'd have v(t) is constant as solution. Second is that this doesn't depend on g and allows for infinitely many solutions. For instance, v(t)=0 is clearly a solution. This isn't THE optimal one, just a constraint on possible optimal solutions. But it's a very useful constraint. Taking time derivative, we get the actual differential equation for velocity.

2k₀ exp(-h(t)/H) h''(t) = k₀ exp(-h(t)/H) h'²(t)/H

Which, thankfully, simplifies.

2H v'(t) = v²(t)

That isn't so bad, and general solution is easy to find.

v(t) = -2H/(t+C)

Where C can be any constant. On one hand, this is a solution. Give me v(t=0) and I can find C that will make this work. But that's not the interesting part. Consider the velocity that satisfies the following equation.

v(t) = v₀ exp(h(t)/(2H))

You will realize that if v₀ is terminal velocity at the surface, the above is terminal velocity at altitude h(t). This differential equation can also be solved, and with h(t=0) = 0 it yields the following.

h(t) = -2H ln((2H-tv₀)/(2H))

What's more interesting is that I can find the velocity of the above profile.

v(t) = 2H v₀/(2H - tv₀)

And this, in turn, satisfies the earlier equation with C = -1/v₀.

In other words, ascent at terminal velocity through exponentially thinning atmosphere is locally optimal. And therefore, being the optimal solution for constant density is the overall optimal way to ascend through exponentially thinning atmo.

Edit: Fixed ∂F/∂h term.

Edited August 20, 2013 by K^2

[tavert]

Interesting, thanks. One quick clarifying question, it looks like your F(t) is thrust divided by rocket mass? So by integrating that, you're minimizing something slightly different than total fuel use, are you not? I suppose you're minimizing delta-V which has a one-to-one correspondence with fuel, so you'll get an equivalent minimizer.

Edit: I think you may have missed a factor of two (or cancelled only on one side of an equation instead of both):

(d/dt) 2 k(t)h'(t) = 0

(d/dt) 2 k0 exp(-h(t)/H) h'(t) = 2 k0 exp(-h(t)/H) h''(t) - 2 k0 exp(-h(t)/H) h'(t)^2 / H

Edited August 20, 2013 by tavert

[K^2]

Heh. Yeah, I missed something. Sorry, I'm actually doing all of this in Mathematica, which lets me skip steps. What's missing is the ∂F/∂h term. I'll make an edit.

Hm. This is a much messier thing to differentiate than it looked. I'm really glad I wasn't doing this on paper. I'm feeling a bit lazy to fill in the steps, but hopefully, you can figure it out. The equations in the main post should be right now.

Edited August 20, 2013 by K^2

[tavert]

Looks right now. It's not too terrible to differentiate by hand, just very easy to miss a term, a sign, or a factor of 2 somewhere. I've got Mathematica open as well and now trying to play with the variable-gravity solution, things didn't cancel as nicely there and DSolve couldn't get it. Only a couple percent variation on Kerbin over the altitude ranges this analysis matters for, so constant's good enough here.

And figures that Euler-Lagrange does the trick here and you get a nice little differential equation. I haven't used that method since the last dynamics course I took about 7 years ago. Always with the brachistochrone problems and involute pendulums, not nearly enough rockets in those textbooks that I can remember.

Edited August 20, 2013 by tavert

[tavert]

Thinking about this a little more, you've shown that terminal velocity satisfies the Euler-Lagrange equation for this problem. But so does any other velocity profile of the same form with a different v0, right? We need another few steps to show that the particular best v0 we want is actually terminal velocity at 0 altitude, sqrt(g / k0), don't we?

[K^2]

Technically, yes. I was a bit lazy about doing this strictly analytically, and I'm not sure how hard it'd be to actually do. If you take h to infinity, even thought time required for ascent is finite, fuel requirement diverges because acceleration goes to infinity as h does. So you probably need to set up cutoff height, optimize for v₀, and only then take the limit.

Alternatively, one can be satisfied with the fact that numerically v₀ = v_t works. That's the attitude I chose. Call me lazy if you must.

[tavert]

Well if I'm going to do the problem numerically, may as well go all the way and do the whole problem with a gravity turn, Isp varying with altitude, hard atmosphere cutoff, the whole deal.

I was trying basically what you describe (and in KSP there's a nice meaningful cutoff altitude at H*ln(10^6)), though the results don't make much sense. I integrated F(t) from 0 to tf after substituting h(t) = -2H ln((2H - t v0)/(2H)).

Then substituted tf = 2H (1 - exp(-hf / (2H))) / v0, differentiated the result wrt v0, and found a zero at:

v0 = sqrt(2 g H / (exp(hf / (2H)) + 2 H k0))

I'm not sure that's right...

To put different v0 solutions on an equal footing for comparison purposes, I think you'd need to include a coast phase and solve for h_meco(v0, h_apo).

[K^2]

I tried to get a numerical solution with gravity turn. I wasn't able to find anything resembling an optimum. If you have good ideas on that, I'd be interested.

[tavert]

You're not the only one, the discussion of getting the full nonlinear gravity turn problem with drag and variable throttle and everything has come up several times. I think AlterBaron was going to have a look at using PSOPT (http://www.psopt.org/) to try a pseudospectral method, they have a rocket-launch example problem that I think he was going to start by trying to adapt for KSP fizziks.

I was going to try a local discretization method, simple RK4 or maybe something implicit with a short time step and suitable initial and final conditions applied as equality constraints (or upper bound = lower bound for those variables). I was going to formulate the nonlinear programming problem including dynamics and all input variables for every time step, using some optimization modeling language that would do automatic differentiation for me, like AMPL or OpenOpt, or a new Matlab toolbox that I've been developing in my research. I'd use a large-scale sparse nonlinear interior point solver, most likely Ipopt, to find local minima, and if it converges to anything interesting then I'll start playing around with multi-start to see how many different local minima it can get stuck at, and how different in performance they are.

I've started and made a bit of progress (I can send you my Matlab stuff and explain how to use it if you want to look at it), but my advisor got back in town last week and I'm behind in other areas of my research, I have code for linear systems I need to write that won't be useful for KSP. And I haven't attempted to convince my advisor that I find KSP-related research more interesting than that other code...

Edited August 21, 2013 by tavert

[K^2]

Ha, ha. Yeah, I know the feeling. Well, I've been wanting to try something along the lines of a genetic algorithm and see if I can "evolve" a near-optimal solution. The "genome" would consist of adjustments to attitude and throttle position. It should work, but time...

As for the RK method, you shouldn't need to go for RK4. Just do a Velocity Verlet. It's equivalent to a mid-point RK2, much easier to implement, and should be sufficiently precise for the given problem.

[tavert]

Are you gmcnew on Github by any chance? I'm not a fan of genetic algorithms, but they have their uses. Here I think it would be great to use a genetic algorithm or some other randomized algorithm (whether it be simulated annealing, particle swarm, neural networks, whatever you feel like) as an outer-level iteration to generate initial guesses for a Newton-type local nonlinear programming solver.

Can you do velocity Verlet with drag? Acceleration depends on velocity, so the velocity update becomes implicit. Implicit's not a big deal since I'm using an optimization algorithm that can handle nonlinear equality constraints, but then I may as well just use trapezoidal. I have no idea what kind of sensitivity to expect with respect to time step length, discretization method, or number of time steps, so I'll play around with things once I get something working. RK4's not that much more difficult than anything else in principle, give me a Butcher tableau and I can do it... I just never feel confident if I start with Euler.

[K^2]

You are technically half a time step off with drag in Verlet, but in practice, it's not a big deal unless you have huge swings in velocity. And no, I'm not on Github. I don't usually go for genetic algorithms. It just somehow feels right for this one.

[SaturnV]

I can't follow you the whole way, but your calculation looks rigorous. I was lost in the Euler-Lagrange part, other derivation seems reasonable. Thanks for your work, I should noticed that this is far beyond me, even when the answer is in front of me, I can't understand it fully.

[alterbaron]

Can you do velocity Verlet with drag?

Here's the code I use to do velocity verlet with velocity dependent forces.

function [r v t] = vel_vert_term( F, m, r0, v0, dt, term_cond )
% VEL_VERT_TERM Integrate the equations of motion for a particle of mass
%     m being acted upon by a vector force F. Terminate when term_cond
%     is true.
%     
%     F - Function which returns the vector force acting on the particle,
%         as a function of position and velocity.
%     m - Particle mass
%     r0 - Initial position
%     v0 - Initial velocity
%     dt - Integration time step

t = 0;
r = r0;
v = v0;

a = @(rin, vin) F(rin, vin)./m;

firstrun = true;

while firstrun || (term_cond(r,v,t) == 0)
    r_old = r;
    v_old = v;
    a_t = a(r_old, v_old);
    r = r_old + v_old*dt + 0.5*a_t*dt.^2;
    v_est = v_old + 0.5*dt*(a_t + a(r,v_old+dt*a_t));
    v = v_old + 0.5*(a_t + a(r,v_est))*dt;
    t = t + dt;
    firstrun = false;
end

end

Edited August 23, 2013 by alterbaron

[tavert]

Interesting, though I'd also want to include a mass state as it simplifies the thrust constraint. I could normalize things in terms of delta-V and TWR, but then the max thrust constraint is messy.

Any luck with PSOPT?

Edit: Have you done an error analysis of that method (or did you get it from some reference that did)? You could do RK3 for the same number of function evaluations. In terms of your code, RK3 would look like:

while firstrun || (term_cond(r,v,t) == 0)
    r_old = r;
    v_old = v;
    a_t = a(r_old, v_old);
    r_half = r_old + 0.5*v_old*dt;
    v_half = v_old + 0.5*a_t*dt;
    a_half = a(r_half, v_half);
    r_est = r_old - v_old*dt + 2*v_half*dt;
    v_est = v_old - a_t*dt + 2*a_half*dt;
    a_est = a(r_est, v_est);
    r = r_old + (v_old + 4*v_half + v_est)*dt/6;
    v = v_old + (a_t + 4*a_half + a_est)*dt/6;
    t = t + dt;
    firstrun = false;
end

Edited August 23, 2013 by tavert