Accuracy test of the Orbiter integration engine

[Keithth G]

Hi, all

I've been conducting a number of simple tests of the accuracy of Orbiter's underlying numerical integration scheme and I thought I would share the results.

I've written a small program that calculates the solution of the parabolic Lambert problem with essentially arbitrary precision. The parabolic Lambert problem can be stated as: calculate both the initial velocity vector at and time of flight needed to travel along a ballistic parabolic trajectory from one point, A, to another point, B, under the influence of a single Newtonian gravitating body (in this case, the Sun).

The 'test', then, is to enter manually the 'firing solution' into Orbiter and then 'time-warp' forward to see how close to the destination point, B, Orbiter's numerical integration of the ballistic trajectory goes. If Orbiter's numerical integration is typically 'close' to the destination point, B, it means two things: first, n-body planetary perturbations for parabolic motion are, for the most part, not particularly significant; and second, Orbiter's underlying numerical integration engine is accurate - i.e., it isn't really a random number generator in disguise. Furthermore, if one gets essentially the same results for different 'time-warp' values (e.g., 1,000x, 10,000x and 100,000x), then the integrator is probably robust to changes in 'time-warp.

To cut a long story short, Orbiter's numerical integration engine seems to be surprisingly accurate. At 10,000x and 100,000x, the integrated ballistic trajectory seems to be essentially identical (to within a few km over flight path distances of around 5 AU).

As for the accuracy of the 2-body parabolic ballistic trajectory, this too seems to provide a very good 'first guess' of the (parabolic) ballistic trajectory between two points in the solar system's n-body gravitational field. As an example, I placed a ship close to Earth's orbit (but not close to Earth) and calculated the ballistic trajectory to Olympus Base on the surface of Mars. The parabolic trajectory was directed 'the long way' around the Sun, with a solar periapsis of 0.765 AU. Duration of flight was 111.4 days. At a time-warp of 10,000x, and following the ballistic trajectory (without any form of mid course correction, of course), Orbiter then had the ship pass overhead Olympus Base at an altitude of 300km. The Mars encounter velocity was around 60+ km/sec which is the sort of speed one would expect for a comet-like, parabolic trajectory.

(It is likely that that the the residual 'error' of 300km is due to n-body perturbations, but I will need to work a little harder to find an accurate 'n-body' firing solution rather that the simple analytical, '2-body' solution.)

For those interested, the analytical 2-body parabolic trajectory firing solution for this example is as follows:

Firing date: MJD 51870.619542570596

Firing co-ordinates:
X 149597870700.0 m
Y 0.0 m
Z 14959787070.0 m

dX/dt -16773.33218530575 m/s
dY/dt 595.872620354075 m/s
dZ/dt -38519.48368154611 m/s

If one wishes to enter this 'firing solution' into Orbiter, the procedure for doing this is as follows:
1. Have Orbiter 'focus' on an undocked Delta Glider
2. 'Pause' Orbiter by hitting Ctrl-P ttemprarily suspend the Orbiter integration engine. While paused make a series of changes using the scenario editor as follows:
a. Open the Scenario Editor and in the 'and amend the current date to the 'Firing Date' given above by using 'copy & paste'. Remember to 'apply' the date before moving on
b. Now use the Scenario Editor to update the state vectors of the ship as per the above. Before entering these changes, change the orbit reference to "Sun / Ecliptic / Cartesian". (This is on the right-hand side of the same panel for changing the state vectors)
c. Use scenario editor to 'Kill Rotation' for the ship by setting all of the angular velocity components to zero (This doesn't affect the trajectory, but simply stops the ship from making arbitrary rotations - a nuisance at high time-warp)

6. Hit 'Ctrl-P' again to start the Orbiter integration engine
7. Time-warp at, say, 10,000x or 100,000x. Open MFDs of your choice and watch the solar system go by until your reach Mars.

 

[Zatnikitelman]

Was this done under fixed timesteps? If not, what was your framerate from the Orbiter-side of things?

 

[Keithth G]

In truth, I don't know. But if you can can tell me where I can find these numbers in Orbiter, I will let you know

 

[Zatnikitelman]

Well, you'd know if you were running fixed timestep as you have to go deep into the launchpad to enable it. Framerate can be quickly displayed by pressing F and there's a framerate grapher somwhere, I'll have to look this weekend if you don't find it before then.

 

[Keithth G]

I'm running a version of Orbiter 2010 in its default configuration on a MacBook Pro within a Wine wrapper. I haven't made any changes to the configuration files that would affect the integration time step. The frame rate is 30.

 

[Keithth G]

Updated analysis

Just as an adjunct to the above, I've continued to d a few accuracy checks on Orbiter's numerical integration scheme. The headline news is that the after a 111 day flight (and circa 4 AU travelled ), Orbiter's integration engine is accurate to within about 5km - or roughly 1 in 10^8. This seems pretty reasonable to me.

The details of new tests:

The original test outlined in an earlier post was between the theoretical solution of a 2-body boundary value problem connecting two points in space, 'A' and 'B' with a high-speed parabolic trajectory, whereas Orbiter uses an n-body numerical integrator to integrate the trajectory from the starting point 'A' forward in time using the theoretical value of the starting velocity of the parabolic trajectory at the point 'A'. The test was to see how close Orbiter's integration engine's estimate of the trajectory came to the end-point 'B'. In principe, the theoretical solution acts as a benchmark against which to gauge the accuracy of Orbiter's integration engine. If Orbiter integrates the trajectory so that it is 'close' to the point 'B' then, the integration engine is working well. If not, then it isn't.

Earlier, I reported that for at least one parabolic trajectory ending at Olympus Base on Mars (with a trajectory duration of 111.4 days, with a specified starting point and velocity vector, and with a total flight distance of around 4 AU), Orbiter projected that the trajectory would pass within about 300 km of Olympus Base.

This seemed pretty reasonable. However...

This wasn't quite a like-for-like comparison. First, Orbiter's integration engine takes into account the perturbation of other planetary bodies - e.g., Jupiter and Saturn. The theoretical solution is a two-body solution and does not take into account the gravitational influence of any bodies in the Solar System other than the Sun. Second, Orbiter's integration engine also takes into account the motion of the Sun around the Solar System's barycentre (or Centre of Mass). Jupiter and Saturn are so massive that they influence significantly the position of the Sun, and the Sun rotates around the Solar System barycentre with a radius of a few hundred thousand km. Over the course of 110 days, the Sun moves appreciably and this motion is not taken into account the theoretical two-body dynamics.

Now, to make a like-for-like comparison one can either make the theoretical benchmark look more like the real Solar System - introducing a more complex system of n-body dynamics into the theoretical trajectory planner - or, one can make the Orbiter's representation of the Solar System look more like the theoretical 2-body system. The former is hard; the latter is easy - so, naturally, I opted for the latter. There are two adjustments to Orbiter's Solar System that need to be made to do this. First, in the planetary configuration files, the masses of the planets (and moons) were set to zero. This left the Sun as the only gravitating body in the Solar System - in perfect alignment with the theoretical 2-body benchmark. Second, the VSOP87 data file for the Sun was amended. In Orbiter, the Sun and the planets 'run on rails' defined by a set of coefficients used to enumerate the widely used VSOP87 ephemeris. This ephemeris engine is used to calculate where the planets (and the Sun) will be located at any given time. By amending those coefficients, one can trick Orbiter into thinking that it should place the Sun at the Solar System barycentre.

Having made these changes, one now has to re-work the theoretical firing solution. These changes give rise to a small change in the 2-body theoretical initial conditions as follows for arriving at Olympus Bas, Mars on MJD 51982.029235:

Firing date: MJD 51870.622998861545 (formerly 51870.619542570596)

Firing co-ordinates:
X 149597870700.0 m
Y 0.0 m
Z 14959787070.0 m

dX/dt -16773.770770805266 m/s (formerly -16773.33218530575 m/s )
dY/dt 595.0251466500241 m/s (formerly 595.872620354075 m/s )
dZ/dt -38519.30685783877 m/s (formerly -38519.48368154611 m/s )

Running these numbers though the amended Orbiter Solar System model then has the ship overfly Olympus Base at an altitude of 5 km (formerly 300 km). In other words, the bulk of the initial discrepancy between theory and Orbiter's integration of the parabolic trajectory has been explained and is now sufficiently small that there seems little benefit in attempting to resolve it further.

Personally, I am satisfied that Orbiter's engine is sufficiently accurate to support quite detailed mission planning.

P.S. The integration was carried out at 100,000x using 8th order symplectic Runge Kutta integration, without fixed time warp and with a frame rate of 30.

---------- Post added at 07:02 AM ---------- Previous post was at 05:16 AM ----------

And one last update on this. In the above, I adjusted the VSOP87 files so that the centre of the Sun was close to (~100 km away from) the Solar System barycentre. I would have set it to zero, but apparently Orbiter crashes if I try to do this.

The closest I can shift the Sun to the barycentre is about 25 km. With this minor update, Orbiter's integrated trajectory takes the ship within 1 km of Olympus Base (the target point 'B'). Or an accuracy of around 1 in 10^9. (It is a little hard to tell, though exactly how close the ship goes since its encounter velocity with Mars is around 60 km/sec - so it overflies the base in a fe hundredths of a second, which is not really enough time for MFDs to update.)

Anyway, this is good enough, methinks.

 

[boogabooga]

Orbiter is excellent at what it does, but see for example GMAT:
http://gmat.gsfc.nasa.gov/

has more options for fancier numerics such as Prince Dormand, time step control, higher order gravity potential models, etc.

 

[Keithth G]

Thanks.

 

[Keithth G]

OK - one final addition to this 'Orbiter accuracy test' thread.

The story so far: The motivation for this thread was an interest in understanding how well Orbiter 2010's integration engine respected Newtonian gravity. And the specific test that I wished to apply was a variant on the Lambert boundary value problem: if I am at a point A now and I wish to get to point B in time dt, what velocity vector do I need now to do this? For the standard 2-body gravitational problem, there are analytical results (i.e., formulas and simple algorithms) that one can use to answer this question. And since I was going through an exercise of working out those analytical results in any case, and the first analytical results to emerge were for parabolic orbits (rather than the more conventional elliptical orbits), I thought to myself "I can use this result to benchmark Orbiter's integration engine".

And so I did.

And It turned out that if one simplified Orbiter's representation of the Solar System so that the only gravitating body was in it was the Sun sitting at the centre of an inertial coordinate system then over a parabolic trajectory of around 4 AU lasting 110 days from a point near Earth's orbit to a point on Mars, it was difficult to identify any material difference between the analytical result and Orbiter's numerical integration of the trajectory. To be sure, there would have been differences, but the were too small for me to easily measure.

All well and good, of course, but what about the real system in which there are material perturbations to the idealised 2-body trajectories of ellipses, parabolas and hyperbolas? How well does Orbiter's integration engine perform then? (A good question.)

To answer this question requires putting perturbative forces back into Orbiter's representation of the Solar System and building a carefully calibrated and tested n-body numerical integrator - i.e., benchmarking one integration engine (Orbiter's) with another (mine). So, I reinstated the barycentric motion of the Sun together with Mars' gravitational field - thereby introducing two significant perturbations of 2-body motion.

For this test, I also changed the target point 'B' from the Olympus base on Mars to the centre of Phobos. As it turned out, the parabolic trajectory from the starting point 'A' to Phobos then required a fast close encounter with Mars at an altitude of 2163 km above the surface of Mars, before reaching (OK, 'crashing into') Phobos a few minutes later. This close encounter with Mars introduced a kind of 'sling shot' effect: as it passed, the ship received a momentum 'push' from Mars and slowed the ship by about 400 m/s and, moreover, introduced a change in direction of about 2 degrees (which, at the altitude of Phobos above Mars amounted to a trajectory shift a few hundred kilometres). And why is all this change in destination important? Well, for 99.999% of the flight the trajectory was almost entirely determined by the gravitational influence of the Sun. But for the last 0.001% of the trajectory before reaching Phobos - i.e. about 10 minutes of the flight, the ship encountered a sharp 'spike' in the gravitational field that altered its trajectory significantly. This close encounter with Mars rather serendipitously provided a good test for an integration engine since the engine has to model two processes occurring on vastly different timescales - the first being the sedate 'encounter' with the Sun as the ship made its parabolic way from point 'A' to Phobos; and the second being a very brief, near hyperbolic encounter with Mars.

And what was the integration engine that I used to benchmark Orbiter's? Well, I built three actually - all of which yielded same trajectory within a few metres. The first was an 10th order symplectic partitioned Runge Kutta integrator based on 'NDSolve' - a suite of algorithms native to Mathematica - using a working precision of in excess of 30 significant digits, well beyond that of normal floating-point arithmetic. The second was a very high order implicit symplectic discrete variational integrator based on Gauss-Lobatto quadrature - also with a working precision of 30 significant digits. And the third was a low (4th order) explicit, exactly time-reversible lattice symplectic integrator - with a lattice spacing of about 1 micron. With each of these integration engine, a newton iteration scheme was devised to identify the required velocity vector at the starting point 'A' to reach (OK, 'hit') Phobos on the target date. The Newton root-finding iteration scheme was applied until the trajectory target point was reached to a tolerance of less than 1 metre. The time step and working precision of the three integrators was varied until each method yielded a consistent trajectory on the < 1 metre level. The point 'A' starting velocity was compared across the three models. As it turned out, the three models yielded the same starting velocity to within 10^-6 m/s. Good enough - and well beyond the level of accuracy with which one could ever hope to measure velocity of a spacecraft on a real interplanetary mission.

So, what happens when these 'theoretical' starting points were put through Orbiter? Well, unfortunately the ship consistently arrived at the orbit of Phobos 1.7 seconds too late (on a journey lasting roughly 100 million seconds) and, consequently, striking Phobos 8.1 km off-centre. Not bad, but materially different from the 'theoretical' result. After having laboriously checked all model parameters - e.g. gravitational constants, planetary and solar masses, J2 fields and the like, I have run out of possible sources of explanation for the difference other than to attribute the difference to rounding error accumulation (or similar) in Orbiter's integration engine. To materially change Orbiter's trajectory forecast so that one strikes Phobos 'dead centre' would require a change in the initial velocity of around 10^-3 m/s - about three orders of magnitude more than the accuracy with which the reference integrators determined the initial velocity vector.

Not, of course, that this round-off accumulation (or whatever it is) is large enough to matter - or indeed be measurable in the real world. The fact of the matter is that Orbiter's integration engine does a very good job at simulating Newtonian gravity and achieves a level of accuracy far greater than realised by either of the commonly used MFD planning tools, IMFD and TransX. For all intents and purposes, other than perhaps planning a real space mission, Orbiter's integration engine does not appear to be materially in error.

 

[Zatnikitelman]

Was your framerate still around 30 when you got the imperfect results? Have you tried either running Orbiter_ng headless (really really really high fps) or using D3D9 which depending on the computer, should drastically improve the framerate?

 

[Keithth G]

I've played around with the size of Orbiter's time-steps to some extent. So far, it hasn't made an appreciable difference to the outcome.

Currently running D3D9.

I will try with Orbiter_ng (headless) to see if that makes a difference.

 

[Zatnikitelman]

Are you getting the same 30fps with D3D9?

 

[Keithth G]

Yes. But I'm running Orbiter 2010 on a MacBook Pro in a Wine wrapper. The Wine wrapper emulates a Windows machine so, inevitably, the frame rate is going to be slower than one would expect.

Just out of interest, though, I looked at the frame rate for Orbiter_ng running in server mode. The frame rate there was reported to be circa 3,000.

 

[martins]

Quite an interesting test. It would be nice to find out where the residual discrepancy is coming from. If you did try varying the time step parameter without significant change in the result, then I guess the source of the difference is probably not the accuracy of the integrator. I assume you were setting Orbiter to use its highest order available propagators (RK8 or SYM8)? Over what range did you vary the time step intervals?

A few other possible sources of error come to mind:

1. Orbiter doesn't compute the VSOP series solution at every time step for performance reasons. It computes it at regular intervals, and then interpolates for a given intermediate time point (it does a linear interpolation of the spherical coordinates). The time interval between computations of the VSOP solution is defined by the "SamplingInterval" parameter in a planet's config file. You could try to reduce this parameter for the relevant planets and the sun, and see if that makes a difference.
2. Orbiter doesn't sum the VSOP solution over all available terms, again for performance reasons. The log file will show the total number of terms, and the number used by Orbiter. The number of terms Orbiter uses is governed by the "ErrorLimit" parameter in the planet's config file. A value <= 1e-8 should ensure that all terms are included for a given body. You could try to reduce that value for the affected celestial bodies to see if it makes a difference.
3. Orbiter doesn't include all celestial bodies present in the simulation in the dynamic state update of a vessel, but only the most significant ones at the given location. If you display the vessel's info sheet (Ctrl-I) it should give a list of celestial bodies included in the current update (at least it does for the Orbiter beta. I can't remember for sure if this is already available in the 2010 release, but it probably is). Have a look at the list and check if it contains all the celestial bodies you would expect at all stages of the trajectory. Currently there is no way for the user to modify the threshold of inclusion/exclusion of gravitational bodies, but if this turns out to be a major limitation of accuracy, I could certainly add it.

Let us know if any of those measures has an effect on the outcome.

Edit:

Yes. But I'm running Orbiter 2010 on a MacBook Pro in a Wine wrapper. The Wine wrapper emulates a Windows machine so, inevitably, the frame rate is going to be slower than one would expect.

You should really make sure to activate the "fixed time steps" debugging feature for your tests, to make the results independent of frame rate and computing platform, otherwise it's too difficult to interpret the results.

 

[martins]

Another question: From your previous post, I gather that you have created a quite sophisticated numerical n-body propagator engine. Would you be happy to make it available to the Orbiter community? This would come in handy in various places, e.g. for mission planning tools, or for creating our own power series ephemeris solutions for celestial bodies (I'm still looking for a solution for Hyperion).

 

[Keithth G]

I quite agree that it would be nice find out the source of the discrepancy. Although an offset of a few seconds lasting 111 days may no sound much, it is, perhaps, larger than one would expect.

In response to your specific questions:

1. I have set the Orbiter's integrator to SYM8 and this integrator has been used consistently throughout these tests.

2. wrt to time-steps, I used three integrators. The first was Mathematica's native 'NDSolve' function. This allows setting of accuracy and precision targets (as well as the working precision. The accuracy and precision targets were set to 20 (digits, I presume), and the working precision was set to 30 (digits). Mathematica then did the rest to adjust time-steps to achieve these accuracy and precision targets. The accuracy and precision targets were set at these levels so as to ensure that the integrated trajectory had converged at the < 1m level. For the other two integrators, I used time steps of 8 minutes, 4 minutes, 2 minutes, 1 minutes and 30 seconds. Trajectory convergence at the < 1m level was reached with around the 2 minute level and confirmed by carrying out the tests again at the 1 minute time step level.

3. I was aware that Orbiter didn't automatically use the full VSOP87 solution: my first line of attack was to build the full VSOP87 planetary function calls using an on-line VSOP87 code generation tool, but spot checks confirmed against a series of test dates confirmed that for the Sun, say, there was a substantial discrepancy in Orbiter's reported positions from VSOP87 of circa 20km - 40km depending on the date. Rather than attempt to decompose the specific algorithm used by Orbiter to determine its planetary and solar positions, I decided that most expedient way forward was to use Orbiter as an ephemeris generator itself - not of VSOP87, per se, but of whatever planetary state vectors that Orbiter used itself (given the default values used in the configuration files). So, for the period of the anticipated trajectory, I built two sets of interpolation functions - one for the Sun and one for Mars (the only two gravitating bodies in my minimalist Solar System). These interpolation functions were based on 100 and 200 sampling points over the 111 day period of the trajectory by interrogating Orbiter for planetary position using a lua script. Two kinds of interpolation functions were built: one based on Mathematica's native 'Interpolation' function; and the other based on Chebyshev interpolation. Spot checks of planetary positions against those reported by Orbiter via lua scripting indciated that the error in planetary and Solar positions was of the order of a few metres. Using all four variants of the interpolation functions produced essentially the same trajectory when run through my integrators.

4. I hadn't realised that Orbiter was a bit selective about which gravitating bodies to use in its integration. Good to know. In my integration, I have assumed that all gravitating bodies contributed to the local gravitational field no matter the location of the ship. Having checked in Orbiter 2010, I see that at the start of my test trajectory, 'the Sun' is the only gravitating source - but at some point during the trajectory, this switches from only 'the Sun' to 'the Sun' and 'Mars'. At a guess, I would say that the difference between Orbiter's integrated trajectory and my own is largely due to the different gravity model used by Orbiter's integrators and mine.

5. wrt 'fixed time steps', noted. However, I did do a few tests withs fixed time steps and found that for the values I used, there was no appreciable change in Orbiter's integrated trajectory.

6. Yes, happy to make my numerical integrators available - and will do so in due course.

 

[Siliconaut]

I'm just going to take a moment here off topic for a moment so please forgive me. But as a guy with a 7th grade math education when I read these threads where a guy whips up 3 separate advanced math engines with offhand ease just to test the mathematic accuracy of orbiter just leaves me stunned. I love Orbiter, but 90% of the numbers on the MFD's are incomprehensible gibberish to me even after a decade or more of play. It's almost like magic how you guys understand all this stuff on such a deep level.

 

[martins]

Applied maths is fun! :thumbup:

 

[Lucy]

Applied maths is fun! :thumbup:

Orbiter's entire reason for being, if I'm not mistaken...

 

[Siliconaut]

Applied maths is fun! :thumbup:

Well Martin, your applied math as provided me with so many hundreds of hours of fun! More than any other application I have. In 2001 I found this applied math project and it has never left my hard drive since. I'll never forget how breathtaking it was for me when Atlantis achieved orbit for the first time and I looked out over the Earth horizion! Or the first time I came into the Jovian system watching the Galliean moons circle the great planet. Thank you!