Feature solo_mgar

pykep/trajopt/CMakeLists.txt

@@ -9,5 +9,6 @@ INSTALL(FILES _pl2pl_N_impulses.py DESTINATION ${PYKEP_INSTALL_PATH}/trajopt)
 INSTALL(FILES _indirect.py DESTINATION ${PYKEP_INSTALL_PATH}/trajopt)
 INSTALL(FILES _direct.py DESTINATION ${PYKEP_INSTALL_PATH}/trajopt)
 INSTALL(FILES _lambert.py DESTINATION ${PYKEP_INSTALL_PATH}/trajopt)
-
+INSTALL(FILES _rvt.py DESTINATION ${PYKEP_INSTALL_PATH}/trajopt)
+INSTALL(FILES _resonance.py DESTINATION ${PYKEP_INSTALL_PATH}/trajopt)
 ADD_SUBDIRECTORY(gym)

pykep/trajopt/__init__.py

@@ -17,4 +17,7 @@
 from . import gym 
 from pykep.trajopt._indirect import indirect_or2or, indirect_pt2or, indirect_pt2pt, indirect_pt2pl
 from pykep.trajopt._direct import direct_pl2pl
-from pykep.trajopt._lambert import lambert_problem_multirev, lambert_problem_stochastic
+from pykep.trajopt._lambert import lambert_problem_multirev, lambert_problem_multirev_ga, lambert_problem_stochastic
+from pykep.trajopt._rvt import rvt, rvt_planet, rotate_vector
+from pykep.trajopt._resonance import resonance
+

pykep/trajopt/_lambert.py

@@ -1,6 +1,7 @@
 import numpy as np
 import math
 import random
+from pykep.core import fb_vel
 
 class lambert_problem_multirev:
 
@@ -53,6 +54,58 @@ def get_tof(self):
     def get_Nmax(self):
         return self.lambert_problem.get_Nmax()
 
+class lambert_problem_multirev_ga:    
+    r"""
+    This class converts a lambert_problem instance - a number of solutions to a multi revolution Lambert problem
+    to an instance representing only "the best" solution. Criteria is the delta velocity of the incoming
+    velocity v_in to the outgoing solution velocity get_v1() considering a GA maneuver at a planet.
+    Can be used as a replacement for lambert_problem in _mga.py and _mga_1dsm.py. lambert_problem_multirev delivers 
+    an optimal unique solution comparing all Lambert solutions of a multi revolution lambert problem which can 
+    improve optimization results for trajectories including inner planets.
+    """
+
+    def __init__(self, v_in, lp, planet, v_planet):
+        best_i = 0        
+        n = len(lp.get_v1())
+        if n > 0: 
+            best_dv = math.inf           
+            for i in range(n):
+                vin = [a - b for a, b in zip(v_in, v_planet)]
+                vout = [a - b for a, b in zip(lp.get_v1()[i], v_planet)]
+                dv = fb_vel(vin, vout, planet)
+                if dv < best_dv:
+                    best_dv = dv
+                    best_i = i 
+        self.best_i = best_i
+        self.lambert_problem = lp
+
+    def get_v1(self):
+        return [self.lambert_problem.get_v1()[self.best_i]]
+
+    def get_v2(self):
+        return [self.lambert_problem.get_v2()[self.best_i]]
+
+    def get_r1(self):
+        return self.lambert_problem.get_r1()  
+
+    def get_r2(self):
+        return self.lambert_problem.get_r2()  
+
+    def get_mu(self):
+        return self.lambert_problem.get_mu()
+
+    def get_x(self):
+        return [self.lambert_problem.get_x()[self.best_i]]
+
+    def get_iters(self):
+        return [self.lambert_problem.get_iters()[self.best_i]]
+
+    def get_tof(self):
+        return self.lambert_problem.get_tof()
+
+    def get_Nmax(self):
+        return self.lambert_problem.get_Nmax()
+
 class lambert_problem_stochastic:
     r"""
     This class converts a lambert_problem instance - a number of solutions to a multi revolution Lambert problem -
@@ -81,7 +134,7 @@ def __init__(self, v_in, lambert_problem):
                 dv = 0.001 * np.linalg.norm([a - b for a, b in zip(lambert_problem.get_v1()[i], v_in)])
                 if dv < best_dv:
                     # skip improvement randomly if dv is only slightly better
-                    if i > 0 and random.random() < 1.0 / (1.0 + 2*(best_dv - dv)**2):
+                    if i > 0 and random.random() < 1.0 / (1.0 + 2*(0.001*best_dv - dv)**2):
                         continue
                     best_dv = dv
                     best_i = i

pykep/trajopt/_resonance.py

@@ -0,0 +1,57 @@
+import math
+
+from pykep.core import epoch, fb_prop, SEC2DAY
+from pykep.trajopt import _resonance
+from pykep.trajopt._rvt import rvt
+
+
+class resonance: 
+    """
+    Determines the best "fitting" resonance orbit.
+    """
+
+    def __init__(self, planet, rvt_in, rvt_pl,
+                 resonances=[[1, 1], [5, 4], [4, 3], [3, 2], [5, 3]]
+        ):
+        """
+        Args:
+            - planet (``planet``): resonance planet. 
+            - rvt_in: (``rvt``): incoming orbit.
+            - rvt_pl: (``rvt``): planet orbit.
+            - resonances (``list`` of ``int``): resonance options. 
+        """
+        assert rvt_in._t == rvt_pl._t  # timing must be consistent
+
+        self._planet = planet
+        self._rvt_in = rvt_in
+        self._rvt_pl = rvt_pl       
+        self._time = rvt_in._t
+        self._resonances = resonances
+        self._period = planet.compute_period(epoch(self._time * SEC2DAY))
+        self._mu = planet.mu_self
+        self._timing_error = -1
+        self._rvt_out = None
+        self._resonance = None
+
+    # useful for debugging
+    def __str__(self):
+        return str(_resonance) + " " + str(self._timing_error * SEC2DAY) + " " + str(self._rvt_out)
+
+    # select a resonance option minimizing the timing error  
+    def select_resonance(self, beta, safe_distance):
+        v_out = fb_prop(self._rvt_in._v, self._rvt_pl._v,
+                        self._planet.radius + safe_distance, beta, self._mu)
+        self._rvt_out = rvt(self._rvt_in._r, v_out, self._time, self._rvt_in._mu)
+        period = self._rvt_out.period()
+        self._timing_error = math.inf
+        for resonance in self._resonances:
+            target = self._period * resonance[1] / resonance[0];
+            dt = abs(period - target)
+            if dt < self._timing_error:
+                self._resonance = resonance
+                self._timing_error = dt
+        return self._timing_error, self._resonance
+
+    # time of flight of the resonance transfer
+    def tof(self):
+        return self._resonance[1] * self._period

pykep/trajopt/_rvt.py

@@ -0,0 +1,94 @@
+from math import pi, sqrt, cos, sin
+import math
+
+from pykep.core import AU, RAD2DEG, SEC2DAY
+from pykep.core.core import propagate_lagrangian, ic2par, epoch, \
+    propagate_taylor
+
+import numpy as np
+
+
+class rvt:
+
+    """
+    Keplerian orbit represented by radius, velocity, time and mu.    
+    """
+
+    def __init__(self, r, v, time, mu):
+        """
+        Args:
+            - r (``tuple`` of ``float``): cartesian position in m.
+            - v: (``tuple`` of ``float``): velocity in m/s.
+            - time: (``float``): time in seconds.
+            - mu (`float``): gravity parameter of the central body.
+        """
+
+        self._r = r
+        self._v = v
+        self._t = time  # in seconds
+        self._mu = mu
+
+    # useful for debugging        
+    def __str__(self):
+        a, e, i, _, _, _ = self.kepler()
+        period = 2 * pi * sqrt(a ** 3 / self._mu)
+        apo = a * (1 + e) / AU
+        per = a * (1 - e) / AU
+        return str(self._r) + " " + str(self._v) + " " + str(self._t * SEC2DAY) + " " + \
+                str(apo) + " " + str(per) + " " + \
+                str(e) + " " + str(i * RAD2DEG) + " " + str(period * SEC2DAY)
+
+    def propagate_lagrangian(self, tof):
+        orb = rvt(self._r, self._v, self._t + tof, self._mu)
+        orb._r, orb._v = propagate_lagrangian(orb._r, orb._v, tof, self._mu)
+        return orb
+
+    def propagate_taylor(self, tof, m0, thrust, veff=1, log10tol=-15, log10rtol=-15):
+        orb = rvt(self._r, self._v, self._t + tof, self._mu)
+        orb._r, orb._v, m = propagate_taylor(orb._r, orb._v, m0, thrust, tof, self._mu,
+                                             veff, log10tol, log10rtol)
+        return orb, m
+
+    # keplarian parameters a, e, i, W, w, E       
+    def kepler(self):
+        return ic2par(self._r, self._v, self._mu) 
+
+    # plots orbit from current time up to time + tof    
+    def plot(self, tof, N=60, units=AU, color="b", label=None, axes=None):
+        from pykep.orbit_plots import plot_kepler
+        plot_kepler(r0=self._r, v0=self._v, tof=tof,
+                    mu=self._mu, N=N, units=units, color=color,
+                    label=label, axes=axes)
+
+    def period(self):
+        kep = ic2par(self._r, self._v, self._mu) 
+        a = kep[0]
+        meanMotion = sqrt(self._mu / (a ** 3))
+        return 2.0 * math.pi / meanMotion;  # in seconds
+
+    def rotate(self, k, theta):
+        orb = rvt(self._r, self._v, self._t, self._mu)
+        orb._r = rotate_vector(self._r, k, theta)
+        orb._v = rotate_vector(self._v, k, theta)
+        return orb 
+
+    def tof(self, rvt2):
+        return rvt2._t - self._t
+
+
+def rvt_planet(pl, time):
+    r, v = pl.eph(epoch(time * SEC2DAY))
+    return rvt(r, v, time, pl.mu_central_body)
+
+
+def rotate_vector(v, k, theta):
+    dP = np.dot(k, v)
+    cosTheta = cos(theta)
+    sinTheta = sin(theta)
+    # rotate using Rodrigues rotation formula
+    r_rot = [
+        a * cosTheta + b * sinTheta + c * (1 - cosTheta) * dP
+        for a, b, c in zip(v, np.cross(k, v), k)
+    ]
+    return r_rot
+

pykep/trajopt/gym/CMakeLists.txt

@@ -7,4 +7,5 @@ INSTALL(FILES _emNimp.py DESTINATION ${PYKEP_INSTALL_PATH}/trajopt/gym)
 INSTALL(FILES _eve_mga1dsm.py DESTINATION ${PYKEP_INSTALL_PATH}/trajopt/gym)
 INSTALL(FILES _cassini2.py DESTINATION ${PYKEP_INSTALL_PATH}/trajopt/gym)
 INSTALL(FILES _messenger.py DESTINATION ${PYKEP_INSTALL_PATH}/trajopt/gym)
+INSTALL(FILES _solo_mgar.py DESTINATION ${PYKEP_INSTALL_PATH}/trajopt/gym)
 

pykep/trajopt/gym/__init__.py

@@ -6,4 +6,5 @@
 from ._emNimp import em3imp, em5imp, em7imp
 from ._eve_mga1dsm import eve_mga1dsm, eve_mga1dsm_a, eve_mga1dsm_n
 from ._messenger import messenger
+from ._solo_mgar import solo_mgar, solo_mgar_unconstrained, solo_mgar_multi_objective