Degenerate Conic

Earth-Mars Free Return

Let's try using the Fortran Astrodynamics Toolkit and Pikaia to solve a real-world orbital mechanics problem. In this case, computing the optimal Earth-Mars free return trajectory in the year 2018. This is a trajectory that departs Earth, and then with no subsequent maneuvers, flies by Mars and returns to Earth. The building blocks we need to do this are:

1. An ephemeris of the Earth and Mars [DE405]
2. A Lambert solver [1]
3. Some equations for hyperbolic orbits [given below]
4. An optimizer [Pikaia]

The \(\mathbf{v}_{\infty}\) vector of a hyperbolic planetary flyby (using patched conic assumptions) is simply the difference between the spacecraft's heliocentric velocity and the planet's velocity:

• \(\mathbf{v}_{\infty} = \mathbf{v}_{heliocentric} - \mathbf{v}_{planet}\)

The hyperbolic turning angle \(\Delta\) is the angle between \(\mathbf{v}_{\infty}^{-}\) (the \(\mathbf{v}_{\infty}\) before the flyby) and the \(\mathbf{v}_{\infty}^{+}\) (the \(\mathbf{v}_{\infty}\) vector after the flyby). The turning angle can be used to compute the flyby periapsis radius \(r_p\) using the equations:

• \(e = \frac{1}{\sin(\Delta/2)}\)
• \(r_p = \frac{\mu}{v_{\infty}^2}(e-1)\)

So, from the Lambert solver, we can obtain Earth to Mars, and Mars to Earth trajectories (and the heliocentric velocities at each end). These are used to compute the \(\mathbf{v}_{\infty}\) vectors at: Earth departure, Mars flyby (both incoming and outgoing vectors), and Earth return.

There are three optimization variables:

• The epoch of Earth departure [modified Julian date]
• The outbound flight time from Earth to Mars [days]
• The return flyby time from Mars to Earth [days]

There are also two constraints on the heliocentric velocity vectors before and after the flyby at Mars:

• The Mars flyby must be unpowered (i.e., no propulsive maneuver is performed). This means that the magnitude of the incoming (\(\mathbf{v}_{\infty}^-\)) and outgoing (\(\mathbf{v}_{\infty}^+\)) vectors must be equal.
• The periapsis flyby radius (\(r_p\)) must be greater than some minimum value (say, 160 km above the Martian surface).

For the objective function (the quantity that is being minimized), we will use the sum of the Earth departure and return \(\mathbf{v}_{\infty}\) vector magnitudes. To solve this problem using Pikaia, we need to create a fitness function, which is given below:

subroutine obj_func(me,x,f)

!Pikaia fitness function for an Earth-Mars free return

use fortran_astrodynamics_toolkit

implicit none

class(pikaia_class),intent(inout) :: me !pikaia class
real(wp),dimension(:),intent(in) :: x !optimization variable vector
real(wp),intent(out) :: f !fitness value

real(wp),dimension(6) :: rv_earth_0, rv_mars, rv_earth_f
real(wp),dimension(3) :: v1,v2,v3,v4,vinf_0,vinf_f,vinfm,vinfp
real(wp) :: jd_earth_0,jd_mars,jd_earth_f,&
dt_earth_out,dt_earth_return,&
rp_penalty,vinf_penalty,e,delta,rp,vinfmag

real(wp),parameter :: mu_sun = 132712440018.0_wp !sun grav param [km3/s2]
real(wp),parameter :: mu_mars = 42828.0_wp !mars grav param [km3/s2]
real(wp),parameter :: rp_mars_min = 3390.0_wp+160.0_wp !min flyby radius at mars [km]
real(wp),parameter :: vinf_weight = 1000.0_wp !weights for the
real(wp),parameter :: rp_weight = 10.0_wp ! penalty functions

!get times:
jd_earth_0 = mjd_to_jd(x(1)) !julian date of earth departure
dt_earth_out = x(2) !outbound flyby time [days]
dt_earth_return = x(3) !return flyby time [days]
jd_mars = jd_earth_0 + dt_earth_out !julian date of mars flyby
jd_earth_f = jd_mars + dt_earth_return !julian date of earth arrival

!get earth/mars locations (wrt Sun):
call get_state(jd_earth_0,3,11,rv_earth_0) !earth at departure
call get_state(jd_mars, 4,11,rv_mars) !mars at flyby
call get_state(jd_earth_f,3,11,rv_earth_f) !earth at arrival

!compute lambert maneuvers:
call lambert(rv_earth_0,rv_mars, dt_earth_out*day2sec, mu_sun,v1,v2) !outbound
call lambert(rv_mars, rv_earth_f,dt_earth_return*day2sec,mu_sun,v3,v4) !return

!compute v-inf vectors:
vinf_0 = v1 - rv_earth_0(4:6) !earth departure
vinfm = v2 - rv_mars(4:6) !mars flyby (incoming)
vinfp = v3 - rv_mars(4:6) !mars flyby (outgoing)
vinf_f = v4 - rv_earth_f(4:6) !earth return

!the turning angle is the angle between vinf- and vinf+
delta = angle_between_vectors(vinfm,vinfp) !turning angle [rad]
vinfmag = norm2(vinfm) !vinf vector mag (incoming) [km/s]
e = one/sin(delta/two) !eccentricity
rp = mu_mars/(vinfmag*vinfmag)*(e-one) !periapsis radius [km]

!the constraints are added to the fitness function as penalties:
if (rp>=rp_mars_min) then
    rp_penalty = zero !not active
else
    rp_penalty = rp_mars_min - rp
end if
vinf_penalty = abs(norm2(vinfm) - norm2(vinfp))

!fitness function (to maximize):
f = - ( norm2(vinf_0) + &
    norm2(vinf_f) + &
    vinf_weight*vinf_penalty + &
    rp_weight*rp_penalty )

end subroutine obj_func

The lambert routine called here is simply a wrapper to solve_lambert_gooding that computes both the "short way" and "long way" transfers, and returns the one with the lowest total \(\Delta v\). The two constraints are added to the objective function as penalties (to be driven to zero). Pikaia is then initialized and called by:

program flyby

use pikaia_module
use fortran_astrodynamics_toolkit

implicit none

integer,parameter :: n = 3

type(pikaia_class) :: p
real(wp),dimension(n) :: x,xl,xu
integer :: status,iseed
real(wp) :: f

!set random number seed:
iseed = 371

!bounds:
xl = [jd_to_mjd(julian_date(2018,1,1,0,0,0)), 100.0_wp,100.0_wp]
xu = [jd_to_mjd(julian_date(2018,12,31,0,0,0)),400.0_wp,400.0_wp]

call p%init(n,xl,xu,obj_func,status,&
            ngen = 1000, &
            np = 200, &
            nd = 9, &
            ivrb = 0, &
            convergence_tol = 1.0e-6_wp, &
            convergence_window = 200, &
            initial_guess_frac = 0.0_wp, &
            iseed = iseed )

call p%solve(x,f,status)

end program flyby

The three optimization variables are bounded: the Earth departure epoch must occur sometime in 2018, and the outbound and return flight times must each be between 100 and 400 days.

Running this program, we can get different solutions (depending on the value set for the random seed, and also the Pikaia settings). The best one I've managed to get with minimal tweaking is this:

• Earth departure epoch: MJD 58120 (Jan 2, 2018)
• Outbound flight time: 229 days
• Return flight time: 271 days

Which is similar to the solution shown in [2].

References

1. R.H. Gooding, "A procedure for the solution of Lambert's orbital boundary-value problem", Celestial Mechanics and Dynamical Astronomy, Vol. 48, No. 2, 1990.
2. Dennis Tito's 2021 Human Mars Flyby Mission Explained. SPACE.com, February 27, 2014.
3. J.E. Prussing, B.A. Conway, "Orbital Mechanics", Oxford University Press, 1993.

Pikaia

I've started a new project on Github for a new modern object-oriented version of the Pikaia genetic algorithm code. This is a refactoring and upgrade of the Pikaia code from the High Altitude Observatory. The original code is public domain and was written by Paul Charbonneau & Barry Knapp (version 1.2 was released in 2003). The new code differs from the old code in the following respects:

• The original fixed-form source (FORTRAN 77) was converted to free-form source.
• The code is now object-oriented Fortran 2003/2008. All user interaction is now through a class, which is the only public entity in the module.
• All real variables are now 8 bytes (a.k.a. "double precision"). The precision can also be changed by the user by modifying the wp parameter.
• The random number generator was replaced with calls to the intrinsic Fortran random_number function.
• There are various new options (e.g., a convergence window with a tolerance can be specified as a stopping condition, and the user can specify a subroutine for reporting the best solution at each generation).
• Mapping the variables to be between 0 and 1 now occurs internally, rather than requiring the user to do it.
• Can now include an initial guess in the initial population.

This is another example of how older Fortran codes can be refactored to bring them up to date with modern computing concepts without doing that much work. In my opinion, the object-oriented Fortran 2003 interface is a lot better than the old FORTRAN 77 version. For one, it's easier to follow. Also, data can now easily be passed into the user fitness function by extending the pikaia class and putting the data there. No common blocks necessary!

References

1. Original Pikaia code and documentation [High Altitude Observatory]