Optimal Earth-Venus-Mercury transfer - PyKEP / PyGMO

Keithth G

New member
Joined
Nov 20, 2014
Messages
272
Reaction score
0
Points
0
This is just a quick note on the use of a PyKEP / PyGMO to the design of an optimal ballistic transfer from Earth to Mercury with one gravity assist at Venus thrown in as well.


PyKEP / PyGMO
PyKEP and PyGMO is a trajectory planning modelling suite built by the Advanced Concepts Team of the European Space Agency. These are powerful, open-sourced (and free!) tools that are available for download onto your computer and for those interested, can be used to supplement the capabilities of TransX and IMFD enabling the design and implementation within Orbiter of some complex, real-world interplanetary trajectories.

OK, PyKEP and PyGMO don't operate within Orbiter as MFDs, they are more difficult to learn how to use and operate, and you will probably have to learn how to do some simple coding in Python to make them do useful things - but for those with an interest in mission design, and coupled with Orbiter, they offer broad scope for designing efficient, very realistic and complex missions.

In this post, I'm just going to work my way through an example by choosing something that is reasonably simple to do in PyKEP / PyGMO but quite hard to do with TransX / IMFD.


The sample problem - Earth - Venus - Mercury transfer
In this note, I'm going to focus on is an Earth - Venus - Mercury ballistic transfer. Although it isn't hard to set up this transfer in Orbiter using TransX, it is quite hard to find a transfer solution that minimise the fuel costs (Earth escape and Mercury orbit insertion). This is where PyKEP / PyGMO comes in: it turns out that this problem can be solved very efficiently and quickly using these tools and once a solution is found, one can have a high degree of confidence that it is the best possible solution.

To be a bit more specific, let's suppose that we want to find the lowest delta-V cost transfer from Earth - Venus - Mercury in, say, a 5,000 day windows from 1 Jan 2000 to 9 Sep 2013. Let's also suppose that, aside from the initial Earth escape burn and the Mercury orbit insertion burn, the trajectory is almost entirely ballistic - aside, that is, from a few minor mid-course corrections to keep the vessel 'on track'. Finally, lets suppose that you want to craft to leave Earth from a 300 x 300 km circular orbit; and end up at Mercury in a 300 x 300 km circular orbit.

This problem can easily be set up in PyKEP / PyGMO as an optimisation problem but, before doing that, let's simply state what the optimal trajectory is calculated to be:

Code:
Date of Earth departue:   2009-Jun-30 10:02:51.202508
Date of Venus encounter:  2009-Nov-25 10:43:25.318647
Date of Mercury arrival:  2010-Mar-29 00:44:14.568205

Transfer time from Earth to Venus:   148.03  days
Transfer time from Venus to Mercury: 123.58  days
Total mission duration:              271.61  days


TransX escape plan - Earth escape
--------------------------------------
MJD:                 55012.4186 
Prograde:             -7158.004 m/s
Outward:               -366.330 m/s
Plane:                  962.705 m/s
Hyp. excess velocity:  7231.737 m/s
Earth escape burn:     5376.649 m/s


Venus encounter
--------------------------------------
MJD:                 55160.4468 
Approach velocity:    15426.835 m/s
Departure velocity:   15426.835 m/s
Turning angle:           18.702 deg
Periapsis altitude:     984.058 km 
dV needed:                0.000 m/s


Mercury arrival
--------------------------------------
MJD:                 55284.0307    
Hyp. excess velocity:  6510.404 m/s
Orbit insertion burn   4810.736 m/s


Total fuel cost:      10187.385 m/s

It took PyKEP / PyGMO less than 1 minute to come up with this solution and I believe it to be optimal for this time period. As you can see from the bottom line, the total dealt-V fuel cost for this mission is 10,187 m/s. This is split into 5,376 m/s for Earth escape; and 4,810 m/s for Mercury orbit insertion. If, over the time period in question, any one can do better than this, I would be interested in hearing about it.

And this is a 'top-down' view of the optimal trajectory:



(The orbits of Earth, Venus and Mercury are shown in light blue. The optimal trajectory is shown in dark blue. The scale on the axes is Astronomical Units (AU). Note that arrival at Mercury is close to Mercury's perihelion.)

How does this compare with a Hohmann transfer? Over the same time period, the optimal PyKEP / PyGMO solution for a direct Earth - Mercury transfer has a delta-a fuel cost of around 12,800 m/s. So, the transfer via Venus (and the ensuing gravity de-assist) reduces mission costs by around 2,600 m/s. A significant fuel saving.

And is the very cheapest way of getting to Mercury? No, absolutely not. As per http://www.orbiter-forum.com/showthread.php?t=36584 one can throw in a few VILT (V-infinity Leveraged Transfers) sequences at the start and at the end to further reduce delta-V requirements by something around a further 4,000 m/s to 5,000 m/s. However, getting PyKEP / PyGMO to optimise the design of that trajectory is well beyond the scope of this note.


TransX escape plan
For convenience, I've converted the PyKEP / PyGMO Earth escape into numbers that TransX can use:

Code:
TransX escape plan - Earth escape
--------------------------------------
MJD:                 55012.4186 
Prograde:             -7158.004 m/s
Outward:               -366.330 m/s
Plane:                  962.705 m/s

Earth escape burn:     5376.649 m/s

For those familiar with TransX, the first for numbers are those required to set up the Earth escape plan on the "Escape Plan" page. If one does so, one should end up with something that looks like this:



with the encounter view looking something like this:



Th point to note here is that this Earth-Venus transfer was not worked out in Orbiter. It was worked out in PyKEP / PyGMO using an optimisation routine, then transcribed to TransX's coordinate system, and finally copied into TransX. And yet, even without any further adjustment, the closest approach to Venus is reported as a mere 13,000 km. What this says is that although I know that the underlying ephemerides and co-ordinate systems are not the same as Orbiter's, they are clearly very close cousins and, at the level of setting up manoeuvres within TransX, they can more or less be used interchangeably.

Once the transfer to Venus is nearly complete, one can then set up the Venus sling-shot and the sling-direct screen should look something like this:



OK, all well and good, but how does one use PyKEP / PyGMO?


The PyKEP / PyGMO code
To use PyKEP and PyGMO you have to have Python installed on your computer. You also need to separately install PyKEP and PyGMO. Once you have done so, then one opens Python and executes the following script:

Code:
from PyGMO.problem import base
from PyKEP         import planet, epoch, ic2par, fb_vel, lambert_problem, DAY2SEC, AU, MU_SUN
from PyKEP.planet  import jpl_lp
from math          import *
from numpy         import *


earth   = jpl_lp('earth')
venus   = jpl_lp('venus')
mercury = jpl_lp('mercury')

muE   = earth.mu_self
muV   = venus.mu_self
muM   = mercury.mu_self

radE  = earth.radius + 300000
radV  = venus.radius
radM  = mercury.radius + 300000


class my_problem_min(base):
    """
    Analytical function to minimize - a direct, Hohmann-like transfer from Earth to Mercury. 
    Transfer to be completed between 1 Jan, 2000 and 9 Sep, 2013 - a span of 5,000 days.  
    The goal is to minimise the total dV cost of escape from a 300 x 300 orbit around Earth
    to a 300 x 300 orbit around Mercury.

    USAGE: my_problem_min()
    """

    def __init__(self):
        super(my_problem_min, self).__init__(3)  # a problem in three variables
        self.set_bounds(0.001, 0.999)            # place bounds on those variables


    # Defines the objective function
    def _objfun_impl(self, x):
        t1     = x[0] *  5000             # departure time from Earth
        t2     = x[1] * (5000 - t1) + t1  # arrival tme at Venus
        t3     = x[2] * (5000 - t2) + t2  # arrival time at Mercury
        
        # calculate the state vectors of Earth, Venus & Mercury using JPL's Low Precision
        # ephemeris
        r1, v1 = earth.eph(epoch(t1))
        r2, v2 = venus.eph(epoch(t2))
        r3, v3 = mercury.eph(epoch(t3))
        
        # calculate the solutions of the two Lambert transfers
        l1     = lambert_problem(r1, r2, (t2 - t1)*DAY2SEC, MU_SUN)
        l2     = lambert_problem(r2, r3, (t3 - t2)*DAY2SEC, MU_SUN)
        
        # perform the dV calculations
        dv1    = array(l1.get_v1()[0]) - v1

        vin    = array(l1.get_v2()[0]) - v2
        vout   = array(l2.get_v1()[0]) - v2
        
        dv3    = array(l2.get_v2()[0]) - v3
        
        # sum the total fuel costs
        f      = sqrt(dot(dv1, dv1) + 2 * muE / radE) - sqrt(muE / radE )
        f      = fb_vel(vin, vout, venus) + f
        f      = sqrt(dot(dv3, dv3) + 2 * muM / radM) - sqrt(muM / radM ) + f
        
        # returb the total fuel cost
        return (f, )

    # Reimplement the virtual method that compares fitnesses - a minimisation problem
    def _compare_fitness_impl(self, f1, f2):
        return f1[0] < f2[0]

    # Add some output to __repr__
    def human_readable_extra(self):
        return "\n\tMinimisation problem - Earth-Venus-Mercury transfer"

This code sets up the optimisation problem in Python. It makes calls to the Kepler toolbox, PyKEP to calculate the mission fuel cost. Once this problem description has been set up, we use PyGMO to run the optimisation using a Differential Evolution algorithm, Here is the code:

Code:
prob = my_problem_min()
algo = algorithm.de(gen = 1000)

l = list()
for i in range(300):
    archi = archipelago(algo,prob, 8, 10, topology = topology.ring())
    for j in range(5):
        archi.evolve(5)
        print min([isl.population.champion.f[0] for isl in archi])
    tmp = [isl for isl in archi]; tmp.sort(key = lambda x: x.population.champion.f[0]);
    l.append(tmp[0].population.champion)
print "Results: \n"
print [ch.f[0] for ch in l]

and to print out the solution, one runs the following:
Code:
l.sort(key = lambda x: x.f[0])

y = array(l[0].x)

# calculate the times
t1     = y[0] *  5000             # departure time from Earth
t2     = y[1] * (5000 - t1) + t1  # arrival tme at Venus
t3     = y[2] * (5000 - t2) + t2  # arrival time at Mercury

# use the JPL low precision ephemeris to calculate the position of Earth, Venus and Mercury
r1, v1 = earth.eph(epoch(t1))
r2, v2 = venus.eph(epoch(t2))
r3, v3 = mercury.eph(epoch(t3))
l1     = lambert_problem(r1, r2, (t2 - t1)*DAY2SEC, MU_SUN)
l2     = lambert_problem(r2, r3, (t3 - t2)*DAY2SEC, MU_SUN)
        
# perform the dV calculations for Earth dearture, Venus encounter and Mercury arrival
dv1    = array(l1.get_v1()[0]) - v1
vin    = array(l1.get_v2()[0]) - v2
vout   = array(l2.get_v1()[0]) - v2        
dv3    = array(l2.get_v2()[0]) - v3
        
# sum the total fuel costs - f2 is the dV required to make the Venus flyby work
f1      = sqrt(dot(dv1, dv1) + 2 * muE / radE) - sqrt(muE / radE )
f2      = fb_vel(vin, vout, venus)
f3      = sqrt(dot(dv3, dv3) + 2 * muM / radM) - sqrt(muM / radM )

# calculate the TransX cordinates
fward   = v1 / linalg.norm(v1)
plane   = cross(v1, r1)
plane   = plane / linalg.norm(plane)
oward   = cross(plane, fward)

# print out the results
print "Date of Earth departue:  ", epoch(t1)
print "Date of Venus encounter: ", epoch(t2)
print "Date of Mercury arrival: ", epoch(t3)
print
print "Transfer time from Earth to Venus:  ", round(t2 - t1, 2), " days"
print "Transfer time from Venus to Mercury:", round(t3 - t2, 2), " days"
print "Total mission duration:             ", round(t3 - t1, 2), " days"
print
print
print "TransX escape plan - Earth escape"
print "--------------------------------------"
print("MJD:                 %10.4f "    % round(epoch(t1).mjd, 4))
print("Prograde:            %10.3f m/s" % round(dot(fward, dv1), 3))
print("Outward:             %10.3f m/s" % round(dot(oward, dv1), 3))
print("Plane:               %10.3f m/s" % round(dot(plane, dv1), 3))
print("Hyp. excess velocity:%10.3f m/s" % round(sqrt(dot(dv1, dv1)), 3))
print("Earth escape burn:   %10.3f m/s" % round(f1, 3))
print
print 
print "Venus encounter"
print "--------------------------------------"
print("MJD:                 %10.4f "    % round(epoch(t2).mjd, 4))
print("Approach velocity:   %10.3f m/s" % round(sqrt(dot(vin,vin)), 3))
print("Departure velocity:  %10.3f m/s" % round(sqrt(dot(vout,vout)), 3))  
ta  = acos(dot(vin, vout)/sqrt(dot(vin,vin))/sqrt(dot(vout,vout)))
print("Turning angle:       %10.3f deg" % round(ta*180/pi, 3))
alt = (muV/dot(vin,vin)*(1/sin(ta/2)-1) - venus.radius)/1000
print("Periapsis altitude:  %10.3f km " % round(alt, 3))
print("dV needed:           %10.3f m/s" % round(f2, 3))
print
print 
print "Mercury arrival"
print "--------------------------------------"
print("MJD:                 %10.4f    " % round(epoch(t3).mjd, 4))
print("Hyp. excess velocity:%10.3f m/s" % round(sqrt(dot(dv3, dv3)), 3))
print("Orbit insertion burn %10.3f m/s" % round(f3, 3))
print
print
print("Total fuel cost:     %10.3f m/s" % round(f1 + f2 + f3, 3))

This provides a basic template for designing a range of other missions in PyKEP / PyGMO.
 
Last edited:

boogabooga

Bug Crusher
Joined
Apr 16, 2011
Messages
2,999
Reaction score
1
Points
0
The following is an online database of optimum trajectories from 2010-2040:
http://trajbrowser.arc.nasa.gov/traj_browser.php

How do the PyKEP / PyGMO and online solutions compare for a 2017 launch?

I notice that the online optimum solution uses a powered rather than ballistic Venus flyby.
 

Keithth G

New member
Joined
Nov 20, 2014
Messages
272
Reaction score
0
Points
0
Hi, 'boogabooga'

I assume you are talking about this?

 
Last edited:

Hlynkacg

Aspiring rocket scientist
Addon Developer
Tutorial Publisher
Donator
Joined
Dec 27, 2010
Messages
1,870
Reaction score
3
Points
0
Location
San Diego
This looks really cool.

Edit: After some googling; here do I get this thing? and how do I install/use it?
 
Last edited:

Keithth G

New member
Joined
Nov 20, 2014
Messages
272
Reaction score
0
Points
0
The following is an online database of optimum trajectories from 2010-2040:
http://trajbrowser.arc.nasa.gov/traj_browser.php

How do the PyKEP / PyGMO and online solutions compare for a 2017 launch?

I notice that the online optimum solution uses a powered rather than ballistic Venus flyby.

Ok, assuming that the database trajectory that you have in mind is:



here are a few comments in response:

First, we have to be careful to make sure that we are comparing like-for-like. There are a number of things that we need to focus on:

1. The radius of the departure parking orbit around Earth

2. Powered flybys around Venus

3. Single vs multi-revolution solutions to the Lambert problem

4. Radius of the parking orbit around Venus

Let's look at each of these in turn. The radius of the departure parking orbit around Earth affects the magnitude of the escape burn needed to achieve the right Venus-bound departure trajectory. Working backwards from the C3 value quoted in the database entry, I conclude that the altitude of the low Earth parking orbit is somewhere around 200 - 300 km. The level of precision quoted in the database entry isn't sufficient to be more precise than this but, overall, I would think it reasonable that there is little material difference between the value of 300 km that I have used in my Python optimisation script and that assumed in the database entry. So, we can reasonably assume that on this measure, at least, we are compare like-for-like.

As you have noted, the database's optimal trajectory to Mercury for 2017 requires a powered flyby of Venus. As it turns out, my script also takes into account of an possible delta-V needed to effect the flyby. In some cases, this could be 0 m/s, but the possibility exists nonetheless for a powered flyby. In this respect, too, the database entry and the Python script are alike.

The database optimal trajectory solution clearly focuses on a solution for which, on the Venus-Mercury leg, the spacecraft completes at least one revolution around the Sun. In my Python script, I make use only of the single-revolution Lambert solution. I have endeavoured to augment my script to include must-revolution solutions but, as yet, I am not convinced that it is working properly.

Finally, the database entry states that the delta-V required for Mercury orbit insertion is the 'C3 = 0' value. This means that it is the fuel needed to achieve capture by Mercury, but not orbit circularisation. My original Python script assumed that the goal was to achieve a low, 300 km circular orbit around Mercury. I have updated my script to allow for Mercury capture (C3 = 0) where I have assumed that capture is achieved by executing a retrograde burn at Mercury periapsis which I set to a value of 300 km.

By making these modifications, I believe that my (updated) Python script and the database entry can be evaluated on a level playing field. Running, then, my Python script through PyKEP / PyGMO for a 5,000 day period that covers 2013 to 2026 (give or take). The optimal solution found by PyKEP / PyGMO is as follows:

Code:
Date of Earth departue:   2017-Jan-03 03:20:49.746921
Date of Venus encounter:  2017-Apr-07 10:26:54.186238
Date of Mercury arrival:  2017-May-26 00:00:23.006635

Transfer time from Earth to Venus:   94.3  days
Transfer time from Venus to Mercury: 48.56  days
Total mission duration:              142.86  days


TransX escape plan - Earth escape
--------------------------------------
MJD:                 57756.1395 
Prograde:             -4097.239 m/s
Outward:                559.155 m/s
Plane:                -1237.272 m/s
Hyp. excess velocity:  4316.349 m/s
Earth escape burn:     4021.844 m/s


Venus encounter
--------------------------------------
MJD:                 57850.4353 
Approach velocity:     8412.598 m/s
Departure velocity:    8683.800 m/s
Turning angle:           50.495 deg
Periapsis altitude:     119.678 km 
dV needed:              439.877 m/s


Mercury arrival
--------------------------------------
MJD:                 57899.0003    
Hyp. excess velocity:  7747.429 m/s
Orbit insertion burn   4713.578 m/s - C3 = 0


Total fuel cost:       9175.298 m/s

This solution is clearly different from the database entry - being some 300 m/s cheaper - but it also requires a powered Venus flyby.

Does this answer your question?

---------- Post added 02-28-16 at 01:55 AM ---------- Previous post was 02-27-16 at 11:07 AM ----------

This looks really cool.

Edit: After some googling; here do I get this thing? and how do I install/use it?

Installation of PyKEP / PyGMO takes a bit of effort, I'm afraid. Although it is a well-written, powerful piece of software provided free of charge by the European Space Agency, it would appear that that venerable institution does not subscribe to the 'Plug-n-Play' concept of software installation. Rather the ESA, appears to subscribe to the 'Let's-make-it-really-difficult-to-install' concept instead.

The first thing to do is install Python on your computer, if it is not already installed. (On Macs, Python is usually pre-installed, for example). Here is a link to the Python download page:
https://www.python.org/downloads/
Python is like Orbiter's Lua scripting language: it is the scripting language that is used to fire-up PyKEP and PyGMO and without it, and unless you are a good C++ programmer, the functionality of PyKEP and PyGMO remains inaccessible.

Next, you will need to install PyGMO. PyGMO is the optimisation engine used by PyKEP. The download page for PyGMO is here: http://esa.github.io/pygmo/install.html
Before downloading PyGMO, you will need to install a number of dependencies: 'git', 'CMake' and the 'Boost' libraries. Links to those are given on the PyGMO installation page. You will need to follow the download and installation instructions carefully. Part of the installation process involves compiling the PyGMO source code on your computer, so you will have to ensure that your operating system has access to a suitable C++ compiler. The compilation process is automated via a 'makefile' script - part of the PyGMO download, but it still needs access to a compiler.

Finally, you will need to install PyKEP. PyKEP is the Kepler toolbox. The download page is here:
https://esa.github.io/pykep/downloading.html
and the installation page is here:
https://esa.github.io/pykep/installation.html
Again, installation of PyKEP can be a somewhat frustrating experience. But persevere.

Once, all of this has been installed you can then run PyKEP / PyGMO. To do this, you need to get Python running and then execute a Python script (see previous posts) that call PyKEP and PyGMO. To do this, it may make sense to install Jupyter notebook interface for Python - which provides an altogether better interface for using Python. The download page for Jupyter is:
http://jupyter.readthedocs.org/en/latest/install.html

The whole download and installation experience is not straightforward, but if you manage to navigate your way through it, you will end up a powerful piece of trajectory planning and optimisation software that is largely compatible with Orbiter.

Good luck!
 
Last edited:

boogabooga

Bug Crusher
Joined
Apr 16, 2011
Messages
2,999
Reaction score
1
Points
0
The radius of the departure parking orbit around Earth affects the magnitude of the escape burn needed to achieve the right Venus-bound departure trajectory. Working backwards from the C3 value quoted in the database entry, I conclude that the altitude of the low Earth parking orbit is somewhere around 200 - 300 km. The level of precision quoted in the database entry isn't sufficient to be more precise than this...

FYI, according to the documentation, the altitude of the earth parking orbit is 200 km.

I find it intriguing that you seem to have found a local optimum solution that is apparently not in the database.

Still, for easy installation I recommend GMAT. GMAT should have a lot of the same capabilities. However, optimization may require access to MATLAB, which is not free, but many people with engineering inclinations keep a copy on hand anyway.
 
Last edited:

Keithth G

New member
Joined
Nov 20, 2014
Messages
272
Reaction score
0
Points
0
Still, for easy installation I recommend GMAT. GMAT should have a lot of the same capabilities. However, optimization may require access to MATLAB, which is not free, but many people with engineering inclinations keep a copy on hand anyway.

I have tried twice to install GMAT but failed on both occasions. I'm a Mac user and have a MacBook Pro running the latest OSX El Capitan operating system. One of the GMAT dependencies is 'wxwidgets'. The installer for 'wxwidgets' on Mac OSX El Capitan appears to be broken and I am waiting for this to be 'fixed' before I can complete the GMAT installation process.

Sadly, I am not a MATLAB user so, for me, trajectory optimisation possibilities using GMAT seem a little limited.
 
Last edited:
Top