python tests on mga

doc/api.rst

@@ -17,8 +17,10 @@ pykep API is designed to maximize its usability. Let us know what you think abou
    flyby
    leg
    planet
-   plot
+   udpla
    propagation
    taylor_adaptive 
-   udpla
    trajopt
+   gym
+   plot
+

doc/gym.rst

@@ -0,0 +1,17 @@
+.. _gym:
+
+Trajectory Optimization Gym
+###########################
+
+.. currentmodule:: pykep.trajopt.gym
+
+MGA problems
+************************
+
+.. autoattribute:: pykep.trajopt.gym.cassini1
+
+This is an MGA problem inspired to the Cassini spacecraft interplanetary transfer to Saturn. 
+The objective of this mission is to reach Saturn and to be captured by its gravity into an orbit having pericenter radius :math:`r_p=108950` km, 
+and eccentricity :math:`e=0.98`. The planetary fly-by sequence considered is E-VVEJ-S (as the one used by Cassini spacecraft). 
+As objective function we use the total :math:`\Delta V` accumulated during the mission, including the launch :math:`\Delta V` and the various :math:`\Delta V` one 
+needs to give at the planets and upon arrival to perform the final orbit injection. 

doc/index.md

@@ -29,4 +29,10 @@ bibliography
 
 tut_basic
 tut_trajopt
-```
+```
+
+```{image} _static/esa_logo.png
+:alt: ESA patch
+:width: 400px
+:align: center
+``` 

doc/udpla.rst

@@ -1,6 +1,6 @@
 .. _udpla:
 
-List of user implemented planets (UDPLAs)
+User defined planets (UDPLAs)
 #########################################
 
 .. currentmodule:: pykep.udpla

pykep/CMakeLists.txt

@@ -2,7 +2,7 @@
 configure_file("${CMAKE_CURRENT_SOURCE_DIR}/_version.py.in" "${CMAKE_CURRENT_SOURCE_DIR}/_version.py" @ONLY)   
 
 # The list of pykep's Python files.
-set(PYKEP_PYTHON_FILES __init__.py _patch_planet.py _version.py test.py)
+set(PYKEP_PYTHON_FILES __init__.py _patch_planet.py _version.py test.py test_trajopt.py)
 
 # Core module.
 Python3_add_library(core MODULE WITH_SOABI 

pykep/docstrings.cpp

@@ -1688,7 +1688,7 @@ final augmented state with :math:`\mathbf x_f = [\mathbf r_f, \mathbf v_f, m_f]`
   \frac{\partial \mathbf {mc}}{\partial \mathbf u}
 
 Returns:
-    :class:`tuple` [:class:`numpy.ndarray`, :class:`numpy.ndarray`, :class:`numpy.ndarray`]: The three gradients. sizes will be (7,7), (7,7) and (7,nseg*3)
+    :class:`tuple` [:class:`numpy.ndarray`, :class:`numpy.ndarray`, :class:`numpy.ndarray`]: The three gradients. sizes will be (7,7), (7,7) and (7,nseg*3 + 1)
 
 Examples:
   >>> import pykep as pk

pykep/test.py

@@ -8,7 +8,7 @@
 ## with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
 
 import unittest as _ut
-
+from .test_trajopt import mga_tests
 
 def float_abs_error(a: float, b: float):
     return abs(a - b)
@@ -313,6 +313,9 @@ def run_test_suite():
     suite.addTest(propagate_test("test_stark"))
     suite.addTest(py_udplas_test("test_tle"))
     suite.addTest(py_udplas_test("test_spice"))
+    suite.addTest(mga_tests("test_construction"))
+    suite.addTest(mga_tests("test_encoding_to_encoding"))
+
     suite.addTest(tl.loadTestsFromTestCase(vsop2013_test))
 
 

pykep/test_trajopt.py

@@ -0,0 +1,56 @@
+## Copyright 2023, 2024 Dario Izzo ([email protected]), Francesco Biscani
+## ([email protected])
+##
+## This file is part of the pykep library.
+##
+## This Source Code Form is subject to the terms of the Mozilla
+## Public License v. 2.0. If a copy of the MPL was not distributed
+## with this file, You can obtain one at http://mozilla.org/MPL/2.0/.
+
+import pykep as pk
+import pygmo as pg
+
+import unittest as _ut
+
+def float_abs_error(a: float, b: float):
+    return abs(a - b)
+
+def float_rel_error(a: float, b: float):
+    return abs(a - b) / abs(a)
+
+class mga_tests(_ut.TestCase):
+    def test_construction(self):
+        import pykep as pk
+        earth = pk.planet(pk.udpla.jpl_lp("earth"))
+        venus = pk.planet(pk.udpla.jpl_lp("venus"))
+        udp = pk.trajopt.mga(
+            seq=[
+                earth,
+                venus,
+                earth,
+                venus,
+                earth
+            ],
+            tof_encoding = "direct",
+            t0=[0, 1000],
+            tof=[[30, 200], [30, 300], [30, 300], [30, 300]],
+            vinf=2.5,
+        )
+        prob = pg.problem(udp)
+        pop = pg.population(prob, 100)
+
+    def test_encoding_to_encoding(self):
+        import pykep as pk
+        udp_direct = pk.trajopt.mga(tof_encoding="direct", tof = [[30, 200], [200, 300]])
+        udp_alpha = pk.trajopt.mga(tof_encoding="alpha", tof = [230, 500])
+        udp_eta = pk.trajopt.mga(tof_encoding="eta", tof = 500)
+        prob = pg.problem(udp_direct)
+        pop = pg.population(prob, 100)
+        x_direct = pop.champion_x
+        gt = udp_direct.fitness(x_direct)[0]
+        x_alpha = udp_alpha.direct2alpha(x_direct)
+        x_eta = udp_eta.direct2eta(x_direct)
+        self.assertTrue(float_rel_error(gt, udp_alpha.fitness(x_alpha)[0]) < 1e-14)
+        self.assertTrue(float_rel_error(gt, udp_eta.fitness(x_eta)[0]) < 1e-14)
+        self.assertTrue(float_rel_error(gt, udp_direct.fitness(udp_direct.alpha2direct(x_alpha))[0]) < 1e-14)
+        self.assertTrue(float_rel_error(gt, udp_direct.fitness(udp_eta.eta2direct(x_eta))[0]) < 1e-14)

pykep/trajopt/CMakeLists.txt

@@ -1,2 +1,5 @@
 set(PYKEP_TRAJOPT_PYTHON_FILES __init__.py _direct_point2point.py _direct_pl2pl.py _mga.py)
 install(FILES ${PYKEP_TRAJOPT_PYTHON_FILES} DESTINATION ${_PYKEP_INSTALL_DIR}/trajopt)
+
+ADD_SUBDIRECTORY(gym)
+

pykep/trajopt/__init__.py

@@ -14,4 +14,7 @@
 from ._direct_pl2pl import direct_pl2pl
 
 # Evolutionary encodings for high energy transfers (chemical propulsion)
-from ._mga import mga
+from ._mga import mga
+
+# The interplanetary trajectory gym
+from . import gym 

pykep/trajopt/_mga.py

@@ -99,6 +99,10 @@ def __init__(
 
         # 3 - Check the tof bounds
         if tof_encoding == "alpha":
+            if type(tof) is not type([]):
+                raise TypeError(
+                    r"When the tof_encoding is 'alpha', the tof must be in the form [lb, ub]."
+                )
             if len(tof) != 2:
                 raise TypeError(
                     r"When the tof_encoding is 'alpha', the tof must be in the form [lb, ub]."
@@ -199,7 +203,8 @@ def _decode_tofs(self, x: List[float]) -> List[float]:
                 T[i] = (dt - sum(T[:i])) * x[i + 1]
             return T
 
-    def alpha2direct(self, x):
+    @staticmethod
+    def alpha2direct(x):
         """alpha2direct(x)
 
         Args:
@@ -213,7 +218,8 @@ def alpha2direct(self, x):
         retval = _np.insert(retval, 0, x[0])
         return retval
 
-    def direct2alpha(self, x):
+    @staticmethod
+    def direct2alpha(x):
         """direct2alpha(x)
 
         Args:
@@ -324,12 +330,13 @@ def _compute_dvs(self, x: List[float]) -> Tuple[
         if self.orbit_insertion:
             # In this case we compute the insertion DV as a single pericenter
             # burn
+            MU_SELF = self.seq[-1].get_mu_self()
             DVper = _np.sqrt(
-                DVarrival * DVarrival + 2 * self.seq[-1].mu_self / self.rp_target
+                DVarrival * DVarrival + 2 * MU_SELF / self.rp_target
             )
             DVper2 = _np.sqrt(
-                2 * self.seq[-1].mu_self / self.rp_target
-                - self.seq[-1].mu_self / self.rp_target * (1.0 - self.e_target)
+                2 * MU_SELF / self.rp_target
+                - MU_SELF / self.rp_target * (1.0 - self.e_target)
             )
             DVarrival = _np.abs(DVper - DVper2)
         return (DVlaunch, DVfb, DVarrival, l, DVlaunch_tot, ep, T)

pykep/trajopt/gym/CMakeLists.txt

@@ -0,0 +1,2 @@
+set(PYKEP_GYM_PYTHON_FILES __init__.py _cassini_mga.py)
+install(FILES ${PYKEP_GYM_PYTHON_FILES} DESTINATION ${_PYKEP_INSTALL_DIR}/trajopt/gym)

pykep/trajopt/gym/__init__.py

@@ -0,0 +1 @@
+from ._cassini_mga import cassini1, cassini1_a, cassini1_n

pykep/trajopt/gym/_cassini_mga.py

@@ -0,0 +1,91 @@
+import pykep as _pk
+from pykep.trajopt import mga as _mga
+
+# CASSINI
+_seq_cassini = [_pk.planet(_pk.udpla.jpl_lp('earth')), 
+                _pk.planet(_pk.udpla.jpl_lp('venus')), 
+                _pk.planet(_pk.udpla.jpl_lp('venus')), 
+                _pk.planet(_pk.udpla.jpl_lp('earth')), 
+                _pk.planet(_pk.udpla.jpl_lp('jupiter')), 
+                _pk.planet(_pk.udpla.jpl_lp('saturn'))]
+
+class _cassini1_udp(_mga):
+    def __init__(self):
+        super().__init__(
+            seq=_seq_cassini,
+            t0=[-1000., 0.],
+            tof=[[30, 400], [100, 470], [30, 400], [400, 2000], [1000, 6000]],
+            vinf=3.,
+            tof_encoding='direct',
+            orbit_insertion=True,
+            e_target=0.98,
+            rp_target=108950000)
+
+    def get_name(self):
+        return "Cassini MGA direct tof encoding (Trajectory Optimisation Gym P1)"
+
+    def get_extra_info(self):
+        retval = "\tTrajectory Optimisation Gym problem (P1): Cassini MGA, single objective, direct encoding\n"
+        retval += "\tPlanetary sequence" + \
+            str([pl.get_name() for pl in _seq_cassini1])
+        return retval
+
+    def __repr__(self):
+        return self.get_name()
+
+
+class _cassini1a_udp(_mga):
+    def __init__(self):
+        super().__init__(
+            seq=_seq_cassini,
+            t0=[-1000., 0.],
+            tof=[4000., 7000.],
+            vinf=3.,
+            tof_encoding='alpha',
+            orbit_insertion=True,
+            e_target=0.98,
+            rp_target=108950000)
+
+    def get_name(self):
+        return "Cassini MGA alpha tof encoding (Trajectory Optimisation Gym P2)"
+
+    def get_extra_info(self):
+        retval = "\tTrajectory Optimisation Gym problem (P2): Cassini MGA, single objective, alpha encoding\n"
+        retval += "\tPlanetary sequence" + \
+            str([pl.get_name() for pl in _seq_cassini1])
+        return retval
+
+    def __repr__(self):
+        return self.get_name()
+
+
+class _cassini1n_udp(_mga):
+    def __init__(self):
+        super().__init__(
+            seq=_seq_cassini,
+            t0=[-1000., 0.],
+            tof=7000.,
+            vinf=3.,
+            tof_encoding='eta',
+            orbit_insertion=True,
+            e_target=0.98,
+            rp_target=108950000)
+
+    def get_name(self):
+        return "Cassini1 MGA eta tof encoding (Trajectory Optimisation Gym P3)"
+
+    def get_extra_info(self):
+        retval = "\tTrajectory Optimisation Gym problem (P3): Cassini MGA, single objective, eta encoding\n"
+        retval += "\tPlanetary sequence" + \
+            str([pl.get_name() for pl in _seq_cassini1])
+        return retval
+
+    def __repr__(self):
+        return self.get_name()
+
+# Problem P1: Cassini MGA, single objective, direct encoding
+cassini1 = _cassini1_udp()
+# Problem P2: Cassini MGA, single objective, alpha encoding
+cassini1_a = _cassini1a_udp()
+# Problem P3: Cassini MGA, single objective, eta encoding
+cassini1_n = _cassini1n_udp()