Lambert Solver for Lagrange Points

[MontBlanc2012]

For a while now, I've been thinking about the kinds of tools that one needs to navigate efficiently and effectivel in the vicinity of the Lagrange points of, say, the Earth-Moon system. To be sure, Orbiter has a number of lesser-known (but quite sophisticated) tools in the Orbiter MFD stable - notably Lagrange MFD v1.5 for Orbiter 2016 - which permit some basic navigation relative to the Lagrange points of the Earth-Moon system and the Sun-Earth system. But none of these tools allow one to design optimal transfers in 'cis-lunar space'.

Surely we can do better than this.


A comparison with two-body motion

As most users already know, Orbiter has been very successful in enabling tools that allow for the design of efficient interplanetary transfers using the now very familiar 'patched conics' approximation . In particular, IMFD's Course Program can be used to design efficient Hohmann-like transfers between planets. It is able to do this because it makes use of a Lambert Solver to construct a Keplerian arc between a trajectory start-point and a trajectory end-point with a specified time-of-flight. A Lambert Solver solves the two-point boundary value problem "Given a starting point A, an end point B and a transfer time dt, what initial velocity is required at the point A to make this transfer". For the Kepler problem efficient solutions to this boundary value problem can be computed. And ny judicious adjustment of the departure date and arrival date, one can construct an optimal, two-impulse transfer manoeuvre that minimises the amount of dV required to make the transfer.

But what about motion in the vicinity of a Lagrange Point? In the vicinity of the Lagrange points of the Earth-Moon system say, IMFD's Keplerian Lambert Solver built into its Course Program no longer works (because the gravitational field in the vicinity here is no longer Keplerian). And, yet, we do not have an equivalent Lambert Solver for cis-lunar space that we can use to carry out the same sort of trajectory optimisation that we would for the patched-conics model. For example, let's suppose that we are at some arbitrary position reasonably close to the L1 Lagrange point and travelling reasonably slowly relative to it. Let's also suppose that we wish to make a transfer to the Lagrange point that takes a specific time, dt. What velocity should our satellite have now in order to effect this transfer to L1? And what velocity will we have relative to the Lagrange point when we arrive? When is the optimal time to leave; and what is the optimal duration of our transfer? In the absence of a Lambert Solver for the restricted three-body problem, we struggle to answer basic questions of this kind.


A different kind of Lamber Solver

Evidently, the only way of answering questions of this kind is to have a Lambert-like Solver that has been designed to work in the vicinity of Lagrange points. Needless to say, I wouldn't be posting this note if a tool of this kind didn't exist.

Let's make things a little more concrete by considering a simple (but arbitrary) example. And because we are thinking about navigation in the vicinity of Lagrange points, we are going to set this example up in the rotating reference frame of the three-body problem such that the location of the Lagrange points in that reference frame is fixed. In this rotating reference frame, imagine that we have a satellite located 4,000 km from L1 towards the Moon (the 'x' coordinate in our rotating reference frame); 4,000 km vertically above the orbital plane of the Moon in its orbit around the Earth (the 'z' coordinate in our rotating reference frame); and 4,000 km 'ahead' of the Lagrange point (the 'y' coordinate in our rotating reference frame). In other words, the satellite starts off at some arbitrary point about 7,000 km from L1. Let's also suppose, for sake of example only, that we wish our transfer to L1 to take a rather leasurely 71.3 days. Of course, there is no particular reason for having a transfer of this kind, but we choose these number somewhat arbitrarily simply to work through an example in this note. In the context of this scenario, we can ask a number of questions:

* what is our transfer trajectory?

* what velocity do we need to have now to send us on this trajectory? and

* what velocity will we have when we finally arrive at L1?

* if we vary our time-of-flight, how do these things change?


These questions parallel those that we would ask of a standard Lambert Solver. And we use answers to those questions to drive the search for optimality. So, if we can answer these questions for the three-body problem, we should be able to carry out the same sort of optimsisation tasks that we would normally apply to a standard Keplerian (two-body) Lambert Solver. And since I do have a solver of this kind I can, of course, use it to answer these questions.


What is our transfer trajectory?

First of all, what is the transfer trajectory that we need to effect the transfer from our (arbitrary) starting position to L1 such that the time-of-flight is 71 days. Well, if we apply our new Lambert Solver to this problem, we find (without proof) that the three-dimensional transfer trajectory looks like this:

Clearly, the trajectory has a complicated structure. But, then again, this is a transfer orbit that takes 71.3 days (~ 2.5 lunar orbits of the Earth), so one might have expected that the transfer wouldn't be straightforward. So what's going on here in this trajectory. Roughly speaking, from the start point, the trajectory follows the stable manifold of Lissajous orbit around L1. Once 'on' the Lissajous orbit, the satellite circumnavigates the L1 Lagrange point six times before finally exiting the Lissajous orbit on the unstable manifold of that Lisaajous orbit that happens to pass through the Lagrange point arriving at L1 on the correct date. Of course, there is nothing particularly significant about the starting point or the duration of the transfer. They have merely been used here to demonstrate the capability of the Lambert Solver. As for the Lambert Solver itself, I'll focus on the mathematical details elsewhere but its foundations rest on sold maths and physics and, in effect, reduced the problem to an exercise in linear algebra. Although the end result is a trajectory solution that makes use of the stable, unstable and centre manifolds in the vicinity of the Lagrange point, the Lambert Solver does not solve the problem in terms of those structure.


What velocity do we need to have now to send us on this trajectory?

Naturally, we can find the velocity vector (in rotating coordinates) that we need to start ourselves on this trajectory. For our sample solution the x-y-z components of the velocity vector (in the rotating coordinate system) are as follows:

* x-direction - -34.65 m/s

* y-direction - +3.03 m/s

* z-direction - +3.19 m/s

Just as with a standard Keplerian Lambert Solver, we are able to calculate the initial velocity vector with ease. However, whereas Keplerian trajectories are meta-stable so that once on a Keplerian arc, and absent perturbations, one stays on the same Keplerian arc, for the Lambert Solver for Lagrange points, we have a situation where we will need to undertake periodic corrections in order to remain on the prescribed transfer trajectory. However, this kind of 'station keeping' doesn't present any particular problem since we can periodically re-run our solver during the flight to calculate the minor correction needed to remain 'on track' for the entire 71.3 days.


What velocity will we have when we finally arrive at L1?

It is just as straightforward to calculate the velocity components upon arrival at the L1 Lagrange point. For our example, for our test problem, we find that the speed of the vessel upon arrival is:

* x-direction - +8.75 m/s

* y-direction - +19.02 m/s

* z-direction - +18.31 m/s

Upon arrival at L1, and assuming one wishes to 'park' at L1, this is the velocity vector that would need to be nulled out. If we take the length of the velocity vector, we find that the dV required to park at L1 is 27.8 m/s.


If we vary our time-of-flight, how does this change?

We can also ask the question: is this solution optimal? Well, let's suppose that we only care about minimising the dV requirements at L1, given our starting point, what is the time-of-flight that minimises this dV? At this point, we can use the same kind of techniques that we would use when using IMFD's Course Program to optimise the transfer trajectory. And without going into the details, but simply using the time-of-flight as a control parameter, we find that if we target a time-of-flight of 69.5 days rather than 71.3 days, we can reduce the dV requirements to 'park' at L1 from 27.8 m/s down to just 15.7 m/s. For those interested, the optimal trajectory is shown below:

Of course, this optimality search was not intended to be exhaustive - but the point remains valid: once one has a three-body Lambert Solver, one has one of the basic instrument for finding optimal transfers in the vicinity of Lagrange points.


Where to from here?

By virtue of the trajectory planning tools that have been created thus far for Orbiter, Orbiter has largely been a vehicle for designing, optimising and executing patched conic trajectories. As such, it has been remarkably successful. But even though Orbiter includes a faithful rendition of Newtonian gravity for the Solar System, very few tools have been created that allow one to design, optimise and execute trajectories in the vicinity of Lagrange points. The Lambert Solver briefly showcased in those note goes some way to addressing that issue.

As for the Lambert Solver, itself. At the moment it exists, and it works - but there is a substantial amount of work to be done to make this into a tool that is useful by other in the Orbiter community. Nonetheless, this is the long-term goal of this endeavour and this note is but the first small step towards achieving that end.

 

[Thorsten]

Not sure how your solver looks like in practice, but feel free to use the numerics of LEO Targeting for your needs (or to contribute to the code) - I'm currently at the point where I toyed some with the Earth/Moon Lagrange point dynamics and plan to expand on it. Basically it does fits over brute-force numerical trajectory computations with arbitrary perturbations (3rd body gravity, J2 and J3,...)

 

[MontBlanc2012]

Thorsten

I haven't described the solver in detail in part because the maths is probably a little too much for most on this forum; and partly because setting out the details will take considerable effort.

But, in brief:

The solver is configured as a single-step very high-order symplectic discrete variational integrator using Gauss-Lobatto quadrature as per the procedure described by Marsden & West (Marsden & West 2001) reconfigured as the solution of a boundary value problem.

In effect, this fits a high-order (near mini-max polynomial) to the motion of the satellite and, as such is capable of calculating the trajectory of a body in a arbitrarily complex gravitational field with machine precision over multiple orbits in a single time step. These integrators are fearsome and powerful beasts.

The physical problem that is being modelled is the elliptical three-body problem in rotating, pulsating coordinates. In the vicinity of the Lagrange points, the equations of motion can be linearised and, in this case, and using the Marsden & West variational integrators, one can solve the two-point boundary problem in rotating coordinates using linear algebra after having specified the transfer time.

For the purpose of this note, I have restricted attention to the linear approximation. But this linear solution can be used as starting-point for finding highly accurate solutions to the full, non-linear elliptical three-body boundary value problem.

At the moment, the modelling and calculations are done in Mathematica. Although these can be ported to c/c++ without much effort, I'm not sure how this endeavour fits with your LEO targeting tool - although I'm sure that one could solve a multi-period perturbed Keplerian variant of the Lambert Solver for that case if one wanted to.

 

[Thorsten]

Okay, then it really seems you don't need any of my code. Good job accessing an interesting problem!

 

[ADSWNJ]

It would be interesting to see your solver mated with my LagrangeMFD. The code in LagrangeMFD is a 4th order symplectic integrator, well able to plot really eccentric orbits that you find around LP's, and able to do hypothetical burn predictions with really good accuracy over several thousand seconds.

What it lacks is the ability to find burns into really beautiful orbits.

An interesting solution would be to code your solver into an MFD, and then pipe results via ModuleMessaging into my utility to visualize the orbit and to execute the burns. It's non-trivial work though.

 

[MontBlanc2012]

ADSWNJ

I'm broadly familiar with your Lagrange MFD - I've even used it a couple of times to construct ad hoc, low-energy transfers from Each to the Moon and back again.

Although I'm quite pleased with the Lambert Solver I've introduced in this thread, I'm conscious that I have still a fair bit of work to do on the mathematical side of things before I would contemplate housing the thing in any form of MFD.

Having said that, I have in mind using this kind of tool to do just the thing that you have highlighted - i.e., "to find burns into really beautiful orbits." Two things are needed to do this. The first is to build a Lambert Solver of the kind described above; the second is to build a widget that describes mathematically the various kinds of Lagrange point orbits that one can have.

On the plus side, the mathematical foundations of the latter component has already been largely done by Lei et al (High-order solutions of invariant manifolds associated with libration point orbits in the elliptic restricted three-body system. Still, coding of Lei's work into a useful tool is still a fairly formidable (and time consuming) challenge.

All that I can see is that I'm moving in the direction that you have outlined - but it is going to take a while to get a point where it makes sense to have meaningful discussions about how best to house in the Orbiter MFD stables.

---------- Post added 04-14-18 at 05:14 AM ---------- Previous post was 04-13-18 at 10:24 AM ----------

To further illustrate what one can do with a Lambert Solver as outlined above (and to expand on comments in reply to ADSWNJ) consider the following scenario in the Earth-Moon system:

Imagine that we have a vessel initial at the point labelled 'Start Point' in the above diagram and it wishes (for some reason) to transfer to the planar Lagrange point orbit highlighted in red and finally achieving a final rendezvous with the target Lagrange point orbit at the point labelled 'A'. First of all, we need a way of describing the red target orbit; and second we need a way of describing the transfer from the Start Point to the rendezvous point 'A'.

For simplicity for this illustration, we can work with the linear approximation to motion in the vicinity of the Lagrange point, in which case the motion in the x-direction (i.e., in the direction pointing from the Earth to the Moon) is given by:

[MATH] x(t) = 2000.0 \,\cos (\omega\, t )+334766.0[/MATH] km

and in the y-direction:

[MATH] y[t] = -7173.0 \sin (\omega\, t )[/MATH] km

where [MATH]\omega = 0.523835\,\text{day}^-1 [/MATH] and where [MATH]t[/MATH] is measured in days. From this, description we can determine any point on the target Lagrange point orbit and, moreover, the velocity that we need to have in order to stay on the target orbit.

Lambert Solver solution for transfer time = 28.075 days
Given this description of the start and end-points, we can now feed this description into our Lambert Solver. To get this solver to work, we need to choose a transfer time. Let's try, for sake of argument, a transfer time from the Start Point to the rendezvous point 'A' of 28.07 days. The Lambert Solver solution (the think black line) is shown below:

This shows that the transfer trajectory dropping down to a lower amplitude Lagrange point orbit before exiting to meet the rendezvous point with an approach speed of 13.18 m/s.

Lambert Solver solution for transfer time = 28.966 days
Let's try a different transfer time of 28.966 days, say, and re-run our Lambert Solver to calculate a new transfer trajectory with the new transfer time. The Lambert Solver solution (again, the think black line) is shown below:

This shows that the transfer trajectory dropping down to a higher amplitude Lagrange point orbit before exiting to meet the rendezvous point with an approach speed of 14.78 m/s.

Lambert Solver solution for transfer time = 28.623 days
Now, we can play around with our transfer time to try and minimise our approach speed at the rendezvous point 'A'. It doesn't take long using a simple process of trial and error to estimate that if our transfer time is 28.623 days, then we can reduce our approach speed to a trivial 0.06 m/s. And with further tinkering we could probably reduce it even further.

In this particular case, the Lambert Solver transfer trajectory solution is:

What is happening here, is that because we have sought the minimum approach speed to the target Lagrange point orbit, we have indirectly found the stable manifold trajectory that connects our Start Point to the target Lagrange point orbit. Once we are on it, we will glide smoothy towards the reference orbit. And, once on it, we should need only minor station-keeping corrections to keep us n the target orbit.

Extensions
Although the above has focused on a very simple Lagrange point orbit, we can apply exactly the same kind of technique to approach other, more complicated, Lagrange point orbits -e.g., halo orbits. And although the above analysis has focused on the simple linear approximation, incorporating the various non-linear corrections to calculating target and transfer orbits is a just a matter or incorporating the appropriate detail and complexity of the full elliptical three-body model.

 

[ADSWNJ]

It's such a shame that Dr Keith Gelling is no longer active in this forum. He would love to have discussed this topic.

In terms of the mechanics of coding such a solution, you basically have two choices: code the Lambert Solver as a free-standing C++ executable, or code it as an MFD. In both cases, you need the solver to run in multi-threaded code, to unhook the deep calculation engine from the user interface. This technique is used in LagrangeMFD, so you can dial in a target calc time for its integrator - e.g. 2 seconds - without blocking the responsiveness of the Orbiter simulator - e.g. refreshing each 0.1 sec.

Interested to see how you get on with this project!

 

[cristiapi]

At the moment, the modelling and calculations are done in Mathematica. Although these can be ported to c/c++ without much effort, [...]

My biggest effort, at the moment, is for the calculation of the Lagrange points.

After the thread: https://www.orbiter-forum.com/showthread.php?t=38788 I found the NAIF files for the Lagrange points of the Sun-Earth system: https://naif.jpl.nasa.gov/pub/naif/generic_kernels/spk/lagrange_point/
If I do all the calculations told me by perseus in that thread, I get very different results from NAIF.
Is there any “standard” way to obtain an accurate position of L1, L2 and L3 for any system?

 

[MontBlanc2012]

Is there any “standard” way to obtain an accurate position of L1, L2 and L3 for any system?

Not sure I followed perseus' comments fully, but the link referenced by dgatsoulis made a lot more sense to me. At the risk of a little plagiarism, the core of the calculation of the position of L1, say, is to solve the equation:

[MATH]\alpha=\frac{\mu_{1}}{\left(\alpha+\mu_{2}\right){}^{2}}-\frac{\mu_{2}}{\left(\alpha-\mu_{1}\right){}^{2}}[/MATH]
for [MATH]\alpha[/MATH] where:

[MATH]\mu_1 = \frac{m_1}{m_1+m_2}[/MATH]
and

[MATH]\mu_2 = \frac{m_2}{m_1+m_2}[/MATH].

and where [MATH]m_1[/MATH] is the mass of the first gravitating body; and where [MATH]m_2[/MATH] is the mass of the first gravitating body.

This expression is quite general and works for all restricted three-body systems. The rest of the note referred to by dgatsoulis uses this value to calculate the position of the L1 point in Orbiter's inertial reference frame.

First, you have to calculate:

[MATH]\boldsymbol {e} = \frac {v^{2}} {gm}\, \boldsymbol {r} - \frac{\boldsymbol {r}.\boldsymbol {v}} {gm}\, \boldsymbol {v} - \frac{\boldsymbol {r}} {r}[/MATH]
[MATH]a = \frac {gm} {2\, gm/r - v^{2}}[/MATH]
[MATH]\nu = \arccos (\frac {\boldsymbol {e}.\boldsymbol {r}} {e\, r})[/MATH]
The first of these is the eccentricity vector (and the length of which is the orbital eccentricity, [MATH]e[/MATH]); the second is the set-major axis; and the third is the true anomaly of the elliptical orbit of the second gravitating body around the first gravitating body. In turn these quantities depend on the state vectors of the primary and secondary gravitating bodies - but you can always get this information by using an appropriate ephemeris calculator.

If you work through calculations in dgatsoulis' referenced thread, it seems to me that you should be able to calculate the position of L1 (or L2 and L3) in the context of the elliptical restricted three-body problem with some precision.

 

[cristiapi]

I tried to use the Keithth’s method, but it seemed overly complex to me and not good for the heavily perturbed systems (like Pluto and Jupiter).
But now I coded the Keithth’s Lua snippet for the SPICE library; the result is good for the Sun-Earth system (the ERTBP L1 point is no more than 11000 km away from the L1 point obtained from the SPICE file). But what about ERTBP method in other systems?

I thought to use the gravitational-centripetal accelerations method explained in the other thread because it seems an “always good” method, but it also seems that I don’t understand what happens in a Lagrange point, because the error in my calculations is too big.
This is a very simple method because I just need to find a point where centrifugal acceleration = gravitational acceleration, is that correct?

 

[MontBlanc2012]

Cristiapi

Strictly speaking, Lagrange points are only well-defined 'points' in the CR3BP and the ER3BP - i.e., only in specific classes of restricted three-body mathematical models. In a general n-body system, the Lagrange points do not (strictly speaking) exist since you will never get a stable point where gravitational and centripetal forces exactly cancel.

However, if perturbations are weak enough we can still make use of notional Lagrange points. Take, for example, the Earth-Moon system which is quite heavily perturbed by solar tidal forces. Tidal forces fall of as 1/r^3 rather than the usual 1/r^2 - so, at the Earth-Moon system, solar tidal forces are quite weak but still strong enough to distort noticeably the Moon's orbit around the Earth. Here, it seems to me that mission designers use three-body concepts from the ER3BP to construct 'reference orbits'. In principle, if the Earth-Moon system were a true restricted three-body system then tiny amounts of dV would be required to remain 'on station' on the reference orbit. In practice, mission designers have to allow for the fact that the perturbations mean that the satellite has to do more work to stay on the reference orbit to counter the effect of the perturbations and so more dV is used. As the magnitude of the perturbations increases, the satellite has to do more and more additional work to stay on the reference orbit. And, at some point if perturbations are high enough, the concept of a Lagrange point and a reference point breaks down since the dV cost of staying on it is simply too high. At this point mission designers need to go back to the drawing board and work, I presume, with more general n-body calculations.

wrt the position of L1 for the Sun-Earth system, are you using the state vectors of the Earth for your Lagrange point calculations or the state vectors of the centre of mass of the Earth and the Moon? Because the Sun-Earth L1 point is 1.5 million km from Earth, it is better to model the combined gravitational contribution of the Earth and the Moon as if it were emanating from a point source located at the centre of mass of the Earth-Moon system. Including the mass of the Moon in the calculations increases the mass of the secondary by about 1% and so will change the location of the L1 point by around 1% - i.e., by around 15,000 km. My guess is that you haven't included the Moon in your calculations and that doing so may well eliminate the bulk of the gap between your calculations and the spice file locations.

You mentioned 'the heavily perturbed systems of Jupiter and Pluto'. In what respect do you regard these as heavily perturbed systems?

 

[cristiapi]

My guess is that you haven't included the Moon in your calculations and that doing so may well eliminate the bulk of the gap between your calculations and the spice file locations.

You're right; if I use EMB instead of just Earth I get exactly the same JPL result (the difference is on the order of 10^-4 m).

You mentioned 'the heavily perturbed systems of Jupiter and Pluto'. In what respect do you regard these as heavily perturbed systems?

I mean that if I calculate Pluto-Charon L2 with ER3BP, the true equilibrium point will be very distant because of all the other satellites perturbations. The same for Pluto-Hydra and for the other planetary systems.
I did many simulations for the stability of L1 and L2, but it was 1 year ago...
Anyway, the result was as you say:

In a general n-body system, the Lagrange points do not (strictly speaking) exist since you will never get a stable point where gravitational and centripetal forces exactly cancel.

 

[MontBlanc2012]

For Pluto-Charon system, as you say, L2 is indeed heavily perturbed. Basically, trying to calculate the location of L2 is fairly meaningless here. L1 for Pluto-Charon ought to behave a little better, though. But still, because the mass of the other moons aren't negligible, a satellite parked in the vicinity of the notional ER3BP L1 point is still going to get kicked around a lot.

wrt the Lagrange points of the moons of Jupiter - Io, Europa, Ganymede, and Callisto. Yes, perturbations (due in large measure to Jupiter's J2 field) are likely to force higher station-keeping costs than one might ideally like. Callisto, in particular, ought to be reasonably well behaved. It's the furthest out so the J2 perturbation should be minimal. And if one replaces the state vectors of Jupiter with the state vectors of the centre of mass of Jupiter + Io + Europa +Ganymede, that also ought to improve things a little more.

wrt "if I use EMB instead of just Earth I get exactly the same JPL result (the difference is on the order of 10-4 m)", that's probably good enough