Goddard Problem Solver

[Racer319]

Hey everyone! I have some free time, so I'm writing a program that will take parameters for a ship, planet, and atmosphere and return the ideal launch path. Now that I'm starting to code the physics for it, I have a few questions about the drag force. As far as I know, in the real world drag varies with velocity squared at high speeds and just with velocity at low speeds. I'm guessing in KSP it's just v^2 the whole time, but does anyone know for sure? Also if you have any advice or suggestions I'd be happy to hear it.

[tavert]

Will this be a purely vertical ascent, or a more complicated version including a gravity turn to orbit, specific impulse varying with atmospheric pressure, hard atmosphere cutoff, variable throttle with an upper limit, etc? How do you plan to do the optimization?

There have been a number of discussions on this topic, you'll have to be very careful how you go about trying this and it will require some sophisticated numerical algorithms if you're trying to tackle the full version of the problem. Read through this recent thread for some guidance: http://forum.kerbalspaceprogram.com/showthread.php/46194-I-need-someone-help-me-do-some-math-for-launch-optimization

The real-world drag equation says drag is proportional to 0.5 times ambient density times speed squared times cross-sectional area, but the proportionality coefficient (called the coefficient of drag) is not actually constant, it varies quite a bit as a function of mach number, angle of attack, and potentially several other factors. Anything you've seen that says drag is linear with speed at low speeds is not referring to fluid drag, it's referring to mechanical losses or viscous damping or rolling resistance for a car, which are approximately linear trends that dominate at low speeds for something like a car. For a flying vehicle at low speeds (like a hovering quadcopter or Grasshopper), disturbances like wind gusts dominate and drag is not terribly significant.

In KSP, the drag equation is a little different. It's still proportional to 0.5 times ambient density times speed squared, but cross-sectional area is replaced by mass so the drag coefficient has a different meaning. The drag coefficient is the mass-weighted average of the drag coefficients of the parts that make up a craft - these coefficients are 0.2 for most parts, but slightly higher for solid rocket boosters, rover wheels, docking ports, and a few other things, and much higher for deployed parachutes. There's also an extra factor of 0.008 in the drag equation for some strange reason. See http://wiki.kerbalspaceprogram.com/wiki/Atmosphere for more information, although don't trust the 1.2002 number for density at sea level. As far as I can tell that appears to be the density at 76ish meters above sea level on the pad at KSC, but if you go all the way down to actual sea level then density's a little higher, around 1.223 kg/m^3.

Edited August 22, 2013 by tavert

[Racer319]

I'm planning on finding a 2D solution that accounts for the gravity turn, but right now it's 1D. I'm using 4th Order Runge-Kutta for integration, specific impulse, gravity as a function of altitude, atmospheric density, etc... basically I want to fully solve the problem.

I haven't completely thought it out yet, but I think I'll use some form of a gradient descent algorithm for optimization as it scales nicely with lots of variables.

About the drag, I was referring to Stokes Drag at low velocity. From past research, I believe drag transitions to velocity squared somewhere around 10-20 m/s, for example: a basketball's drag will probably vary with velocity, while a baseball is velocity squared. I know how to handle these equations pretty well, but supersonic and hypersonic effects are beyond me. It sounds like the game currently ignores effects at transonic speeds, thank god for me. That will make things much simpler.

It sounds like the KSP drag equation is F_d = (C_d)(v²)(p)(some constant)

where:

C_d = drag coefficient

v = velocity

p = atmospheric density

and it just ignores mach effects, cross section, and angle of attack? If so, I will need a value for C_d for each test ship, or at least a decent estimate. It sounds like a list of parts might be enough to calculate this? Let me know if you think that could work.

Eventually I think I can make this work with multi-stage rockets, but that's probably quite a long way off.

[Duxwing]

You may want to look into MechJeb, which can take a rocket up to an arbitrary orbit or down to an arbitrary landing site provided sufficient delta-V.

-Duxwing

[tavert]

Ah yeah, when things are dominated by viscosity and laminar skin friction rather than dynamic pressure I suppose it does end up linear. That's at pretty low Reynolds numbers, below the velocities we tend to care about. And in KSP you only have drag from dynamic pressure, no compressibility effects or angle of attack dependence (as long as the craft doesn't have any wings or control surfaces, and you aren't using the FAR mod).

It's not really Cd, it's more like Cd*m. Mass is used in place of cross-sectional area. Most ships have drag coefficient of 0.2, unless you have a lot of solid boosters or aerospikes or aircraft cockpits (the latter two are lower-drag than average).

Gradient descent doesn't handle poorly-scaled problems. It also isn't the cleanest way to approach problems with inequality constraints - you can use projection which is simple enough when your inequalities are just upper and lower bounds, but otherwise projection can be nontrivial. Applying constraints on the final conditions will also be problematic for a simple gradient method. Convergence of gradient methods isn't great, I recommend a Newton-based interior-point method instead. These are complicated enough to get right that it's best to use existing libraries that are out there.

And with all gradient methods you need derivatives of your dynamics (plus second derivatives for Newton methods). It's simple enough to write a simulation program given a pitch and throttle profile, but much more time-consuming to also return the gradient, Jacobian, and Hessian matrices that an efficient optimization algorithm requires. Finite differencing is terrible for any problem with more than a handful of variables, so don't do that. There are optimization modeling languages out there that are much more suited to this task than just writing up a forward simulation in C or Python. You ideally want the modeling language to have automatic differentiation capability so it can generate all the information required by the optimization algorithm without you having to code the derivatives up by hand. AMPL or OpenOpt are probably the best options, though AMPL isn't free and I'm not sure how well OpenOpt scales to large problems.

Edited August 22, 2013 by tavert

[magnemoe]

You may want to look into MechJeb, which can take a rocket up to an arbitrary orbit or down to an arbitrary landing site provided sufficient delta-V.

-Duxwing

Mechjeb just launch along an set path, 5 km start turn 70 km end turn and reduce trust so speed stays terminal. If you launch from other places you need to change the launch profile for better performance.

[Racer319]

tavert thanks for all the helpful suggestions! and Duxwing i obviously know what mechjeb is, but it's ascent profiles are not perfectly optimized, and depending on your ship design they could be quite bad, so I want to try to improve on this.

tavert: I have been writing the program in Java, and I'm not sure why you recommend other languages. It seems like the most complicated math I will have to do is to solve differential equations and find derivatives, which I can estimate very closely with simple algorithms, so I don't really see the need the need for a higher level language. Although this might just be because I don't have any experience with them.

[tavert]

Java's not that commonly used for low-level numerically intensive work. If it's what you're most comfortable working in go for it, there are Java interfaces to efficient optimization solvers, but none of the serious ones are written in Java themselves. Just like with Python or Matlab, the most important numerical kernels should be written in a compiled language and interfaced with a Python/Matlab/Java wrapper. The languages that I recommended are specifically designed for optimization as opposed to general-purpose programming or technical computing, so they handle all of the derivatives and the optimization problem formulation and interface with solvers automatically.

If you try a naive simple algorithm like gradient descent on a difficult large-scale sparse nonlinear optimization problem, you'll have to be extremely lucky for it to work well. The optimization algorithms that I recommend which work well for these problems are nonlinear interior-point methods. I can give you a few references to read about the basics of the algorithm, from the simple statement of the method you'll quickly realize you need a numerical linear algebra library that can solve large sparse symmetric indefinite linear systems at each iteration of the optimization algorithm. You'll find a few of these, most of which are actually still written in Fortran. But if you're not careful about the details of the algorithm in terms of inertia perturbation, primal-dual line search, and barrier parameter update, it won't work well. Use a proven library like Ipopt (https://projects.coin-or.org/Ipopt, it even has a Java interface if you must) that has already done the hard work, been extremely careful about the details, and been demonstrated to work reliably on a wide variety of problems.

Getting the derivatives right is the time-consuming part. Not necessarily difficult, but annoying to code properly. You can easily write up a finite differences approximation, but you'll have numerical accuracy problems and finite differencing is not viable past a handful of variables (pitch and thrust multiplied by number of time steps, plus intermediate states like mass, altitude, horizontal and vertical velocities), as I mentioned. You need the derivative of your objective function and each of your constraints with respect to each individual variable, and the second derivatives with respect to every pair of variables. Looks like there are a couple of Java automatic differentiation libraries you could use that might work here, JAutoDiff or ADiJaC or Leibniz or JAP etc, never used any of them though so can't say which or any to recommend.

Or have a look into this tool I just ran across, it's an optimization modeling language for Java, might be your best bet: http://ait.upct.es/~ppavon/jom/index.php. It doesn't appear to have an exponential function yet, but oddly it has ln so you could accomplish the same thing with an equality constraint. I bet one message to the author could see exp added pretty easily.

Edited August 24, 2013 by tavert

[tavert]

Did this really need to be moved to Add-on Development? There's no add-on here, just mathematical discussion quite similar to about 3 or 4 other recent threads in Science Labs. This isn't the kind of thing you would need to run real-time in an add-on anyway...

[Mekan1k]

Actually, this would be an incredibly useful addon. If this program can be written, (since we can ignore the 3-body gravity problem and since the in-game physics are perfectly understood) incredibly efficent flight-paths can then be produced.

Think of all the fuel you will save!

[tavert]

Well once anyone comes up with a piece of code that can actually solve the difficult mathematical problem, then we can talk about using it in the game. It's entirely possible we could look at an optimization-based ascent planner and integrate it into MechJeb for example, but we're not at that stage yet.

[Vector]

I'm taking a crack at this problem.

So far I've got the forward modeling working, meaning simulation is pretty close to experimental results. For now I've left out Coriolis force and the rotation of Kerbin but everything else is in place.

The representation of the problem is a huge factor in the optimization, more important than (or really an integral part of) the optimization algorithm. I've got two or three failures so far where my representation/optimization failed to converge or where the discrete approximation made too much of a difference.

All this is in a brand new C# program so hopefully there will be avenues to construct a new plugin or integrate with MechJeb.

[tavert]

Problem formulation is important, but don't underestimate the importance of the optimization algorithm (or try to throw something together haphazardly instead of using proven existing tools). You can formulate a problem perfectly but if you try to solve it with a terrible optimization algorithm it won't work properly.

Anyway, until anyone comes up with working code that can solve the problem reliably, this is more of a mathematical discussion than a plugin discussion. I'm still waiting for an explanation why this thread was moved out of Science Labs, but there's a great deal of similar discussion happening in this thread: http://forum.kerbalspaceprogram.com/showthread.php/46194-I-need-someone-help-me-do-some-math-for-launch-optimization

I think we should keep all the discussion in one place. There are 3 big families of approaches to this problem, and people have been trying each of them: pseudospectral methods, local discretization (as in forward simulation, but with additional gradient information required by an optimization algorithm), and optimal control Pontryagin methods that express the solution as a boundary value problem.

Edited August 29, 2013 by tavert