Mini-challenge: max altitude with this supplied spacecraft

[Kosmo-not]

Ok, now for explaining how the throttle is adjusted:

I have one column that calculates the velocity needed for drag to equal weight. A second column calculates the acceleration needed to get to that speed within the next time interval (it also adds on acceleration due to the change in needed velocity from the previous time interval). The third column calculates the thrust required to achieve the needed acceleration. Finally, the last column sets the rocket\'s thrust (within bounds of max thrust and zero thrust).

Formulas:

V_needed = SQRT(2*F_g/(?*m*C_D))

a_needed = ((V_needed - V) + (V_needed - V_{needed,previous}))/?t

thrust_needed = W + F_D + a_needed*?t

To further improve this, a predictor equation for the change in needed velocity will have to be employed.

@ PakledHostage

Can you provide some data to prove that C_d is not used?

To support that it is used, I\'ll refer you to my drop test with parachute deployed.

http://kerbalspaceprogram.com/forum/index.php?topic=5623.msg92588#msg92588

[closette]

Seconded - the maximum_drag factors (let\'s not call them c_D\'s, which are unitless) are definitely used. You can edit the values in the .cfg file for, say, a command pod, and then launch and drop it from a large height. If you make the maximum_drag factor over 100, you won\'t need a parachute for it!

Again, in KSP, the drag force on any given part F_D = 1/2 x maximum_drag x mass x density x v². Add up the force on separate parts (no matter how they are arranged) to get the total drag force.

[closette]

@Kosmo-not,

That\'s a very successful 'shooting' method you are employing, and your first solution is pretty close to optimal. I bet that if you iterated for one more pass, using (TWR)best = 2 +2a/g, you\'d get even closer.

For the 'target' velocity if you stick with TWR=2 on your first pass, as you did, that is already a known function of height (see my PDF above) since it is the same as the local terminal velocity. You could try a first pass with target velocity =v_T=97.8 m/s exp(+height/9891m). (*)

One clever part of your scheme is that it can switch from full throttle at launch to the correct amount once the best speed has been achieved. I was thinking about using piecewise functions to do that but your method is much more flexible.

(*) The sea-level terminal speed is 97.8 m/s is for an object with mass-averaged maximum_drag =0.2. It is probably a bit less for the ship in this challenge, but only depends on the square root of maximum_drag so will be close.

[Kosmo-not]

Ok, I made the change as you suggested. It increases the maximum altitude by a small amount. I\'ll also have to make some changes because the drag coefficient for the rocket varies with time, thanks to that tri-coupler.

[closette]

The drag coefficient for the rocket varies with time, thanks to that tri-coupler.

OK you\'re teasing me, please explain! The multiplier of (density x v²) should be constant for that part, right?

I am pleased that my suggested changes made a small difference in the positive direction.

We really need someone who can do optimal control theory to find the theoretical max. height for this craft. Quite a few people here (and whom I know) tell me it\'s a straightforward problem, but no-one has come out and just solved it!

[PakledHostage]

Can you provide some data to prove that C_d is not used?

To support that it is used, I\'ll refer you to my drop test with parachute deployed.

http://kerbalspaceprogram.com/forum/index.php?topic=5623.msg92588#msg92588

OK. You got me there! I admit I didn\'t think of that obvious example... I only noticed Closette\'s description of how to calculate C_d for an assembly:

'For an assembled spacecraft replace the maximum_drag with a mass-weighted average of all the parts\' maximum_drag factors, which usually results in a value of 0.2-0.3.'

after I made my last post, but I don\'t know how that method was arrived at. Maybe one of you two can fill the rest of us in?

My (apperantly flawed) idea that C_d is not used in the game is based on my own numerical integration of the two 'drop tests' that I described earlier in this thread. As I mentioned above, my model (using the two equations that I posted above, and which are the same equations that you describe) allowed me to make very accurate predictions of the two flight trajectories that I tested.

As unlikely as it would have to be, it may be that the 'mass-weighted average' maximum_drag factor for the two (very different) stacks that I tested was 0.2027. This would have made it appear that the C_d term was not used because C_d would have been the same for both my tests as it was for your test.

Closette\'s post above claims that the 'mass-weighted average' maximum_drag factor typically ranges between 0.2 and 0.3, so this is a possible explanation.

I\'ll trust that that is the case. But please, for the rest of our benefit, can you please explain where using the 'mass-weighted average of all the parts\' maximum_drag factors' comes from?

Thanks.

PH.

[Kosmo-not]

@ closette

drag coefficient at launch: (12.8*0.2 + 0.8*0.3)/13.6 = 0.20588

drag coefficient at 0 fuel: ((12.8-4.4)*0.2 + 0.8*0.3)/(13.6 - 4.4) = 0.20870

It\'s not much of a difference, though.

@ PakledHostage

It comes from summing the drag forces. It\'s easier to work with an equivalent drag coefficient for the mass of the entire rocket.

[PakledHostage]

But please, for the rest of our benefit, can you please explain where using the 'mass-weighted average of all the parts\' maximum_drag factors' comes from?

Never mind. I just figured it out... I\'ll post the explanation if anyone wants it, but I don\'t have time to make an image of the equation right now.

PH.

Edit: I\'m adding the aerodynamic drag and density equations discussed in this thread to this post for the benefit of future readers who, like me, require things to be a bit more obvious:

In-game aerodynamic drag equation:

where:

m is the instantaneous total spacecraft mass

C_d is the 'mass-weighted average of all the parts\' maximum_drag factors'

The 'mass-weighted average of all the parts\' maximum_drag factors' is defined as:

where:

n is the number of individual components in the stack

m_j are the individual masses of each of the spacecraft components (from the VAB)

Cd_j are the individual drag coefficients of the spacecraft components (from the VAB)

Kerbin\'s atmospheric density (rho in the drag equation above) can be approximated as:

where:

Altitude is in metres.

And while solving the Goddard problem is FAR beyond my own mathematical abilities, I have benefited from Kosmo-not, Closette, et al\'s efforts in this and related threads. I hope to use these equations to more accurately predict my re-entry trajectories and for planning more precise aerobraking manoeuvres.

Thank you all for your excellent work!

[closette]

OK I just figured it out as you were replying - the tricoupler is the only part in there with maximum_drag not equal to 0.2. Sorry to take up your time.

(As a favor please don\'t call maximum_drag a 'drag coefficient' as this will confuse other readers who will try to apply the terrestrial, real-world drag equation and we\'ll have to \'splain it all over again! No need to reply). I know what you mean though!

If you can tell the difference between mass-averaged maximum_drag of 0.20588 and 0.20870 you must be very, very close to optimum!

The mass-weighted average maximum_drag for an assembly of parts 1, 2, ...

= (maximum_drag₁ x m₁ + maximum_drag₂ x m₂ + ...) / (m₁ + m₂ + ...)

The numerator x density x v² gives the total drag force. The denominator is just the total mass of the craft.

That way the total drag force on the craft can be written as:

F_D =mass-weighted average maximum_drag x total mass x density x v^2.

[Kosmo-not]

I fixed an error in my spreadsheet. Now it\'s reporting a maximum altitude of about 33,500. Tomorrow, I\'ll add in a column for a changing maximum_drag and an accurate fuel cutoff. I don\'t think it will be possible for anyone to get to 35,000.

[jebbe]

About the theoretical solution: I don\'t think it\'s as simple as it seems to be. While it is fairly obvious that there is a solution, to actually obtain it one has to solve an inverse problem. However, since the problem is non-linear and since we have some mathematically strange constraints (mainly from burnout: there\'s a minimum mass, and the end time isn\'t fixed), this is not straight forward. There\'s a lot of theory about that stuff, in this case one would probably go with an adjoint method, but I don\'t know enough about them to do it for this problem. For problems with a small parameter space one can usually use Monte Carlo (try a lot of random thrust profiles) or brute force (try all possible thrust profiles), but my attempts in that direction have failed so far: Monte Carlo doesn\'t converge, and brute force doesn\'t get to a sufficient resolution in thrust and time. So it seems that there\'s no way around this kind of stuff: http://cas.ensmp.fr/~petit/papers/ifac2008/0637.pdf. Probably one could email the authors, obtain their MATLAB-script and adapt it to kerbal parameters? Their problem is very close to what we want to solve...

@closette,

thanks for the pdf, makes sense now. So the way I understand it: a TWR of 2 is optimal if we are close to the surface (g=const), and if we don\'t accelerate much, which is the case after the initial full throttle ascent and before we reach higher altitudes.

[jebbe]

I just came up with another idea how one could find a theoretical optimum. I\'m still missing one essential idea, but maybe somebody here can help?

The problem we want to solve is something like this:

dt(h) = v

dt(v) = u/m - D(h, v) - b/h^2

dt(m) = -u/c

with the time derivative dt(), distance from planetary center h, velocity v, thrust u, drag D, gravitational parameter b, and specific momentum of the engine c.

The idea is to discretize this by replacing dt(x) by (x_(i+1) - x_i) / dx for some suffiently large time interval [0..n] and set up a linear system of equations Ax=b with the discretized quantities x=(h_1..h_n, v_1..v_n, m_1..m_n, u_1..u_n), a matrix of coefficients A, and some right hand side b. Once we have that system, we could just solve it for x.

Obviously, there\'s a couple of problems with this:

[list type=decimal]

[li]The system of ODEs is nonlinear[/li]

[li]Thrust has a max and a min value[/li]

[li]We need to constrain the integral over the thrust to a certain value (burnout after no more fuel is left)[/li]

[li]The matrix A is not square, so the system has infinitely many solutions[/li]

[li]We need to introduce an equation that optimizes the height[/li]

The first three problems I could deal with:

[list type=decimal]

[li]Nonlinearity: Implement some iterative scheme to approach the solution. The whole thing is highly non-linear, so I don\'t know if it would converge, but with a clever choice for the linear problem I think it should be possible[/li]

[li]Thrust is bounded: Solve for some variable psi(u) that is not bounded and define a map between psi and u like eg tanh[/li]

[li]We have limited fuel: Just add one equation that sums up the thrust at all time steps and set it to a fixed value[/li]

However, the last two issues I\'m not yet sure how they could be solved, but they go hand in hand: I think one could define some penalty functions that pick and optimize the highest altidude achieved, but I\'m not sure how... If someone can come up with an idea I would be willing to go ahead and implement the whole thing.

As of now this is just an idea, and I don\'t know if it would work, but who knows?

[closette]

@jebbe,

I admit I\'m at my limit of math skills for this problem, but it\'s a challenge in itself. If you can go for a solution (sounds like you are close) I\'ll be cheerleading all the way!

If it makes things simpler, you can assume that 'g' is approximately constant at 9.7 m/s/s for this problem and just redefine your 'h' as altitude above Kerbin\'s surface. I know, 'g' does decrease by about 10% at 30000m but for a first attempt at these coupled equations you can probably ignore this.

In the paper you quoted, the bounding function psi(thrust) caused the authors problems where they had to dream up some penalty function later on. Ugh. Again only if it helps, and you are going with iteration, I suggest that thrust could be limited to (say) 6 possible values: 0, 120, 240, 360, 480, or 600 N at each time step. That\'s not much worse than the resolution one can get with the SHIFT/CTRL keys, and you could change the number later.

[jebbe]

@closette:

Sorry, but I think I have to disappoint you, this exceeds my math skills as well. I was hopeful with the Monte Carlo approach (where I used discretized thrust levels), but it seems that the problem is too complex for that. As I said, the standard approach here would probably be an adjoint method, but I\'ve never worked with them myself and unfortunately I don\'t have the time to spend weeks getting into that.

As of the other idea: So far this is only speculation, even if someone comes up with an idea for an optimizing condition I don\'t know if it would actually work. In any case it\'s far from standard, it\'s just an attempt to apply a method I know on a problem that is probably not really well suited for it.

So let\'s hope someone with actual knowledge in control theory picks this up!

[closette]

I\'m taking off for the weekend but will email ask one of the profs I know to send an email the authors of that paper - who knows they may want to help... hopefully their code is readable (sensible variable names and comments) and can be adapted to different spacecraft easily.

[Kosmo-not]

I cleaned up my spreadsheet and added the features I said I would.

The change in the maximum_drag factor didn\'t make much of a difference to the maximum altitude.

Attached is the spreadsheet if you wish to mess around with it.

The 'a factor' at the top is \'x\' in 'x*(a/g)' for determining needed acceleration.

I think it\'s accurate to within a few percentage points.

[closette]

What is the maximum theoretical altitude? I see different numbers in each Sheet of the Excel file.

[Kosmo-not]

My bad. I forgot about the extra sheets I no longer use. I posted the correct excel spreadsheet above.

[boolybooly]

OK two questions, first if KSP applies drag using a stepwise model related to altitude, shouldn\'t we apply thrust to attain max velocities in a stepwise fashion ?

Also, there is a big elephant in the room regarding engine efficiency, is it constant or does it vary at different thrust levels?

My best 33,319

[Causeless]

My brain has exploded.

Too much maths. ???

[jebbe]

OK two questions, first if KSP applies drag using a stepwise model related to altitude, shouldn\'t we apply thrust to attain max velocities in a stepwise fashion ?

Are you thinking about a layered atmosphere, as the height indicator with the different shades of blue seems to suggest? Some people here have figured out that there are no steps, the density of the air decreases exponentially, just like on earth. See e.g.:

density = 0.0098*exp(-altitude/5000)

it\'s not perfect, but it\'s pretty close

Also, there is a big elephant in the room regarding engine efficiency, is it constant or does it vary at different thrust levels?

What do you mean by engine efficiency? Are you talking about the delta-v you get out of a certain amount of fuel? That\'s the specific momentum, and as far as I know it\'s constant, regardless of the throttle setting. So in outer space it doesn\'t matter how vigorously you burn, if you burn a certain amount of fuel, you accelerate by one given velocity.

The efficiency of your burn as a whole on the other hand does depend on the throttle, but only because of gravity and air resistance. You want to escape gravity as fast as you can to avoid hovering (which consumes fuel without bringing you forward), but at the same time want to go as slow as possible to minimize air resistance. Optimizing between these two contradicting goals is known as the Goddard problem, which we\'re discussing here.

[closette]

The problem here among us is 'not quite enough maths'.

@boolybooly - thanks for trying the challenge, 33319 m is a worthy attempt. I note that not many of those who downloaded the spacecraft have actually posted their results.

As for stepwise application of the physics (drag model), I presume it is done with a time step that is 'small enough' to look continuous, otherwise we\'d see our ships rising up jerkily against the background. Kosmo-not\'s spreadsheet approach also uses steps (rows) to move the spacecraft a bit, subtract some fuel, calculate the new weight and drag, change the thrust etc. And the more steps used (subject to rounding errors) the closer we get to the continuous case.

As far as we know, the engine thrust at a given setting is constant with altitude and time, which is one simplification this problem has compared to real-world cases.

[boolybooly]

Well that\'s two illusions shattered in one afternoon!

I think it was because people like HarvesteR talked about the atmosphere cutting out at a certain altitude that I assumed since it had an outer boundary it also had some kind of internal layering to account for the subjectively sudden changes in density.

I had also laboured under the delusion that burning at low throttle ie <20% might be more efficient than burning at higher throttles, again subjectively. I always meant to test it and never got round to it until now, but I just have and you are right, it doesn\'t make a difference. Would prefer it did actually though that does make life simpler.

Learnt something useful anyway.

[closette]

@boolybooly,

Glad you learned something new - one of the reasons I submitted this challenge was so we could learn from each other, and I think this has been a success despite the low participation.

When someone tells us '100% throttle all the way, all the time', now we can point them to this thread.

I had an epiphany as to why Kosmo-not\'s spreadsheet scheme works so well, and I don\'t think we need fancy optimal control software.

Stay with me here!

I. Ascent through a uniform 'slab' atmosphere

Imagine a uniform-density atmosphere in a slab over Kerbin\'s surface. From the equation of motion (refer to the PDF attached to one of my posts above: http://kerbalspaceprogram.com/forum/index.php?topic=7161.msg102364#msg102364) we found that the most fuel-efficient speed up through this slab is equal to the craft\'s terminal speed, which in KSP\'s drag model would be ~ constant. Also TWR=2. So the best strategy to improve your fuel efficiency (in miles per gallon or mpg) would be (1) accelerate as fast as you can up to the terminal speed, then (2) stay at this speed, which would mean that your thrust/weight ratio = 2.

The initial maximum acceleration helps you in two ways:

- you spend less time at low mpg (fuel efficiency)

- at maximum thrust you reduce the craft\'s weight the fastest, leaving less to lift later on and greater fuel efficiency for the rest of the ascent. (As shown in my PDF, the instantaneous efficiency e is inversely proportional to weight W).

- but be sure not go above the terminal speed at any time, because then you\'ll be wasting fuel by fighting the increased air drag,

Still with me?

II. Adding a second, less dense slab

Now put another, less dense slab above the first one. The terminal speed in this slab is higher, so on crossing the boundary what should you do to maximize your mpg? As before, speed up to the new terminal speed as fast as possible, then throttle back to TWR=2 and stay ascending at that constant speed.

III. Putting it all together: ascent through an exponential atmosphere

Adding additional constant-density slabs on top would force you into the same strategy of 'chasing' the terminal speed and staying there with TWR=2 until you need to throttle up above TWR=2 just as you hit the next slab. Now make lots of such slabs and make them thinner and thinner, and you get a close approximation to Kerbin\'s atmosphere.

IV. Optimal fuel-saving ascent strategy for a single stage

So, you need 100% throttle at launch (and some SRBs would come in handy as well) until you reach terminal speed, then throttle down if you have to, to match the terminal speed for as long as you have fuel. To do this, from the equation of motion, will require a TWR= 2 + v_T²/(2H) where H is the scale height of the exponential atmosphere.

For most spacecraft, due to KSP\'s unrealistic drag model, the terminal speed which is the 'target' speed to shoot for as a function of altitude is given roughly by v_T ~ 97.3 exp(+altitude / 10183m), so for example:

Altitude (m) Target speed (m/s) - rounded to nearest 5 m/s

0 97.3- you\'d better get a move on!

500 105 - you should be catching up by now!

1000 110

2000 120

3000 130 - usually the first benchmark I have time to look at

5000 160

8000 215

10000 260

13000 350

15000 425

16000 470 - most of us are at least thinking about staging and pitchover maneuvers by now!

..

32000 2250 - this is equal to the orbital speed, so by now you should be pushing over hard for orbit, and air drag is not as important above this altitude.

The resulting ascent profile looks very similar to those that use optimal control theory to solve the Goddard Problem, e.g. Figure 2 in this paper (thanks jebbe for the link): http://cas.ensmp.fr/~petit/papers/ifac2008/0637.pdfComments welcome, if you\'ve read this far!

[PakledHostage]

Glad you learned something new - one of the reasons I submitted this challenge was so we could learn from each other, and I think this has been a success despite the low participation.

Although I\'ve already extensively used what I learned in this challenge, I\'ve been distracted by other challenges... I\'ll give it a try when I next get a chance to play the game. I might not get that chance until tomorrow evening though. I\'ll post back with my results.