Orbital parameters of more than two bodies I am not even entirely sure that the question makes sense except in certain simple cases, but let me explain what I am trying to do, but with the question first.
The question is essentially:
Given n bodies with known position, mass, velocity vectors, how do I in a meaningful way calculate what they are each orbiting right now? It does not matter than the orbits will commonly change over time for chaotic systems. I just need the current best orbit definition.
To better explain my issue...
Assume two gravitating bodies isolated in space. They will each orbit their common bary center and you can define the position of that barycenter as well as the parameters for each bodys orbit around it.
Now consider a case with three equal massed bodies orbiting their common center of mass in an equilateral triangle. Again you can calculate their barycenter and orbital parameters (however unstable to any pertubation). 
Now consider a case like Sun, Earth, Moon where the three bodies have a common center of mass, which they orbit, but where Earth and Moon more rapidly orbit their common barycenter. How do you calculate, given only position, mass and velocity vectors, that this system is best described as sun and the earthMoon barycenter orbiting their common barycenter while Earth and Moon orbit their barycenter?
I currently attempt to do this by looking at all pairs sunEarth, sunMoon, earthMoon and calculate their barycenter and orbital properties in isolation and then consider "the real" orbit around the real barycenter the one with the shortest period. That gives me the correct answer in the two body case and in the Sun, Earth moon case, where two bodies are very isolated, but it would fail with the tree equal sides bodies.
I also considered defining the gravitational potential equation for the entire system and looking at the gradient to see that in the Sun, Earth, Moon case, Earth and Moon are circling in a little local indentation while still being caught in the larger indentation caused by the sun, but I am not certain that will really help and it seems somewhat complicated to properly model when a body may be in such an indentation but with a velocity direction and magnitude which will let it escape regardless.
If the answer is that for the general case there is no answer, then I guess that is an answer as well and I can stop wasting energy on it. Any other thoughts about the validity of the question, when it certainly doesn't make sense and anything else will also be greatly appreciated.
 A: @DavidZ comment in the suggestions is correct: there is no unique, general solution to this problem.  Depending on the context, dozens of different heuristic approaches are used in simulations. This is part of the general 'N-body' problem.
The most literature is in the context of large N-body, cosmological simulations where the problem is called often referred to as 'group-finding' (or 'Halo-finding').  In the context of (proto)stellar systems, this topic per se doesn't get much attention (to my knowledge) because you don't really need to know secular groupings.
See also the review: Knebe et al. 2013 - Structure Finding in Cosmological Simulations: The State of Affairs
In general, the only way to do this would be to simulate the system over the time period of interest and measure whatever properties you want (i.e. if a body is constrained to a certain region, remains gravitational bound to another, etc). There are any number of approaches to get an effective treatment... but they are, usually, not based on any strong principles --- just what happens to work well.  There was a(n oddly) similar question yesterday which received a good answer that might be helpful.
Likely the best approach is to find systems which are gravitationally bound, perhaps hierarchically (e.g. moon and earth are bound, both are bound with the sun).  Grouping in 'phase-space' will likely also be important, where you consider particles which have similar positions and velocities (phase-space, in this context, is the 6 dimensional space of $\vec{x}$ and $\vec{v}$).
Because there is no universal 'solution', in the end, any practical solution must be tailored to the specific properties a project/simulation is looking for, and over what timescales.
