Finding the best ascent path for rockets

[mhoram]

In the Payload Fraction Challenge Thread came up a discussion about finding the optimal ascent strategy for rockets.

Well here is my take on this problem.

A perl program simulating the physics of KSP (implementation was based on V 0.21.1).

At the moment it is missing documentation, but it consists of different modules like Rocket & Stage design, Physic constants, gravity turn handling, Kepler and Cartesian coordinates.

The basic way it works is that the ascent profile is parameterized by

- pirchover altitude

- pitch angle

- pitchover duration

- twr during downrange

- target orbit radius

The program does not find the best ascent profiles yet, as some details are still missing in the implementation, like elevated pitch during the LV-N-phase.

Other details are simplified like orientation can change instantly without delay. No time-delay at the end of a stage.

Using this simulator with multiple sets of random values the results are compared by

1. Is the orbit reached?

2a. If yes, remaining fuel

2b. If no, highest Periapsis

The best set of paramerers is kept and partially randomized for the next set of random values to be tested.

This procedure is iterated until one is content with the results.

You can find the code (released under GPL V3) here:

https://github.com/mhoram-kerbin/ascent-optimizer

Edited November 17, 2013 by mhoram

[Vanamonde]

(Thread moved to gameplay questions, where it's more likely to receive helpful attention.)

[Kasuha]

I don't know perl so the code is not very readable to me but it seems to me you assume some predefined ascent trajectories and in general you work with very few parameters. If that's true (I may be wrong) it's scanning just a very small subspace of all possible ascents and while it may be near the optimum there's rather high chance the optimum will not lie in the subspace you scan. For instance I am very uncertain whether offcenter thrust will be only limited to the LV-N phase. Also there's no point to make LV-Ns any special, it's always possible certain stage will be underpowered or overpowered so at least each stage could likely use its own approach. Or there may be need for previous stage to prepare conditions for the next stage so the granularity could be even higher.

I also wonder which ascent trajectory you assume for non-equatorial liftoff point. It would probably to make sense to add liftoff azimuzh or inclination as one of initial parameters (unless it's already there, I can't read the code well).

[Seanner]

I'm working on the exact same thing. Currently in the testing phase for 1d motion to confirm that my algorithm comes up with terminal velocity ascent.

[numerobis]

My KSP-scripts have this ascent.py module that is relevant... I haven't touched it in a year.

Perhaps more important than that is the engine.py and planet.py modules. I use those all the time when designing spacecraft.

Edited November 18, 2013 by numerobis

[tavert]

I recommend you folks read through this thread from a few months back: http://forum.kerbalspaceprogram.com/threads/46194-I-need-someone-help-me-do-some-math-for-launch-optimization?highlight=launch+optimization

3 main approaches for this optimization problem: pseudospectral methods, local discretization, or using the Pontryagin principle to express the solution in terms of a boundary value problem. Alterbaron got a pseudospectral method using PSOPT (and possibly local discretization via an option setting with the same software?) working for the single-stage problem, fixed design. Multiple stages are trickier, you have to split up the time horizon for each stage, scale the time variable by a different interval in each stage, and stitch together the dynamic conditions between the end of one stage and the start of the next.

[mhoram]

@Kasuha you are right, that I only look at a small fraction of all possible ascents.

In the current sourcecode the same ship as in http://forum.kerbalspaceprogram.com/threads/39196-Launch-Efficiency-Exercise-Updated-for-0-21-1 is tested for it's optimal ascent and I get a remaining fuel of about 8.1 while the best in the challenge has 10.55 left. So yes there is improvement possible. To test the result, kOS seems to be an interesting addon.

Offcenter thrust is probably necessary in the downrange phase whenever the TWR is low.

Non-equatorial liftoff is an issue that I did not adress yet.

@Seanner I found most information about KSP physics in the Wiki and some forum posts,

For general physics I found https://downloads.rene-schwarz.com/dc/category/19 and https://downloads.rene-schwarz.com/dc/category/20 very helpful.

The only thing that seems to be unknown is the handling of rotation by torque and gimbaling.

@tavert Thanks for the link. The thread is a treasure trove friom before my KSP-time! PSOPT looks very promising. I think I will have sometime soon a deeper look at it.

Edited November 18, 2013 by mhoram

[numerobis]

the handling of rotation by torque and gimbaling.

The gimbals max out at 1 degree. You can largely ignore the effect of off-center thrust in practice, it's less than a part in a thousand. How does command-pod torque matter to the question at hand?

[Seanner]

Thanks I've got most of it now, but I'm ignoring centripetal and Coriolis for the time being.

By the way, I wasted a lot of time figuring out why my algorithm got progressively worse over increasing iterations and generally performing horribly until I realized that I was returning the minimum instead of maximum performer at the end. (WHY DOES MAX HEIGHT ARE ZERO?!?) True story...

edit: Found another problem...incorrect fuel calculation was causing excessive fuel burn. I was using altitude-dependent real g instead of 9.81...survival of the fittest now produces nearly perfect throttle genes. Now that this garbage works, I'm going to try to evolve genes for a neural net to fly it instead of directly encoding throttle. Hopefully this results in a brain that can fly any craft instead of the specifically-solved-for one

Edited November 19, 2013 by Seanner

[mhoram]

PSOPT is a bit tricky to install because of it's dependencies. Will report back, when I got it running.

The gimbals max out at 1 degree. You can largely ignore the effect of off-center thrust in practice, it's less than a part in a thousand. How does command-pod torque matter to the question at hand?

I think there are two different applications of off-center thrust in this discussion: thrust caused by gimbaling and thrust while the orientation is different from prograde. While the effects of gimbaling should be quite low, the other must be taken into account, especially since the angle betweeen orientation and prograde can be large.

Torque is the only physical property, that I was unable to model correctly because there is no information about it anywhere to be found.

[numerobis]

Seanner: g0 is exactly 9.82 in KSP (but gravity at sea level is about 9.81).

[Camacha]

Now that this garbage works, I'm going to try to evolve genes for a neural net to fly it instead of directly encoding throttle. Hopefully this results in a brain that can fly any craft instead of the specifically-solved-for one

I try to code my scripts to be as universal as possible, usings sensor where possible and assuming nothing that does not have to be assumed. But what you are doing is in another league.

<Appropiate pic>

[Seanner]

But why would they make it 9.82?! This throws off everything. I don't even know which way is up. Launch angle = -pi/2 from now on.

[numerobis]

As long as you mount the rockets on the opposite side of the rocket, that launch angle works fine.

[mhoram]

Finally I got PSOPT running thanks to alterbaron.

First I will have a look at alterbaron's formulation and solution - and have a look at how it might be extended: the multistage-problem is very appealing to me.

[tavert]

Cool, did you take notes about what you had to tweak to get it to run? I never actually ran alterbaron's PSOPT code, but I have tried implementing something similar off-and-on in various other tools (ACADO, Pyomo, Mathematica, custom OSiL) and nothing has quite worked yet.

[mhoram]

Sorry no, but feel free to ask if you run into any problems. There were some manual adjustments needed.

[mhoram]

This project kind of resulted in another ...

Since I need the KSP physics as part of the PSOPT source code, I have to implement the physics of KSP ... again in C/C++.

It is time to gather all physics-related information in a single document.

https://www.dropbox.com/sh/jobzpt7feosc3rt/b4PwI_5WXi

And here is the repository for the documentation: https://github.com/mhoram-kerbin/ksp-physics-documentation

Be aware ... it is still work in progress.

[mhoram]

I managed to implement my first testrocket in PSOPT.

In contrast to alterbarons implementation I use a full 3d coordinate representation.

The target was to optimize the thrust in such a way that maximizes the rocket's mass at the point in time where the orbit's periapsis reaches 75km.

Edit: These results are based on a wrong implementation of the KSP physics and should be ignored

Javascript is disabled. View full album

Launches at 0°5'49" South like the Kerbal Space Port and has two stages

According to the theoretical results, it has a final mass of 1880 kg.

Mechjeb (10, 70, 0 40) performed worse with a final mass of 1814 kg. (these are my default values and I did not try to optimize them for this rocket)

One thing I find strange and can not explain is why the gravity turn starts directly after liftoff. Orientation-pitch is near 80° at the time of liftoff.

Source: https://github.com/mhoram-kerbin/ksp-launch-optimizer/tree/38d382fba55b5e54a913c217f9988f9998d2caf1

Edited January 8, 2014 by mhoram
wrong results

[Kasuha]

Interesting.

Regarding the pitch, probably the benefit of extra horizontal speed exceeds the loss of longer travel in thicker atmosphere.

Another thing that came to my mind was that upper atmosphere is moving at higher speed than surface. But the difference is not that big - at 10 km it's just about 3 m/s.

It might be interesting to see the angular difference between motion vector and thrust vector over the ascent.

The thing that I find strange is that it is increasing thrust right before the staging.

Edited December 31, 2013 by Kasuha

[tavert]

Well, for one thing your drag force is being applied in the direction opposite to orbital velocity instead of surface velocity.

[mhoram]

Well, for one thing your drag force is being applied in the direction opposite to orbital velocity instead of surface velocity.

You are so right. Thanks for the review! Now the pitch seems much more reasonable.

Final mass: 1856 kg

I should add that for the implementation I assume that above 70km there is still an atmosphere that causes drag. Using a density of 0 caused inconsistencies in the calculations. Perhaps I will multiply it with a modified arctan(altitude) function - tests will follow.

Javascript is disabled. View full album

It might be interesting to see the angular difference between motion vector and thrust vector over the ascent.

The thing that I find strange is that it is increasing thrust right before the staging.

I added pitch of motion (ground vel and orbit vel) to the graphs - and also the thrust.

And the increase in thrust before staging is quite simple: at that altitude that atmospheric drag decreases and allows a faster velocity, so the thrust increases. At this altitude the terminal velocity increases faster than the rocket can increase it's velocity.

source: https://github.com/mhoram-kerbin/ksp-launch-optimizer/tree/a97d823df7b5ae7ae1bd8c9de0ecf79ed992ef97

Edited January 8, 2014 by mhoram

[tavert]

You are so right. Thanks for the review! Now the pitch seems much more reasonable.

Thanks for posting your source. Linus' law, given enough eyeballs all bugs are shallow. Looks like you're getting closer to something usable and useful here. Just needs some setup guidance for anyone other than you or alterbaron to get PSOPT working from a clean slate. Have you ever used Travis CI? Could be worth setting up here, gives you free testing and helps you establish a blueprint for how to set up the dependencies in a controlled environment.

It's a pretty minor difference, but I recommend using the gravitational parameter (aka mu or G*M) instead of G and M individually, since I believe planet mass is a derived quantity in the game. The product is the value that really matters.

Makes no difference for Kerbin, but the wiki's a bit confusing on atmosphere height - you should use the version that does not depend on sea-level pressure. The second equation on the wiki that depends on sea-level pressure gives the wrong values. On Eve for example, the atmosphere ends just under 97 km, or ln(10^6) = 13.8 times the scale height.

For handling the atmosphere cutoff you could either try a sigmoidal function (arctan, or tanh, or erf, or a logistic function - watch out for overflow though) to smooth out the small discontinuity in pressure, or make 2 assumptions: that the last staging event happens before the craft leaves the atmosphere, and the craft only leaves the atmosphere once (doesn't reach apoapsis then start falling back into the atmosphere before it gets into orbit, might happen for low TWR). With those assumptions, you could add another problem phase for the final portion of the trajectory that happens in vacuum.

And the increase in thrust before staging is quite simple: at that altitude that atmospheric drag decreases and allows a faster velocity, so the thrust increases. At this altitude the terminal velocity increases faster than the rocket can increase it's velocity.

It's started to pitch over by that point so the story's a little more complicated than that I suspect.

Edited January 1, 2014 by tavert

[mhoram]

Just needs some setup guidance for anyone other than you or alterbaron to get PSOPT working from a clean slate. Have you ever used Travis CI? Could be worth setting up here, gives you free testing and helps you establish a blueprint for how to set up the dependencies in a controlled environment.

Have not yet used Travid CI. Looks interesting - but with that much theory during the last few days, I would like to get back playing a little bit ;-)

The second equation on the wiki that depends on sea-level pressure gives the wrong values.

You are right that the wiki contradicts itself on the same page.

For handling the atmosphere cutoff you could either try a sigmoidal function (arctan, or tanh, or erf, or a logistic function - watch out for overflow though) to smooth out the small discontinuity in pressure, or make 2 assumptions: that the last staging event happens before the craft leaves the atmosphere, and the craft only leaves the atmosphere once (doesn't reach apoapsis then start falling back into the atmosphere before it gets into orbit, might happen for low TWR). With those assumptions, you could add another problem phase for the final portion of the trajectory that happens in vacuum.

First I will test the sigmodial functions, because setting up an additional phase in PSOPT is more complicated.

Edited January 2, 2014 by mhoram

[mhoram]

First I will test the sigmodial functions, because setting up an additional phase in PSOPT is more complicated.

As it seems PSOPT has problems finding a solution when the pressure is multiplied with a logistic function.

Multiplying it with (1 / (1 + exp(f * (altitude - AtmosphericHeight)))) makes the solving-algorithm not converge for f > 0.001.

So I will stick with the unmodified (wrong) pressure function, since this modified version causes a considerable amount of approximation error for f = 0.001.

As an input for the calculation one needs mass, propellant mass, ISP_1ATM and ISP_VAC for each stage. Since there is no easy way to get these infos from an existing rocket (MechJeb and Engineer redux lack this functionality for getting ISP), I am currently writing a plugin in the very early stages to get these infos directly from the VAB. It is available under GPLV3 at https://github.com/mhoram-kerbin/psopt-integrator.