mockingbirdnest / Principia

𝑛-Body and Extended Body Gravitation for Kerbal Space Program
MIT License
769 stars 69 forks source link

Principia endpoints for mods implementing closed-loop guidance #1659

Open lamont-granquist opened 6 years ago

lamont-granquist commented 6 years ago

To do accurate state vector projection this would be useful, and would be very useful for things like PEG. The general problem is to take state vectors (r0, v0) at t = now and integrate them forward to (r1, r1) at t = now + deltaT (not under thrust).

This would replace the "PRECISE PREDICT" and "SUPER-G" functions in the PEG-4 routine:

https://ntrs.nasa.gov/archive/nasa/casi.ntrs.nasa.gov/19760024151.pdf

Pages 16-19 and 75-81 (I think you can ignore the yapping in there how they used it with J2=0, I don't need to be able to switch off J2).

Ultimately PEG might call this three times for burn-coast-burn functions (currently it needs to calculate coasting trajectories once per guidance cycle). The burn arcs also need to call this in order to do gravity integral averaging (PEG calculates gravity integrals and thrust integrals separately and approximates gravity integrals over a thrust arc by a "nearby" non-thrusting arc which is a first order approximation). For PEG's usage, the accuracy does not have to be to be to machine-precision, I'm not sure of the order of magnitude of errors that would be acceptable but my guess is that 1e-4 would be sufficient, and should likely be another parameter. PEG is already tunable so that while PEG might make three calls per guidance cycle, the guidance cycles can be run slower than 50 Hz of physics ticks and can be set to run at 1 Hz for users with potatoputers. The order of magnitude of the time of integration is likely to be order of a dozen minutes or so and not multiple orbits. Hyperbolic orbits could conceivably have coast phases to target SOI exit velocities so might have coasts on the order of days (this is more of a nice-to-have though and would be a very negotiable requirement if its hard to do).

On subsequent guidance cycles the values that would be used would be close to the prior values so it might make sense to hold onto objects with internal state of the integrators to use the previous solution as a first guess to the next solution.

egg sez:

egg: lamont: anyway; a  FlowBodyCentric(centre_index, world_body_centred_degrees_of_freedom, t_initial, t_final) seems like a sensible entry point, file a feature request and we'll see what we can do. One problem is that you need to have the plugin object, which is owned by the ksp plugin adapter if it exists and which I'm really reluctant to expose too much... maybe I could have a static variable with the current plugin, and
[4:35pm] egg: have your entry point use that one, but that sounds like it might be brittle?

A complication is that I need to be able to call these routines from below the inverse rotation threshold and input inertial r,v and get out inertial r,v without getting trolled by the inv rotation the way the Orbit class does.

lamont-granquist commented 6 years ago

Per discussion on IRC, the r,v that I need must be ECI with Z being earth north axis up and left handed. Beyond that I get confused, but r = vessel.CoMD - vessel.mainBody.position and v = vessel.obt_velocity should be the r,v of where the active vessel is.

eggrobin commented 6 years ago

(that's ECI but with KSP's ineffable occasionally-rotating World axes, we know how to handle that better than we know how to explain it)

lamont-granquist commented 6 years ago

egg explains to me that those are still in rotating axis world coordinates and that hurts my damn brain a lot, but those are the vectors that MechJeb operates on and everything seems to work.

eggrobin commented 6 years ago

Additional endpoint that @lamont-granquist would like for guiding with respect to the flight plan: Given a position q (in some frame, presumably body-centred inertial), and a vessel with a flight plan, return the degrees of freedom of its first future planned closest approach to q.

Confirmed with @lamont-granquist that he is fine with obtaining the pointer to the Plugin via reflection on the adapter, so no need for static horrors.

lamont-granquist commented 6 years ago

yeah q is body-centered-intertial. and i need the unpowered target orbit without any finite burn effects, and would need the unpowered orbit solution for times potentially before the end of the burn.

and like i yapped about quite a bit on irc, i only need a single 'reasonable' solution. the initial guess for 'q' from PEG could conceivably be quite pathological -- consider the ideal (not-N-body) case of a circular equatorial orbit target with 'q' being over the north pole. the solution is then literally any point on the orbit. but PEG won't care and the t=now solution for the orbit would do nicely to start moving PEG towards that value, from which it should then hunt down the solution.... not that PEG should ever do anything quite that substantially pathological, but the point is more that in any pathological double-valued case, PEG just needs a value to move towards a solution. if there were ever a case where PEG really winds up converging towards something pathological, then that would be a high non-convex optimization problem and i'm okay with PEG doing nonsense for that. need a next generation trajectory planner + burn executor for that.

eggrobin commented 6 years ago

Ah, that makes things a little tricky, since we then need to integrate backward from the end of the burn; possible, but less off-the-shelf.

lamont-granquist commented 6 years ago

yes, so PEG does approximations that may very well converge to times before you've ended your finite burn simulation and started the coasting orbit. PEG will also start with guesses that are just 1 second ahead from where the burn starts (actually i start running PEG a few seconds before the burn is even supposed to start and just let it run). so it needs those negative timescales.

lamont-granquist commented 6 years ago

so i'd like to grab this orbit object sometime early (like whenever the user pokes the 'execute node' button which might be days in advance). i will NOT want backpropagation to that timepoint. at some point i'll fire up PEG about 10 seconds before i want to burn. at that point i'll want backpropagation to "now", and will never need backpropagation to the past i think.

lamont-granquist commented 6 years ago

oh yes, C++ and garbage collection -- forget i mentioned object handles.

lamont-granquist commented 6 years ago

Actually I think I would like to see the API without backtracking -- just give the first point in the trajectory if it is the closest point. And then we can go from there and determine if backtracking is actually necessary or not. That would simplify what you need to do in order to deliver the endpoint. Then I can plug it into PEG and see how it explodes or not. If it doesn't explode then there's obviously no point in doing the extra work. I still suspect that it will have issues in certain circumstances, but it'd be better to actually verify that it does. I've been wrong before.

eggrobin commented 6 years ago

Status update: An untested draft of the ExternalGetNearestPlannedCoastDegreesOfFreedom function is available in my lamont-endpoints branch.

@lamont-granquist has been able to find it by reflection, but has yet to use it for guidance.

eggrobin commented 6 years ago

Longer term feature request: calculus of variation voodoo for arbitrary guidance, see https://logs.tmsp.io/principia/2018-01-30#1517337985-1517337551;

This requires implementing Runge-Kutta methods (probably embedded explicit) with inhomogeneous components (so operating on an std::tuple, or perhaps an std::tuple of std::vector, rather than an std::vector).

eggrobin commented 6 years ago

Note that it may also make sense for us to do the shooting ourselves, in which case we need a slightly better root finder.

Eventually we could even look at fancier methods than shooting.

eggrobin commented 5 years ago

Note: this feature request tracks the Principia-side prerequisites for "tight integration with principia's predictions for N-body accuracy" in Phase II of MuMech/MechJeb2#1061, and possibly "support principia node execution without tight integration" in Phase I (the tap-dance of MuMech/MechJeb2#1062 is not acceptable user experience).

Given that MechJeb is moving away from PEG, we may want to drop the existing API.

Before we start thinking about the actual details of the numerical integration of primer vector guidance with Principia-based shooting (optimized by the caller), the first question is the computation of the residuals (this is relevant to phase I too, since usable Principia node execution would entail having some form of residual from the target Principia trajectory).

For Principia flight plans, only Principia can know about the relevant reference frame and target trajectory, so Principia should compute the residual (getting an orbit that looks the same in ECI when one seeks a specific EMB orbit is pointless).

However, for some manoeuvres, some degrees of freedom should be left to the optimizer; for instance, a circularization should leave mean longitude free, i.e., should let the optimizer choose where it wants to be along the target circular orbit; a naïve integral over time of some metric between the trajectories will not achieve that.

Some metric of the distance between the lines in the reference frame (independently from time) may be a better approach: if the trajectory looks the same in the right reference frame (but is shifted in time), assuming that the correct reference frame was chosen, it likely fulfills the intended purpose.

@lamont-granquist, what do you think?

lamont-granquist commented 5 years ago

Well the 5-constraint free-anomaly residual ("I want to hit this orbit but I don't know where") looks like this:

https://github.com/lamont-granquist/MechJeb2/blob/4b5290b9a3f9390248c746cccec45361107dcd14/MechJeb2/Pontryagin/PontryaginBase.cs#L467-L472

What you wind up doing is computing the miss in the orbital angular momentum then take the x and y vectors of the miss of the ecc vec, then combine that with the transversality condition, which comes from the calculus of variations, and is the value of the free-space/zero-thrust hamiltonian set to zero (all my units are mu = 1, so there's probably factors of mu lurking around in there which are hidden).

Because you drop the z component of the ecc vec you need to worry about equatorial orbits with h_z = 0 and flip the coordinate system if that happens:

https://github.com/lamont-granquist/MechJeb2/blob/4b5290b9a3f9390248c746cccec45361107dcd14/MechJeb2/Pontryagin/PontryaginBase.cs#L449-L458

(I should probably if (Math.Abs(h[1]) < 1e-15 ) or something there as well, since I imagine it doesn't like it much if it just gets very close to that as well.

PEG like residuals for launches are even a bit weirder:

https://github.com/lamont-granquist/MechJeb2/blob/4b5290b9a3f9390248c746cccec45361107dcd14/MechJeb2/Pontryagin/PontryaginLaunch.cs#L74-L80

How that all translates to N-body IDK, particularly since IDK about using calculus of variations in the low-energy regime.