What is the numerical procedure to find the next closest approach between two bodies on keplerian orbits?

Question

Given two bodies on Keplerian orbits, what is the procedure to find the next closest approach of the two bodies? The bodies would be effectively massless objects like two spaceships orbiting around a uniform central body. The orbits can be given by state vectors or the six keplerian elements or any representation of elements that leads to the easiest solution.

Note that by "next closest approach" I mean the immediate next time when the orbits switch from getting closer to getting further apart (which may not actually be very close at all compared to the optimal closest approach of the two orbits).

I expect it is a numerical root finding problem followed by a second derivative check of some sort.

I guess one simple way to do this would be to write a function which takes the state vectors, propagates them to time t and then calculates the distance, then just wrap that with Brent's 1-dimensional minimization algorithm (e.g. Matlab's fminbnd), then have a heuristic to find a bound range that the minimum is in or expand the search radius if it fails. It'd be nice to find a method with better convergence properties than that though.

Answer 1

Partial answer only: I hope a more specific answer will also be posted!

---

I see you've nicely updated your question. It's certainly still on-topic here but I think you can also abstract it just a bit and ask a related question in Math SE or SciComp SE as well.

A math question might look like:

Given two period functions $g_1, g_2$ of the independent variable $t$ how can all minima of some parametric function $f(g_1, g_2)$ be found within a range $t_1 \le t \le t_2$?

Analytical approaches

If you know $g_1, g_2$ and $f$ analytically (in this case $\mathbf{x_1}(t), \ \mathbf{x_2}(t)$ and $|\mathbf{x_2} - \mathbf{x_1}|$) then you can just use the chain rule to differentiate analytically, and then apply all that the fields of mathematics and numerical technique have to offer to look for zeros.

But for Keplerian orbits there are no simple analytical solutions for $\mathbf{x}(t)$.

This answer to What is the analytical closed-form solution of the two-body problem to verify its numerical integration results? explains that there is an analytical solution for the inverse problem, find the time for a given $\theta$ (or the two times for a given $r$) within one period, and that may or may not be helpful.

Infinite series

There are however infinite series solutions for $\mathbf{x}(t)$ and for mildly eccentric orbits you could consider using the first few terms as analytical solutions then proceed to zero finding.

Advanced techniques

• This just in: Kepler's Goat Herd: An Exact Solution to Kepler's Equation for Elliptical Orbits

A fundamental relation in celestial mechanics is Kepler's equation, linking an orbit's mean anomaly to its eccentric anomaly and eccentricity. Being transcendental, the equation cannot be directly solved for eccentric anomaly by conventional treatments; much work has been devoted to approximate methods. Here, we give an explicit integral solution, utilizing methods recently applied to the "geometric goat problem" and to the dynamics of spherical collapse. The solution is given as a ratio of contour integrals; these can be efficiently computed via numerical integration for arbitrary eccentricities. The method is found to be highly accurate in practice, with our C++ implementation outperforming conventional root-finding and series approaches by a factor greater than two.

• Towards Automated Satellite Conjunction Management with Bayesian Deep Learning

[...]For this reason, development of automated tools that predict potential collision events (conjunctions) is critical. We introduce a Bayesian deep learning approach to this problem, and develop recurrent neural network architectures (LSTMs) that work with time series of conjunction data messages (CDMs), a standard data format used by the space community. We show that our method can be used to model all CDM features simultaneously, including the time of arrival of future CDMs, providing predictions of conjunction event evolution with associated uncertainties.

See also:

• How might one approach using AI (convolutional neural network) to predict collisions in orbit?
• What DID NASA do with "the world's largest quantum computer"? (re question from 2013)
• https://www.seertracking.com/

YouTube:

• The Space Debris Apocalypse | Amber Yang | TEDxJacksonville
• This 19-year-old can keep astronauts safe from space junk

Lissajous patterns and repeating/non-repeating patterns

It is notable that the Wikipedia article for parametric equation begins with a discussion of the butterfly curve and has a whole section about Lissajous curves.

@RogerWood points out that it's important to decide if you can treat the ratio of the two periods as a reasonably sized rational number (e.g. 42:137) or if they are unrelated in any simple way and basically independent.

Who does stuff like this?

Space command and any satellite-owning or dependent organization does!

When the minimum distance between two Earth-orbiting satellites gets too close for comfort, meaning several times the uncertainty in each satellite's orbit, it's called a conjunction.

With zillions of satellites, calculating potential conjunctions in the future and determining when action needs to be taken (e.g. generate tweets and perhaps move some satellites around) is big business.

See the following for more on how this is done:

• Celestrak: SOCRATES
• Algorithmic methods or techniques to find conjunctions in large ensembles of state vectors?
• Algorithmic methods or techniques to find conjunctions in high N Keplerian element ensembles?

But what if I want to try solving this myself?

Define your problem in a practical way.

Do you care about every minimum, or only those that are actually close as in conjunction-finding?

If they are on opposite sides of the central body, then a short time later they are not going to be right next to each other.

I think some kind of divide and conquer approach may be helpful. Break up time into boxes small enough that they contain zero or one minimum but not two. Then divide those in half again. Check the derivatives at the edges, if it changes sign, then maybe a minimum is nearby.

What's next?

I still think that a carefully written question in either Math SE or SciComp SE will be quite fruitful.

Answer 2

This is an answer, but definitely not a good answer. If, $d$, is the separation between the two satellites then we can use Newton's method on the time derivatives $d'$ and $d''$ to iteratively find the minimum.

Given the 12 orbital elements, we first find the separation, $d$ and the relative velocity and acceleration between the two satellites. We then resolve the components parallel to the direction of $d$ and orthogonal to $d$ using the notation ($v_p,v_o$) and ($a_p,a_o$).

For the first derivative of $d$, we have simply $d'=v_p$. The second derivative $d''$ obviously includes the acceleration component along $d$, but it is also a function of the orthogonal velocity component that acts like an acceleration in causing the separation to increase quadratically. The expression for the second derivative is thus $d''=a_p+(1/2)v_o^2/d$. If we are approaching a minimum, $v_p$ must be negative and $a_p+(1/2)v_o^2/d$ must be positive. However, $a_p$ is always negative since both satellites are continually accelerating towards the central mass. So the relative speed must be at least $v>\sqrt(2a_pd)$ at the point of closest approach.

Close to a minimum of separation, $d$ is probably a reasonably quadratic function of time and Newton's method will work very well. At other positions, all hell will break loose! Unfortunately, none of this avoids the need for a good method of propagating the orbits in time (per @uhoh comments).

Answer 3

This isn't an actual answer, but you may find it of interest. ;)

As uhoh has mentioned, the central difficulty to an analytical solution of this problem is that there's no way to invert Kepler's equation using elementary functions. The equation is

$$M = E - e\sin E$$

where $M$ is the mean anomaly, $E$ is the eccentric anomaly, and $e$ is the eccentricity.

It's easy enough to solve it using a few rounds of Newton's method (unless $e$ is very close to 1), but that's not much use in an analytical approach, and your problem needs to solve it twice, once for each body. The traditional analytical approach is to approximate it using a few terms of a series; several options are given in the Wikipedia article.

The alternative is to just search for solutions numerically. If you know the time of the last closest approach you can make a reasonable initial estimate of the time of the next one by simply adding the relative synodic period of the two bodies.

Let $T_1 < T_2$ be the periods of the bodies. Then the synodic period is

$$T_s = \frac{T_1T_2}{T_2 - T_1}$$

If both the orbits are circular, this gives the exact solution. For eccentric orbits, it gives you the mean time between closest approaches.

Note that the relative velocity of the two bodies is zero at the closest (and furthest) approach, so the relative distance doesn't change much around that time.

---

I assume that you want a general solution for this problem. OTOH, if you want to solve it for actual Solar System bodies, you can get JPL Horizons to do the orbit calculations. It can easily give you the distance (and its rate of change) between any pair of objects in its system, and it knows a lot of objects: 1130203 asteroids, 3757 comets, 209 planetary satellites, 8 planets, the Sun, L1, L2, select spacecraft, and system barycenters. 

Here's a graph of the distance between Mars and Mercury over the period 2000-Feb-04 to 2003-Dec-30. I chose them because they are the planets with the most eccentric orbits. Mercury's eccentricity is 0.20563, Mars's is 0.0934. Their periods are 87.9691 days and 686.98 days, respectively, so their relative synodic period is ~100.888 days, which is ~1.1469 Mercury periods and ~0.1469 Mars periods. They had a close approach within a day of the starting time of the graph. The interval between the vertical date lines is 100 days.

As you can see, the period between successive close approaches is quite close to 100 days, although there is a bit of drift.

---

Here's a link to a live version of the Sage / Python program I used to create that graph. It can plot in SVG and PNG formats. You can use it on any objects that Horizons knows. It can be used to search "manually" for solutions to your problem, and perhaps it will be useful for validating other solution methods.

Please see the Horizons documentation for full details on specifying body IDs, times, and time steps. Briefly, a number from 1 to 9 identifies a planet barycenter (including Pluto), eg, 1 is Mercury, 4 is Mars. Append 99 to specify the body center. You need to prefix the observation center with @, otherwise the number is treated as the ID of an observatory. If you type a string into the target or center fields, Horizons will respond with a list of IDs that match that string.

Horizons accepts numerous date and time formats. To input a Julian day number, prefix it with jd.

The datestep parameter of my program sets the spacing of the vertical gridlines in the graph. A datestep of 10 means those lines are drawn for every 10 data points.

To reduce the number of Horizons requests, the program caches the last 3 sets of data that it fetches. If you make "cosmetic" changes to the graph without altering the target, center, or time parameters, the old data is reused.

Answer 4

I would try a numerical approach using basic engineering approaches. First set up 2 arbitrary keplerian orbits using r and v data or elemental data. If orbits of actual objects are of interest use accurate solar system data.

Second, pick arbitrary starting positions for each object.

Third, set up computation of the distance between the two objects in cylindrical (or spherical) coordinates or by transformation to cartesian coordinates.

Fourth, Set up procedrue to propagate the objects in orbit for a short time period (most likely a few minutes), then recompute the distance between them.

Fifth. run (propogate) the system for a few computer hours (likely tens of thousands of orbits), developing a log of the time and distance between the objects and other parameters that might be of interest in analyzing the results.

Six. Display the distance parameter vs time and examine for minimum points and time between minimums and other parameters of interest. If there is insufficient data, run for a few hours longer.

Working on a closed form solution for this does not appear to be worth the effort.