Life's a series of trade-offs. Today instead of cleaning the bathrooms I continued to work on optimal rocket trajectory planning. This version represents a switch from Haskell to C++, which is where I'll need it to be ultimately anyway.
New things I'm trying:
- Eigen linear algebra library. It's the greatest! It has plugged in very smoothly and just works, at least so far as I have used it. There are lots of constructor and operator options, so the resulting code is pretty close to the Haskell in conciseness.
- GLFW is also pretty great. I've used GLUT for years as my test application framework. GLFW fills much the same space but is actively developed, and can be statically linked (I hate dangling DLLs). I was able to trivially get input from my 360 controller, which was exciting.
- Learn OpenGL tutorial web page. I've used OpenGL 2.0, more or less, for years; this represents my initial forays into modern OpenGL. The planet is a flat piece of geometry that computes its lighting and anti-aliasing in the pixel shader, for instance. The lighting is in linear space; the crisp terminator on the planet is the give-away. Fun stuff, apart from the runtime compilation of text source code, which seems half-baked to me.
The current approach is to build an initial guess trajectory out of two equal-duration parabolic arcs. I chose two as being the minimum that could work. As it turns out, the nonlinear nature of the gravity field means that two segments do not work very well in some cases; the resulting trajectory does not bear a close resemblance to the trajectory you'd get with more segments. The process of ensuring feasibility tends to make it slew across the gravity field a bunch to get a usable average gravity vector.
The state vector includes the segment duration, the thrust vector for each arc, and the position and velocity after each arc. There is a set of 11-13 constraints. Each segment has four constraints asserting that the position and velocity variables integrate correctly across that interval with respect to each other and the acceleration. There are end-condition constraints, in this case specifying that the rocket is at ground height with zero velocity. Finally there are constraints that the acceleration vectors are less than the maximum rocket thrust. These are inequality constraints so they may or may not be present depending on whether the rocket is thrusting at maximum.
The initial guess is not feasible, in my case; I haven't made all the variables agree such that the constraints are all satisfied. So the first thing to do is to try to get to feasibility. This is done by measuring the error of all the constraints; getting the first derivatives of the constraint functions with respect to all the state variables; and then solving a linear system to find multipliers for the constraint gradient vectors that, when applied to them and summed up, will equal the constraint error vector. The sum-of-scaled-gradient-vectors is the step to restore feasibility. Of course this is a linear approximation, and the constraints are non-linear, so we won't actually be feasible. Rinse and repeat and hope for the best.
What works pretty well is to do a few iterations to get a close-to-feasible trajectory, and then subdivide it into four segments and iterate some more, then subdivide into eight segments and iterate some more. This gets to a nicely feasible trajectory in only a few iterations. The screenshot shows an example result from doing this.
I think next I may try switching to cubic arcs instead of parabolic arcs. Gravity and rocket thrust would be evaluated at the segment endpoints where I have position and velocity defined, which should be a bit simpler. In my parabolic arcs the rocket thrust is assumed to be constant over a segment, and gravity has to be assumed to be constant as well. I'm evaluating gravity at the midpoint between the two end positions, and then using that.
To get optimal trajectories I will need second derivatives of the constraints, and matrices with roughly double the dimensions; the fabled KKT system. I've gotten that working for a simpler problem over in Haskell; just have to get it over here with my planetary gravity field and so forth. More derivatives means more ways for things to explode, though.
