I need someone help me do some math for launch optimization

[alterbaron]

I did some numeric simulation, the result shows that stick to terminal velocity is practically optimal.

Nice simulation. I'm actually quite surprised at how little that thrust oscillation affected the final result.

Another update:

I've been able to solve for an optimal ascent trajectory for a body without an atmosphere, which is kind of cool in its own right.

Here's the results from a test ascent from the surface of Kerbin (with the atmosphere removed):

Altitude wrt. Time

Thrust wrt. Time

Thrust Angle wrt. Altitude

The target orbit was set to 400,000m altitude. The ship model used was the same as in the first test, but the engine Isp was set to 800s so it would have enough fuel to orbit.

There's a lot of noise in the "thrust profile" function whenever the thrust is zero. This is because there's no preferred direction to point your ship when the thrust is zero. (Is there a way to eliminate this noise?)

I've had partial success with ascent in the atmosphere. The program definitely figures out to do a gravity turn (beginning around 10km) and to do a burn at periapsis, but the output is still quite messy:

Atmospheric Ascent Thrust

Atmospheric Ascent Profile

I'll play with the program and model some more to see if I can get a better result.

Edited September 1, 2013 by alterbaron

[K^2]

There's a lot of noise in the "thrust profile" function whenever the thrust is zero. This is because there's no preferred direction to point your ship when the thrust is zero. (Is there a way to eliminate this noise?)

You have to make turns cost something. Like a very small amount of fuel, in this case, or gong for a complete solution with RCS fuel consumed and also optimzied for. That's the only way I know.

[tavert]

Or change your coordinates to horizontal thrust and vertical thrust. Then the origin is unique.

[alterbaron]

Or change your coordinates to horizontal thrust and vertical thrust. Then the origin is unique.

That's a great idea (and a simple change to make).

You have to make turns cost something. Like a very small amount of fuel, in this case, or gong for a complete solution with RCS fuel consumed and also optimzied for.

Might be worth looking into, but might make the already hard to solve problem harder to solve.

[alterbaron]

Taking tavert's advice, I converted over to using horizontal and vertical thrust (rather than thrust and thrust angle).

With a little tweaking, I've been able to coax some nice results out of PSOPT.

I input a model of the SSTO from the launch efficiency challenge to put the program through it's paces.

It was able to solve for a launch profile to orbit that challenging craft. However, it used 0.1 tons more fuel than the craft has available.

Perhaps with more computation time this will improve? It's still a good result regardless.

Here are the results:

Pitch as a function of altitude: Thrust Angle Chart

Thrust as a function of altitude (there's a bit of noise visible, might account for extra fuel consumption.): Thrust Magnitude Chart

Velocity compared to terminal velocity: Low-Altitude Velocity (with Terminal Velocity for Comparison) Chart

Final craft mass: 3.2066 tons (craft dry mass = 3.3 tons).

There's still some improvement possible, and progress is being made!

EDIT: I managed to get a much better result that eliminated the noise in the thrust function above.

Interesting things to note (do these make sense?):

• Ship stays throttled down to maintain a velocity lower than terminal between 10-15km when it could have throttled up again.
  
• Ship throttles up past 15km, then throttles down significantly as it pitches towards the horizon.
  
• Continuous application of (light) thrust in the upper atmosphere.
  

Final craft mass was 3.2148 tons, so still a little short on fuel. Still, a margin of <0.1 ton of fuel is pretty OK.

Optimization took 12 minutes (18k solver iterations) using Legendre collocation method with 40 nodes. NLP tolerance 1e-6.

I've tried local discretizations as well (through PSOPT). They run much more quickly and consistently return reasonable results. There's still some noise in the control output though, but they're worth investigating a bit more.

Edited September 1, 2013 by alterbaron

[tavert]

That looks great! Can you share the source for your model?

It would be interesting to plot thrust divided by mass. In the challenge, we were seeing that applying an acceleration limit of 22 m/s^2 for the entire flight, even after starting the gravity turn, did give better efficiency - so your results do make sense. The challenge craft was using a LV-T30, so the max thrust should be 215 kN, not 200. Maybe with more thrust you'll make orbit with some mass left over.

[Stochasty]

Something I would like to see plotted during the gravity turn is an "effective terminal velocity" computed using gravity times the cosine of the thrust angle - the idea being that, as you begin the gravity turn, a smaller portion of your thrust is fighting gravity losses, and therefore the balance between gravity losses and drag losses should occur at a lower total velocity.

[alterbaron]

The challenge craft was using a LV-T30

Hm, looks like I used the wrong engine. The correct engine has more thrust and less mass.

I took another look at the challenge thread, and re-adjusted my parameters to properly match the challenge craft.

This time, it was able to orbit successfully, with fuel to spare!

Craft parameters:

Isp(1atm) = 320s

Isp(vac) = 370s

Tmax = 215kN

Fully-fueled mass: 11.15 t

Dry mass: 3.15 t

Mass after circularization at 74km: 3.1659 t

EDIT: Version with more collocation points added: http://imgur.com/a/NuHhi

Model source code: https://gist.github.com/alterbaron/6407325

Needs to be tidied up a little, but I hope it's clear enough to follow.

Edited September 1, 2013 by alterbaron

[Stochasty]

Hmmm... it's a much closer fit, but velocity is still well below even my "effective" terminal velocity through the gravity turn, despite having thrust to spare. That's very interesting.

Edit: Are you accounting for centripetal acceleration in your terminal velocity computation?

Edited September 2, 2013 by Stochasty

[alterbaron]

Edit: Are you accounting for centripetal acceleration in your terminal velocity computation?

The above plot doesn't include any centrifugal force due to the rotating frame of reference.

I played around a bit with including fictitious forces in the computation, and arrived at a model that fits the data pretty well up to the 30km mark.

In the force balance used to determine the terminal velocity, I included both the centrifugal force as well as the component of the coriolis force opposite to the direction of gravity.

Then, I applied your thrust angle idea to the new terminal velocity equation:

F_drag = sin(Thrust_angle) * (gravity + centrifugal_force + coriolis_force*cos(Velocity_angle))

This gave quite a nice fit:

Also, I took a look at the thrust angle relative to the orbital velocity direction (prograde indicator on navball):

Seems to suggest that above 20km (and with a craft with a high TWR), one should essentially start following the prograde indicator on the navball. This could be explained as maximizing the useful work done by the engine.

Edited September 3, 2013 by alterbaron

[NikkyD]

Sorry to pull up this old thread but the topic is rather timeless.

I read about half the thread and understood about half of that... my understanding of engl. math terms sucks

So could you help me / tell me if the following description of the problem is somehow accurate ?!

If i have a rocket with max thrust = gravity (TWR 1) it would not be able to lift off. For every bit of thrust above that i could make it to orbit, given infinite fuel (for simplicity lets assume the rocket has constant mass and infinite fuel).

If i go really slow, i dont waste fuel on fighting drag but i fight gravity a lot.

If i go really fast i fight a lot of drag but fight gravity only a little, cause it wont take long.

So if i keep a balance of 1 to 1 of the drag/gravity forces i fight them both equally so it should be best...

BUT because the drag has velocity squared in it, a little less velocity would greatly reduce the drag and only add a little more time = gravity.

I wrote myself a little program that simulates kerbal physics and i had it run vertical ascends to 80 km with throttle limited to x*terminal velocity. For one TWR combination 90% was the best solution (least fuel).

I googled a lot and found some explanation that in reality the best speed is a bit above TV and is limited by max_q...

Im kinda confused. Everyone says its the best and the formulas seems to make sense... but everyone is leaving something out to simplify. Does no one believe that the real solution might be somewhere above/below TV ?! Or is everyone convinced that TV is an ironclad solution ?!

[tavert]

Post the source of your simulation and we can look over it. Were you accounting for the reduction in effective weight due to the rotation of the surface of Kerbin?

[K^2]

The vertical ascent problem has a number of assumptions in it that might not be fulfilled.

1) Atmospheric density is either constant or drops exponentially. (Exponential in KSP, a little different in real world.)

2) Drag is proportional to mass. (True in KSP, but not in real world.)

3) Gravity is constant. (Not exactly true in either KSP or real world, but very close.)

4) There is no rotation. (Not true, but doesn't make much difference either.)

5) The target is "infinitely" removed. So there is no coasting to reach the necessary altitude. (Depends very much on how you run your tests.)

All of these affect a real launch, which is why real TWR for optimal ascent might not be exactly 2. This also doesn't take into account the gravity turn which makes things very different.

These are just some general points which may result in you getting something a bit different. As tavert suggested, if you post your code, we can probably say something more specific. But under above assumptions, optimization TWR=2 is exact. It's very simple to derive for constant atmospheric pressure, but takes a bit of work for exponential atmosphere of KSP. Result is exactly the same, though.

[NikkyD]

Well, basically what i do is

calculate drag and grav with altitude[now] and velocity[now].

Then i run some back'n'forth estimation to find the best thrust factor to achieve a certain drag ratio in the next calcuation.

with this factor i calc the acceleration as accel = thrust*factor -drag -grav

velocity[next] would then be velocity[now] + accel / timedelta

and altitude[next] is velocity[next] / timedelta

1) atmo is exponential KSP

3) grav is KSP

4) no rotation is considered

5) this is more tricky. I do at every step calculate the AP that could/would be reached and stop accelerating. This does not achieve the AP that is calculated at this point, because future drag/grav may lower it, but its a really small deviation. Coasting is then calculated like above until speed <= 1 and the result is "printed" out

I wrote the stuff in android/java because i work with that stuff so its easy for me and ill add visualisation stuff once it works better. Right now the calcuation time for each ascend is about a second (on my Galaxy S2 with 1,4 Ghz).

Edit:

funfact...

With my sim, if i run the test with a rocket with TWR below 2,1 it can never reach terminal velocity, its close but cant reach, because its always gaining just enough altitude to reduce the drag so that it matches its gain in velocity.

(constant mass of 10 tons, 210 kN thrust, starts coasting at 42148 m with speed 784,3 m/s to reach 80 km altitude)

Edited January 27, 2014 by NikkyD

[mhoram]

Well, basically what i do is

calculate drag and grav with altitude[now] and velocity[now].

Velocity should be the the velocity relative to the atmosphere - that was one mistake tavert made me aware of during my implementation.

Then i run some back'n'forth estimation to find the best thrust factor to achieve a certain drag ratio in the next calcuation.

There should be a way to replace the "back'n'forth estimation" by a one-way-back calculation beginning with target altitude and target drag ratio. This will increase the performance.

with this factor i calc the acceleration as accel = thrust*factor -drag -grav

velocity[next] would then be velocity[now] + accel / timedelta

and altitude[next] is velocity[next] / timedelta

You probably mean?!

accel = (thrustforce*factor -dragforce -gravforce) / mass

velocity[next] = velocity[now] + accel * timedelta

altitude[next] = altitude[now] + velocity[next] * timedelta

Edit: And these values are handled as vectors and not scalars.

Depending on the magnigude of timedelta the last two formulas can lead to significant rounding errors.

The larger the timedelta the larger the rounding errors.

The reasong for the issure of 90% terminal velocity is hard to guess without seeing the code.

Edited January 29, 2014 by mhoram

[fibonatic]

I agree with mhoram about your formulation of acceleration. And your integration method is called explicit Euler method, which is probably the most intuitive and easy numerical integration method, but also one of the least accurate with a global error proportional to your time step-size. There are a lot of different numerical integration methods for ordinary differential equations out there, a well known family are the Runge–Kutta methods, which can easily get a global error proportional to your time step-size to the power of 5. So by decreasing your time step by a factor 1/2 will increase your accuracy by a factor 2^5=32. And there are even methods which have dynamic step-sizes, which will reduce your number of calculations while in a relative "smooth" part of your simulation.

With my sim, if i run the test with a rocket with TWR below 2,1 it can never reach terminal velocity, its close but cant reach, because its always gaining just enough altitude to reduce the drag so that it matches its gain in velocity.

It is correct that you need a TWR higher than 2 to maintain TV, since the TV will increase with altitude, which means you need extra acceleration to achieve the new TV. (gravity will also decrease but a lot slower than the increase of TV).

PS: I have tried to tackle the 2D optimal ascent trajectory my self using the toolkit ACADO[\url] within MATLAB, since I also would like to know what good ascent profiles are on other atmosphere bearing celestial bodies. These are the differential equations:

where v_r is the radial velocity, r the radius, v_theta the tangential velocity, omega_s the angular velocity of the surface/atmosphere, mu the gravitational parameter, alpha_r the applied acceleration (force per mass) in the radial direction, alpha_theta the applied acceleration the tangential direction, H the scale height of the atmosphere and C a composed drag parameter (which removed the need of exp(r0/H)).

However every attempt of a script seems the have the problem that it ignores all/most constrains, such as final radius and velocity or just jumps to it on the final step.

[tavert]

I tried ACADO (via C++ https://gist.github.com/tavert/7783938) with KSP physics a few months ago, it couldn't even solve a strictly vertical ascent problem. The authors are also absolutely unresponsive to most issues that get reported on their discussion board.

Not sure trying to eliminate time and solve with radius as the differential variable is a great idea. You'll need to take into account the variation of Isp vs altitude if you want to minimize actual fuel use, and there can be trajectories for which altitude versus time is not an invertible function.

[K^2]

I'm really starting to think that the "simplest" way to solve this problem is a genetic algorithm with simulated annealing. I guess, I'll play around with that.

[NikkyD]

Velocity should be the the velocity relative to the atmosphere

I dont understand.

accel = (thrustforce*factor -dragforce -gravforce) / mass

Makes no difference ?!

And these values are handled as vectors and not scalars.

Please elaborate as well

[mhoram]

About the issue of vectors, I meant, that you have to represent the velocity of the ship as the horizontal and a vertical component.

The forces also have horizontal and vertical components, so

accelX = (thrustforceX*factor -dragforceX -gravforceX) / mass

accelY = (thrustforceY*factor -dragforceY -gravforceY) / mass

and the directions of the three forces are different:

- Gravity points towarts the planets core

- Thrust points towarts the current orientation of the ship

- Drag points in the opposite direction of the current momement of the ship

A simple subtraction of the absolute force values will lead to errors.

If you implemented the status of a rocket with the compoments altitude, horizontal velocity and vertical velocity then the issue of the atmosphere is irrelevant, because you most likely assume already that the ship moves relative to the atmosphere.

[NikkyD]

Its a pure vertical ascent, there is no horizontal component at this time (simplified).

[K^2]

I've been playing around with all of this, and I'm having trouble getting simulation to match the game to a satisfactory precision.

For a ship, I have an FL-T100 with LV-909 underneath and Mk-1 on top. Oh, and Mk-16 parachute for not killing the test pilot. That gives me dry mass of 1.463T and 0.4995T of fuel on the pad.

The engine produces constant 50kN of thrust. I_SP is computed according to atmospheric pressure, and I'm using g0 value of 9.82m/s², same as the game's code, to get the fuel flow.

Speaking of the atmosphere, Kerbin has legacy atmo, meaning pressure as function of altitude is P = exp(-h/5000m) atmospheres. I multiply that by 101325.0/82843.125 kg/m³/atm to get density at altitude h.

This leads me on to drag, which is computed as 0.5*0.008*ÃÂ*d*v². The value for d is computed dynamically as fuel is used up, but it's so close to 0.2 that it might as well be that.

Finally, I add in acceleration due to gravity g = μ/(h+R)² with μ = 3.5316000x10¹² and R = 600km. As well as centrifugal force, ɲ(R+h). Both of these are consistent with how the game computes them. I do ignore Coriolis force, but since the craft lands almost next to where it took off from, it cannot be significant.

I've tried integrating with Euler and Verlet using a 20ms time step, which seems to be default physics time step for Unity, and I think that's what the game uses at Warp x1. The problem is that in my simulation, the craft overshoots by almost 10%. In KSP, I get maximum altitude with throttle to the wall at 5488m, while simulation yields 5840m.

In addition, I've tried comparing other parameters. Fuel cutoff seems to happen at the correct time, so I'm pretty sure fuel usage is fine. However, right before cutoff, my simulated craft is going at 216m/s, and the game tops out at about 206-207m/s. Similarly, if I just drop the ship in my simulation from 5488m/s, the impact happens at about 106m/s, while in KSP I get something clsoer to 104m/s.

So error in drag computation would be my first guess. But I just can't find any problem with it. The equation I'm using for drag computation is consistent with both the wiki and what I've been able to dig up in the game's code. Of course, I don't expect perfect agreement, nor do I need it, but 10% is much too much for something that seems to be using all the same computations. Any suggestions?

[mhoram]

For a ship, I have an FL-T100 with LV-909 underneath and Mk-1 on top. Oh, and Mk-16 parachute for not killing the test pilot. That gives me dry mass of 1.463T and 0.4995T of fuel on the pad.

I come up with 1.4625T dry mass and 0.5T of fuel for these parts. (I assumed Command Pod Mk1)

This leads me on to drag, which is computed as 0.5*0.008*ÃÂ*d*v². The value for d is computed dynamically as fuel is used up, but it's so close to 0.2 that it might as well be that.

Somehow I am missing the ships mass in this formula. (Area = 0.008 * currentMass)

Finally, I add in acceleration due to gravity g = μ/(h+R)² with μ = 3.5316000x10¹² and R = 600km. As well as centrifugal force, ɲ(R+h). Both of these are consistent with how the game computes them. I do ignore Coriolis force, but since the craft lands almost next to where it took off from, it cannot be significant.

What starting altitude do you use? I usually start the ascent around 70-80m

[tavert]

Did you remove the monopropellant from the pod? And what do you mean by "legacy atmo"?

[K^2]

500g definitely won't make a difference, but I'll recheck my numbers. For drag, yeah, all of that times mass. Since everything but thrust is proportional to mass, I've been working with acceleration, but I should have clarified. And I'm starting at 75m.

tavert, game supports two atmospheric models. One is exponentially decaying, referred to as legacy atmosphere, and the other is based on curves. Later can be used to model a more realistic atmosphere. There are some mods out there that already support it, but vanilla planets use the exponential atmosphere.

So I started putting together a genetic optimizer. It's still a bit rough, and might need better annealing, but it looks like a good start. This is 1D, using same simulation as above. (So this isn't a perfect replica of KSP, but qualitatively it should match, at least.) The only optimization parameter is maximum altitude achieved.

Altitude (left) and velocity (right) as function of time. Full throttle ascent in red, optimized in blue.

The thin blue line on the velocity graph indicates terminal velocity at altitude achieved by optimized flight.

And this is the throttle position (left) and resulting TWR (right). Only first 45 seconds shown, because engine cuts out after that.

Yes, it's noisy. Yes, I can do something about the noise, but it's not a priority right now. It's good enough for a test. It does show an interesting feature, though. As expected, the ship takes off with full throttle, and then dials it down as ascent settles on terminal velocity. TWR holds around 2 at that point. However, right before the engine cutoff, throttle goes up again. Apparently, to give the ship that one last push before it begins coasting. That bit was somewhat unexpected.

So, I guess, now I just have to figure out a good way to do this in 2D and give it enough fuel to make orbit.