diff --git a/pykep/trajopt/CMakeLists.txt b/pykep/trajopt/CMakeLists.txt index 6bf336cb..fd522765 100644 --- a/pykep/trajopt/CMakeLists.txt +++ b/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) diff --git a/pykep/trajopt/__init__.py b/pykep/trajopt/__init__.py index 8a4406e1..c2f614e5 100644 --- a/pykep/trajopt/__init__.py +++ b/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 + diff --git a/pykep/trajopt/_lambert.py b/pykep/trajopt/_lambert.py index 34553242..b4d43f8e 100644 --- a/pykep/trajopt/_lambert.py +++ b/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 diff --git a/pykep/trajopt/_resonance.py b/pykep/trajopt/_resonance.py new file mode 100644 index 00000000..aaa17791 --- /dev/null +++ b/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 diff --git a/pykep/trajopt/_rvt.py b/pykep/trajopt/_rvt.py new file mode 100644 index 00000000..09ae8054 --- /dev/null +++ b/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 + diff --git a/pykep/trajopt/gym/CMakeLists.txt b/pykep/trajopt/gym/CMakeLists.txt index 73cc080c..7df55e09 100644 --- a/pykep/trajopt/gym/CMakeLists.txt +++ b/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) diff --git a/pykep/trajopt/gym/__init__.py b/pykep/trajopt/gym/__init__.py index 6092b685..07512d4f 100644 --- a/pykep/trajopt/gym/__init__.py +++ b/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 diff --git a/pykep/trajopt/gym/_solar_orbiter.py b/pykep/trajopt/gym/_solar_orbiter.py new file mode 100644 index 00000000..e419d55f --- /dev/null +++ b/pykep/trajopt/gym/_solar_orbiter.py @@ -0,0 +1,757 @@ +from bisect import bisect_left +from math import cos, pi, sin, sqrt +from typing import Any, List, Tuple + +import numpy as np +from pykep.trajopt._lambert import lambert_problem_multirev + +from pykep import AU, DAY2SEC, DEG2RAD, RAD2DEG, SEC2DAY, epoch, ic2par +from pykep.core import fb_prop, fb_vel, lambert_problem, propagate_lagrangian +from pykep.planet import jpl_lp +from pykep.trajopt import launchers + + +class _solar_orbiter_udp: + + earth = jpl_lp("earth") + venus = jpl_lp("venus") + + def __init__( + self, + t0=[epoch(0), epoch(10000)], + multi_objective=False, + tof_encoding="direct", + max_revs: int = 0, + dummy_DSMs: bool = False, + evolve_rev_count=False, + seq=[earth, venus, venus, earth, venus, venus, venus, venus, venus], + ) -> None: + """ + Args: + - multi_objective (``bool``): when True the problem fitness will return also the time of flight as an added objective + - tof_encoding (``str``): one of 'direct', 'eta' or 'alpha'. Selects the encoding for the time of flights + - tof (``list`` or ``list`` of ``list``): time of flight bounds. As documented in ``pykep.mga_1dsm`` + - max_revs (``int``): maximal number of revolutions for lambert transfer + - dummy_DSMs (``bool``): whether to add deep space maneuvers after flyby + - evolve_rev_count (``bool``): whether to treat the number of revolutions as a evolvable parameter in each leg + - seq (``list``) + + """ + + eta_lb = 0.1 + eta_ub = 0.9 + rp_ub = 30 + safe_distance = 350000 + min_start_mass = 1800 + + # Redefining the planets as to change their safe radius. 350km was given as safe distance. + + # alternative: Launch-Venus-Venus-Earth-Venus + for i in range(len(seq)): + seq[i].safe_radius = (seq[i].radius + safe_distance) / seq[i].radius + + tof = [[10, 900]] * (len(seq) - 1) + + # Sanity checks + # 1 - Planets need to have the same mu_central_body + if [r.mu_central_body for r in seq].count(seq[0].mu_central_body) != len(seq): + raise ValueError( + "All planets in the sequence need to have identical mu_central_body" + ) + # 2 - tof encoding needs to be one of 'alpha', 'eta', 'direct' + if tof_encoding not in ["alpha", "eta", "direct"]: + raise TypeError("tof encoding must be one of 'alpha', 'eta', 'direct'") + # 3 - tof is expected to have different content depending on the tof_encoding + if tof_encoding == "direct": + if np.shape(np.array(tof)) != (len(seq) - 1, 2): + raise TypeError( + "tof_encoding is " + + tof_encoding + + " and tof must be a list of two dimensional lists and with length equal to the number of legs" + ) + if tof_encoding == "alpha": + if np.shape(np.array(tof)) != (2,): + raise TypeError( + "tof_encoding is " + + tof_encoding + + " and tof must be a list of two floats" + ) + if tof_encoding == "eta": + if np.shape(np.array(tof)) != (): + raise TypeError( + "tof_encoding is " + tof_encoding + " and tof must be a float" + ) + # 4 - Check launch window t0. If defined in terms of floats transform into epochs + if len(t0) != 2: + raise TypeError( + "t0 is " + t0 + " while should be a list of two floats or epochs" + ) + if type(t0[0]) is not epoch: + t0[0] = epoch(t0[0]) + if type(t0[1]) is not epoch: + t0[1] = epoch(t0[1]) + + # 5 - Resonant Flybys mess up the lambert arcs. Watch out for planet sequences that could have them: + possibly_resonant: List[bool] = [ + dummy_DSMs and seq[i] == seq[i + 1] for i in range(len(seq) - 1) + ] + if possibly_resonant[0]: + raise NotImplementedError("Resonant flybys not yet supported at launch.") + + self._seq = seq + self._t0 = t0 + self._tof = tof + self._tof_encoding = tof_encoding + self._multi_objective = multi_objective + self._max_revs = max_revs + self._eta_lb = eta_lb + self._eta_ub = eta_ub + self._rp_ub = rp_ub + self._safe_distance = safe_distance + self._min_start_mass = min_start_mass + self._dummy_DSM = possibly_resonant + self._evolve_rev_count = evolve_rev_count + + self._n_legs = len(seq) - 1 + self._common_mu = seq[0].mu_central_body + + # initialize data to compute heliolatitude + t_plane_crossing = epoch(7645) + rotation_axis = seq[0].eph(t_plane_crossing)[0] + self._rotation_axis = rotation_axis / np.linalg.norm(rotation_axis) + self._theta = 7.25 * DEG2RAD + + self._eph_cache = (None, None, None) + + def _decode_tofs(self, x: List[float]) -> List[float]: + tail = 3 * sum(self._dummy_DSM) + self._n_legs * self._evolve_rev_count + 2 + if self._tof_encoding == "alpha": + # decision vector is [t0, T, a1, a2, ....] + T = np.log(x[2:-tail]) + return T / sum(T) * x[1] + elif self._tof_encoding == "direct": + # decision vector is [t0, T1, T2, T3, ... ] + return x[1:-tail] + elif self._tof_encoding == "eta": + # decision vector is [t0, n1, n2, n3, ... ] + dt = self._tof + T = [0] * self._n_legs + T[0] = dt * x[1] + for i in range(1, len(T)): + T[i] = (dt - sum(T[:i])) * x[i + 1] + return T + else: + raise TypeError("tof encoding must be one of 'alpha', 'eta', 'direct'") + + def _compute_dvs( + self, x: List[float] + ) -> Tuple[ + List[float], + List[Any], + List[float], + List[Tuple[Tuple[float, float, float], Tuple[float, float, float]]], + List[float], + ]: + + # 1 - we 'decode' the times of flights and compute epochs (mjd2000) + T = self._decode_tofs(x) # [T1, T2 ...] + ep = np.insert(T, 0, x[0]) # [t0, T1, T2 ...] + ep = np.cumsum(ep) # [t0, t1, t2, ...] + + if not len(ep) == len(self._seq): + raise ValueError( + "Got " + + str(len(ep)) + + " epochs, but sequence of length " + + str(len(self._seq)) + ) + # 2 - we compute the ephemerides + r = [(0, 0, 0)] * len(self._seq) + v = [(0, 0, 0)] * len(self._seq) + for i in range(len(self._seq)): + r[i], v[i] = self._seq[i].eph(ep[i]) + + if self._tof_encoding == "alpha": + tof_offset = 2 + self._n_legs + elif self._tof_encoding in ["direct", "eta"]: + tof_offset = 1 + self._n_legs + else: + assert False + + lambert_indices = [0] * self._n_legs + if self._evolve_rev_count: + lambert_indices = [ + int(elem) for elem in x[tof_offset + sum(self._dummy_DSM) * 3 : -2] + ] + assert len(lambert_indices) == self._n_legs + for index in lambert_indices: + if (not index >= 0) or index >= 2 * self._max_revs + 1: + raise ValueError("Invalid lambert index " + str(index)) + + rf_index = np.cumsum(self._dummy_DSM) + # 3 - we solve the lambert problems + l = list() + DVdsm = list() + ballistic_legs: List[ + Tuple[Tuple[float, float, float], Tuple[float, float, float]] + ] = list() + ballistic_ep: List[float] = list() + + v_probe = v[0] + for i in range(self._n_legs): + # start leg at planet i, before flyby + ri = r[i] + Ti = T[i] + + if np.any(np.isnan(v_probe)): + return ([np.nan], [], ep, [], []) + + if self._dummy_DSM[i]: + # insert dummy DSM + + flyby_param_index = tof_offset + (rf_index[i] - 1) * 3 + + assert flyby_param_index + 5 <= len(x) + + tof_ratio, beta, r_p = x[ + flyby_param_index : flyby_param_index + 3 + ] # TODO: adapt this to different encodings + if tof_ratio < 0 or tof_ratio > 1: + raise ValueError( + "Time of flight ratio of " + str(tof_ratio) + " is invalid." + ) + if beta < -2 * pi or beta > 2 * pi: + raise ValueError("Invalid flyby angle beta: " + str(beta)) + if r_p < self._seq[i].safe_radius / self._seq[i].radius: + raise ValueError("Invalid flyby periapsis: " + str(r_p)) + + # perform unpowered flyby + v_probe = fb_prop( + v_probe, v[i], r_p * self._seq[i].radius, beta, self._seq[i].mu_self + ) # TODO: is seq[i] the correct planet? + ballistic_legs.append((ri, v_probe)) + ballistic_ep.append(ep[i]) + + if np.any(np.isnan(v_probe)): + return ([np.nan], [], ep, [], []) + + if tof_ratio > 0: + # propagate after flyby + try: + ri, v_probe = propagate_lagrangian( + ri, v_probe, T[i] * DAY2SEC * tof_ratio, self._common_mu + ) + except RuntimeError as e: + print(e.args) + return ([np.nan], [], ep, [], []) + + # adapt the remaining time for the lambert leg + Ti = Ti * (1 - tof_ratio) + if Ti == 0: + Ti = SEC2DAY + + if np.any(np.isnan(v_probe)): + return ([np.nan], [], ep, [], []) + + # call lambert solver on remaining leg - either after flyby or DSM + lp = lambert_problem( + ri, r[i + 1], Ti * DAY2SEC, self._common_mu, False, self._max_revs + ) + + # the lambert solver might offer fewer solutions than asked for + lambert_index = min(lp.get_Nmax() * 2, lambert_indices[i]) + if lambert_index < 0 or not lambert_index < len(lp.get_v1()): + raise ValueError( + "Lambert leg has " + + lp.get_Nmax() + + " revolutions but only " + + len(lp.get_v1()) + + " solutions." + ) + + if not self._evolve_rev_count: + lp = lambert_problem_multirev(v_probe, lp) + assert lambert_index == 0 + l.append(lp) + + # add delta v of DSM + if self._dummy_DSM[i]: + DVdsm.append( + np.linalg.norm( + [a - b for a, b in zip(v_probe, lp.get_v1()[lambert_index])] + ) + ) + + v_probe = lp.get_v2()[lambert_index] + ep_start_lambert = ep[i + 1] - Ti + # ri is now either the position of planet i or the position of the DSM + ballistic_legs.append((ri, lp.get_v1()[lambert_index])) + ballistic_ep.append(ep_start_lambert) + + # add ballistic leg after final flyby + final_lambert_index = min(lambert_indices[-1], l[-1].get_Nmax() * 2) + eph = self._seq[-1].eph(ep[-1]) + v_out = fb_prop( + l[-1].get_v2()[final_lambert_index], + eph[1], + x[-1] * self._seq[-1].radius, + x[-2], + self._seq[-1].mu_self, + ) + ballistic_legs.append((eph[0], v_out)) + ballistic_ep.append(ep[-1]) + + # 4 - we compute the various dVs needed at fly-bys to match incoming + # and outcoming, also delta v of deep space maneuvers + assert len(DVdsm) == sum(self._dummy_DSM) + DV = list() + j = 0 # index of DSM + for i in range(len(l) - 1): + if self._dummy_DSM[i + 1]: + # use delta v of deep space maneuver + DV.append(DVdsm[j]) + j += 1 + else: + # use delta v of flyby + # the lambert solver might offer fewer solutions than asked for + lambert_index_incoming = min(l[i].get_Nmax() * 2, lambert_indices[i]) + lambert_index_outgoing = min( + l[i + 1].get_Nmax() * 2, lambert_indices[i + 1] + ) + + vin = [ + a - b + for a, b in zip(l[i].get_v2()[lambert_index_incoming], v[i + 1]) + ] + + vout = [ + a - b + for a, b in zip(l[i + 1].get_v1()[lambert_index_outgoing], v[i + 1]) + ] + DV.append(fb_vel(vin, vout, self._seq[i + 1])) + assert j == len(DVdsm) + assert len(ballistic_legs) == sum(self._dummy_DSM) + len(l) + 1 + assert len(ballistic_ep) == len(ballistic_legs) + assert len(DV) == len(l) - 1 + return (DV, l, ep, ballistic_legs, ballistic_ep) + + def _rotate_vector(self, v, k, theta): + dP = np.dot(k, v) + cosTheta = cos(theta) + sinTheta = sin(theta) + # rotate vector into coordinate system defined by the sun's equatorial plane + # 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 + + # Objective function + def fitness(self, x): + if len(x) != len(self.get_bounds()[0]): + raise ValueError( + "Expected " + + str(len(self.get_bounds()[0])) + + " parameters but got " + + str(len(x)) + ) + + lower_bound, upper_bound = self.get_bounds() + + for i in range(len(x)): + if x[i] < lower_bound[i] or x[i] > upper_bound[i]: + return [np.inf] + [np.inf] * self._multi_objective + [np.nan] * 4 + + DV, lamberts, ep, b_legs, b_ep = self._compute_dvs(x) + T = self._decode_tofs(x) + + if np.any(np.isnan(DV)): + return [np.inf] + [T] * self._multi_objective + [np.nan] * 4 + + # compute launch velocity and declination + Vinfx, Vinfy, Vinfz = [ + a - b for a, b in zip(b_legs[0][1], self._seq[0].eph(ep[0])[1]) + ] + Vinf_launch = np.linalg.norm([Vinfx, Vinfy, Vinfz]) + + # We now have the initial mass of the spacecraft + m_initial = launchers.atlas551(Vinf_launch / 1000.0) + + # compute final flyby and resulting trajectory + eph = self._seq[-1].eph(ep[-1]) + v_out = b_legs[-1][1] + + # rotate to get inclination respective to solar equator. This could be moved to the init phase, as it won't change + r_rot = self._rotate_vector(eph[0], self._rotation_axis, self._theta) + v_rot = self._rotate_vector(v_out, self._rotation_axis, self._theta) + + a, e, i, W, w, E = ic2par(r_rot, v_rot, self._common_mu) + final_perihelion = a * (1 - e) + # orbit should be as polar as possible, but we do not care about prograde/retrograde + corrected_inclination = abs(abs(i) % pi - pi / 2) + + # check perihelion and aphelion bounds during the flight + min_sun_distance = final_perihelion + max_sun_distance = AU + + for l_i in range(len(b_legs) - 1): + # project lambert leg, compute perihelion and aphelion + ri, vi = b_legs[l_i] + + # check transfer points for min and max sun distance + min_sun_distance = min(min_sun_distance, np.linalg.norm(ri)) + max_sun_distance = max(max_sun_distance, np.linalg.norm(ri)) + + transfer_a, transfer_e, _, _, _, E = ic2par(ri, vi, self._common_mu) + transfer_period = 2 * pi * sqrt(transfer_a ** 3 / self._common_mu) + + # check whether extremum happens during this leg + M = E - transfer_e * sin(E) + mean_angle_to_apoapsis = pi - M + if mean_angle_to_apoapsis < 0: + mean_angle_to_apoapsis += 2 * pi + mean_angle_to_periapsis = 2 * pi - M + + # update min and max sun distance + if b_ep[l_i] - b_ep[l_i + 1] > mean_angle_to_apoapsis * transfer_period: + max_sun_distance = max(max_sun_distance, transfer_a * (1 + transfer_e)) + + if b_ep[l_i] - b_ep[l_i + 1] > mean_angle_to_periapsis * transfer_period: + min_sun_distance = min(min_sun_distance, transfer_a * (1 - transfer_e)) + + return ( + [ + corrected_inclination + min_sun_distance / AU + ] # TODO: consider changing this fitness to the one from ESOC + + [sum(T)] * self._multi_objective + + [np.sum(DV) - 0.1] + + [self._min_start_mass - m_initial] + + [0.28 - min_sun_distance / AU] + + [max_sun_distance / AU - 1.2] + ) + + def get_nobj(self): + return self._multi_objective + 1 + + def get_bounds(self): + t0 = self._t0 + tof = self._tof + n_legs = self._n_legs + + if self._tof_encoding == "alpha": + # decision vector is [t0, T, a1, a2, ....] + lb = [t0[0].mjd2000, tof[0]] + [1e-3] * (n_legs) + ub = [t0[1].mjd2000, tof[1]] + [1.0 - 1e-3] * (n_legs) + elif self._tof_encoding == "direct": + # decision vector is [t0, T1, T2, T3, ... ] + lb = [t0[0].mjd2000] + [it[0] for it in self._tof] + ub = [t0[1].mjd2000] + [it[1] for it in self._tof] + elif self._tof_encoding == "eta": + # decision vector is [t0, n1, n2, ....] + lb = [t0[0].mjd2000] + [1e-3] * (n_legs) + ub = [t0[1].mjd2000] + [1.0 - 1e-3] * (n_legs) + + # something of a hack: parameters for dummy DSMs of possibly resonant flybys + for i_fl in range(self._n_legs): + if self._dummy_DSM[i_fl]: + pl = self._seq[i_fl] + lb = lb + [0, -2 * pi, (pl.radius + self._safe_distance) / pl.radius] + ub = ub + [1, 2 * pi, 30] + + # something of a hack: parameters for evolving the lambert index + if self._evolve_rev_count: + lb = lb + [0] * self._n_legs + ub = ub + [2 * self._max_revs] * self._n_legs + + # add final flyby + pl = self._seq[-1] + lb = lb + [-2 * pi, (pl.radius + self._safe_distance) / pl.radius] + ub = ub + [2 * pi, 30] + + if self._tof_encoding == "alpha": + assert ( + len(lb) + == 2 + + n_legs + + 3 * sum(self._dummy_DSM) + + self._n_legs * self._evolve_rev_count + + 2 + ) + elif self._tof_encoding in ["direct", "eta"]: + assert ( + len(lb) + == 1 + + n_legs + + 3 * sum(self._dummy_DSM) + + self._n_legs * self._evolve_rev_count + + 2 + ) + else: + assert False + + assert len(ub) == len(lb) + + return (lb, ub) + + def get_nic(self): + return 4 + + def eph( + self, x: List[float], t: float + ) -> Tuple[Tuple[float, float, float], Tuple[float, float, float]]: + if len(x) != len(self.get_bounds()[0]): + raise ValueError( + "Expected chromosome of length " + + str(len(self.get_bounds()[0])) + + " but got length " + + str(len(x)) + ) + + comparison = x == self._eph_cache[0] + if type(comparison) is not bool: + comparison = all(comparison) + if comparison: + _, b_legs, b_ep = self._eph_cache + else: + _, _, _, b_legs, b_ep = self._compute_dvs(x) + self._eph_cache = (x, b_legs, b_ep) + + if t < b_ep[0]: + raise ValueError( + "Given epoch " + str(t) + " is before launch date " + str(b_ep[0]) + ) + + if t == b_ep[0]: + # exactly at launch + return self._seq[0].eph(t) + + i = bisect_left(b_ep, t) # ballistic leg i goes from planet i to planet i+1 + + assert i >= 1 and i <= len(b_ep) + if i < len(b_ep): + assert t <= b_ep[i] + + # get start of ballistic leg + r_b, v_b = b_legs[i - 1] + + elapsed_seconds = (t - b_ep[i - 1]) * DAY2SEC + assert elapsed_seconds >= 0 + + # propagate the lagrangian + r, v = propagate_lagrangian(r_b, v_b, elapsed_seconds, self._common_mu) + + return r, v + + def pretty(self, x): + """pretty(x) + + Args: + - x (``list``, ``tuple``, ``numpy.ndarray``): Decision chromosome, e.g. (``pygmo.population.champion_x``). + + Prints human readable information on the trajectory represented by the decision vector x + """ + T = self._decode_tofs(x) + DV, lambert_legs, ep, b_legs, b_ep = self._compute_dvs(x) + b_i = 0 + Vinfx, Vinfy, Vinfz = [ + a - b for a, b in zip(b_legs[b_i][1], self._seq[0].eph(ep[0])[1]) + ] + + if self._tof_encoding == "alpha": + tof_offset = 2 + self._n_legs + elif self._tof_encoding in ["direct", "eta"]: + tof_offset = 1 + self._n_legs + else: + assert False + + lambert_indices = [0] * self._n_legs + if self._evolve_rev_count: + lambert_indices = [ + int(elem) for elem in x[tof_offset + sum(self._dummy_DSM) * 3 : -2] + ] + assert len(lambert_indices) == len(lambert_legs) + lambert_indices = [ + min(index, 2 * leg.get_Nmax()) + for index, leg in zip(lambert_indices, lambert_legs) + ] + else: + # we assume that the lambert_problem_multirev class is used + lambert_indices = [lam.best_i for lam in lambert_legs] + assert len(lambert_indices) == self._n_legs + + print("Multiple Gravity Assist (MGA) problem: ") + print("Planet sequence: ", [pl.name for pl in self._seq]) + + print("Departure: ", self._seq[0].name) + print("\tEpoch: ", ep[0], " [mjd2000]") + print("\tSpacecraft velocity: ", b_legs[0][1], "[m/s]") + print("\tLaunch velocity: ", [Vinfx, Vinfy, Vinfz], "[m/s]") + _, _, transfer_i, _, _, _ = ic2par(*(b_legs[0]), self._common_mu) + print("\tOutgoing Inclination:", transfer_i * RAD2DEG, "[deg]") + print("\tNumber of Revolutions:", int((lambert_indices[0] + 1) / 2)) + print("\tLambert Index:", int(lambert_indices[0])) + b_i += ( + 1 + self._dummy_DSM[0] + ) # increasing leg index by one, as the launch contains no DSM + + assert len(DV) == len(self._seq) - 2 + for i in range(1, len(self._seq) - 1): + pl = self._seq[i] + e = ep[i] + dv = DV[i - 1] + leg = b_legs[b_i] + print("Fly-by: ", pl.name) + print("\tEpoch: ", e, " [mjd2000]") + if self._dummy_DSM[i]: + print("\tDSM at ", b_ep[b_i + 1]) + print("\tDV: ", dv, "[m/s]") + eph = pl.eph(e) + assert np.linalg.norm([a - b for a, b in zip(leg[0], eph[0])]) < 0.01 + _, _, transfer_i, _, _, _ = ic2par(eph[0], leg[1], self._common_mu) + print("\tOutgoing Inclination:", transfer_i * RAD2DEG, "[deg]") + print("\tNumber of Revolutions:", int((lambert_indices[i] + 1) / 2)) + print("\tLambert Index:", int(lambert_indices[i])) + b_i += 1 + self._dummy_DSM[i] + + print("Final Fly-by: ", self._seq[-1].name) + print("\tEpoch: ", ep[-1], " [mjd2000]") + print("\tSpacecraft velocity: ", lambert_legs[-1].get_v2()[0], "[m/s]") + print("\tBeta: ", x[-2]) + print("\tr_p: ", x[-1]) + + print("Resulting Solar orbit:") + r_rot = self._rotate_vector(b_legs[-1][0], self._rotation_axis, self._theta) + v_rot = self._rotate_vector(b_legs[-1][1], self._rotation_axis, self._theta) + + a, e, i, W, w, E = ic2par(r_rot, v_rot, self._common_mu) + print("Perihelion: ", (a * (1 - e)) / AU, " AU") + print("Aphelion: ", (a * (1 + e)) / AU, " AU") + print("Inclination: ", i * RAD2DEG, " degrees") + + print("Time of flights: ", T, "[days]") + + def plot(self, x, axes=None, units=AU, N=60): + """plot(self, x, axes=None, units=pk.AU, N=60) + + Plots the spacecraft trajectory. + + Args: + - x (``tuple``, ``list``, ``numpy.ndarray``): Decision chromosome. + - axes (``matplotlib.axes._subplots.Axes3DSubplot``): 3D axes to use for the plot + - units (``float``, ``int``): Length unit by which to normalise data. + - N (``float``): Number of points to plot per leg + """ + import matplotlib as mpl + import matplotlib.pyplot as plt + from mpl_toolkits.mplot3d import Axes3D + + from pykep.orbit_plots import plot_kepler, plot_lambert, plot_planet + + # Creating the axes if necessary + if axes is None: + mpl.rcParams["legend.fontsize"] = 10 + fig = plt.figure() + axes = fig.gca(projection="3d") + + T = self._decode_tofs(x) + ep = np.insert(T, 0, x[0]) # [t0, T1, T2 ...] + ep = np.cumsum(ep) # [t0, t1, t2, ...] + _, l, _, _, _ = self._compute_dvs(x) + for pl, e in zip(self._seq, ep): + plot_planet( + pl, epoch(e), units=units, legend=True, color=(0.7, 0.7, 1), axes=axes + ) + for lamb in l: + plot_lambert( + lamb, + N=N, + sol=0, + units=units, + color="k", + legend=False, + axes=axes, + alpha=0.8, + ) + + # compute final flyby + r_P, v_P = self._seq[-1].eph(ep[-1]) + v_out = fb_prop( + l[-1].get_v2()[0], + v_P, + x[-1] * self._seq[-1].radius, + x[-2], + self._seq[-1].mu_self, + ) + + a, e, i, W, w, E = ic2par(r_P, v_out, self._common_mu) + + # final trajectory + plot_kepler( + r_P, + v_out, + 365 * DAY2SEC, + self._common_mu, + N=100, + color="r", + units=units, + axes=axes, + label="Final Orbit", + ) + + return axes + + def plot_distance_and_flybys(self, x, axes=None, units=AU, N=200, extension=300): + import matplotlib.pyplot as plt + + T = self._decode_tofs(x) + ep = np.insert(T, 0, x[0]) # [t0, T1, T2 ...] + ep = np.cumsum(ep) # [t0, t1, t2, ...] + + timeframe = np.linspace(0, sum(T) + extension, N) + + earth = self._seq[0] + venus = self._seq[-1] + + distances = [] + edistances = [] + vdistances = [] + + for day in timeframe: + t = x[0] + day + pos, vel = self.eph(x, t) + epos, evel = earth.eph(t) + vpos, vvel = venus.eph(t) + distances.append(np.linalg.norm(pos) / AU) + edistances.append(np.linalg.norm(epos) / AU) + vdistances.append(np.linalg.norm(vpos) / AU) + + fl_times = list() + fl_distances = list() + for pl, t in zip(self._seq, ep): + fl_times.append(t - x[0]) + pos, _ = pl.eph(t) + fl_distances.append(np.linalg.norm(pos) / AU) + + if axes is None: + fig, axes = plt.subplots() + axes.plot(list(timeframe), distances, label="Solar Orbiter") + axes.plot(list(timeframe), edistances, label="Earth") + axes.plot(list(timeframe), vdistances, label="Venus") + plt.scatter(fl_times, fl_distances, marker="o", color="r") + axes.set_xlabel("Days") + axes.set_ylabel("AU") + axes.set_title("Distance to Sun") + axes.legend() + + return axes + + +solar_orbiter = _solar_orbiter_udp(max_revs=5, dummy_DSMs=False, evolve_rev_count=False) +solar_orbiter_dsm = _solar_orbiter_udp( + max_revs=5, dummy_DSMs=True, evolve_rev_count=False +) +solar_orbiter_evolve_rev = _solar_orbiter_udp( + max_revs=5, dummy_DSMs=False, evolve_rev_count=True +) diff --git a/pykep/trajopt/gym/_solo_mgar.py b/pykep/trajopt/gym/_solo_mgar.py new file mode 100644 index 00000000..f42f10bc --- /dev/null +++ b/pykep/trajopt/gym/_solo_mgar.py @@ -0,0 +1,476 @@ +from math import pi, sqrt +import math +from typing import Any, List, Tuple + +from pykep import AU, DAY2SEC, SEC2DAY, DEG2RAD, RAD2DEG, epoch, ic2par +from pykep.core import fb_vel, lambert_problem +from pykep.trajopt._lambert import lambert_problem_multirev_ga +from pykep.planet import jpl_lp + +import numpy as np +from pykep.trajopt._resonance import resonance +from pykep.trajopt._rvt import rvt_planet, rotate_vector, rvt + +class _solo_mgar_udp: + """ + See https://www.esa.int/Science_Exploration/Space_Science/Solar_Orbiter . + This UDP represents some of the objectives of the original solar orbiter mission: + - ballistic transfers + - multiple GA manoevers at Venus and Earth + - limited mission time + - inclination > 33 deg at the end of the mission + - low perhelion near 0.28 AU during the EMP phase. + thereby fulfilling its constraints. + To keep it as simple as possible, this UDP currently works with the fixed original planet sequence EVVEVVVVVV, + more flexible variants may be created later. Each EV / VE transfer is performed as MGA, + each VV transfer is performed as resonance. This means, for a free beta parameter the corresponding orbiter + period is computed. From the given resonance alternatives the one is chosen, which minimizes the timing error + - how much earlier / later than venus the orbiter arrives after the resonance transfer. + Instead of "caching" of position / velocity data 'rvt' objects are used: Orbit representations + storing position, velocity, time and mu. + The UDP provides a hand crafted 'unconstrain' activated if use_constraints == False which was used to compute the + reference results at https://github.com/dietmarwo/fast-cma-es/blob/master/examples/data/solo_results.txt + An example for test/optimize is here: https://gist.github.com/dietmarwo/5b6bbadbd696d36806177be2df3c6bb7 + """ + + def __init__( + self, + t0_limits=[7000, 8000], + tof_limits=[[50, 420], [50, 400], [50, 400]], + max_revs: int=2, + resonances= + [[[1, 1], [5, 4], [4, 3]], + [[1, 1], [5, 4], [4, 3]], + [[1, 1], [5, 4], [4, 3]], + [[4, 3], [3, 2], [5, 3]], + [[4, 3], [3, 2], [5, 3]], + [[4, 3], [3, 2], [5, 3]]], + safe_distance=350000, + max_mission_time=11.0 * 365.25, + max_dv0=5600, + multi_objective=False, + use_constraints=True + ): + """ + Args: + - t0_limits (``list`` of ``float``): start time bounds. + - tof_limits (``list`` of ``list`` of ``float``): time of flight bounds. + - max_revs (``int``): maximal number of revolutions for Lambert transfer. + - resonances (``list`` of ``list`` of ``int``): resonance options. + - safe_distance: (``float``): safe distance from planet at GA maneuver in m. + - max_mission_time: (``float``): max mission time in days. + - max_dv0: (``float``): maximal delta velocity at start. + - multi_objective (``bool``): when True the problem fitness will return separate objectives + - use_constraints (``bool``): when False the problem fitness will use a customized "unconstrain" function + """ + + self._safe_distance = safe_distance + self._max_mission_time = max_mission_time + self._max_dv0 = max_dv0 + self._min_beta = -math.pi + self._max_beta = math.pi + + self._earth = jpl_lp("earth") + self._venus = jpl_lp("venus") + self._seq = [self._earth, self._venus, self._venus, self._earth, self._venus, + self._venus, self._venus, self._venus, self._venus, self._venus] + + assert len(self._seq) - 4 == len(resonances) # one resonance option selection for each VV sequence + + self._resonances = resonances + self._t0 = t0_limits + self._tof = tof_limits + self._max_revs = max_revs + self._multi_objective = multi_objective + self._use_constraints = use_constraints + + self._n_legs = len(self._seq) - 1 + + # initialize data to compute heliolatitude + t_plane_crossing = epoch(7645) + rotation_axis = self._seq[0].eph(t_plane_crossing)[0] + self._rotation_axis = rotation_axis / np.linalg.norm(rotation_axis) + self._theta = -7.25 * DEG2RAD # fixed direction of the rotation + + def _compute_dvs(self, x: List[float], + lps=None, + resos=None + ) -> Tuple[ + List[Any], + List[Any], + List[Any], + List[float], + List[float], + ]: + + t0 = x[0] * DAY2SEC # only direct encoding + tof01 = x[1] * DAY2SEC + tof23 = x[2] * DAY2SEC + tof34 = x[3] * DAY2SEC + betas = x[4:] + + rvt_outs = [] + rvt_ins = [None] # no rvt_in at first planet + rvt_pls = [] + dvs = [] + used_resos = [] + reso_dts = [] + + rvt_pls.append(rvt_planet(self._seq[0], t0)) + + _dv_mga(self._seq[0], self._seq[1], tof01, self._max_revs, rvt_outs, rvt_ins, rvt_pls, dvs, lps) + + compute_resonance(self._seq[1], self._resonances[0], betas[0], self._safe_distance, + used_resos, reso_dts, rvt_outs, rvt_ins, rvt_pls, dvs, resos) + + _dv_mga(self._seq[2], self._seq[3], tof23, self._max_revs, rvt_outs, rvt_ins, rvt_pls, dvs, lps) + + _dv_mga(self._seq[3], self._seq[4], tof34, self._max_revs, rvt_outs, rvt_ins, rvt_pls, dvs, lps) + + for i in range(1, 6): + compute_resonance(self._seq[i + 3], self._resonances[i], betas[i], self._safe_distance, + used_resos, reso_dts, rvt_outs, rvt_ins, rvt_pls, dvs, resos) + + n = len(rvt_ins) + assert len(rvt_outs) == n - 1 + assert len(rvt_pls) == n + assert len(rvt_ins) == n + assert len(dvs) == n - 1 + return rvt_outs, rvt_ins, rvt_pls, reso_dts, dvs + + # Objective function + def fitness(self, x): + if len(x) != len(self.get_bounds()[0]): + raise ValueError( + "Expected " + +str(len(self.get_bounds()[0])) + +" parameters but got " + +str(len(x)) + ) + + rvt_outs, rvt_ins, _, reso_dts, dvs = self._compute_dvs(x) + # compute final flyby and resulting trajectory + rvt_out = rvt_outs[-1].rotate(self._rotation_axis, self._theta) # rotate + _, _, incl, _, _, _ = rvt_out.kepler() + # orbit should be as polar as possible, but we do not care about prograde/retrograde + corrected_inclination = abs(abs(incl) % pi - pi / 2) * RAD2DEG + + # check perihelion and aphelion bounds during the flight + emp_perhelion = 2 * AU + min_sun_distance = 2 * AU + max_sun_distance = 0 + + for i in range(len(rvt_outs)): + orb = rvt_outs[i] + tof = orb.tof(rvt_ins[i + 1]) + transfer_a, transfer_e, _, _, _, _ = orb.kepler() + transfer_period = 2 * pi * sqrt(transfer_a ** 3 / orb._mu) + perhelion = transfer_a * (1 - transfer_e) + # update min and max sun distance + if i >= len(rvt_outs) - 3: + emp_perhelion = min(emp_perhelion, perhelion) + min_sun_distance = min(min_sun_distance, perhelion) + if tof > transfer_period: + max_sun_distance = max(max_sun_distance, transfer_a * (1 + transfer_e)) + + time_all = SEC2DAY * (rvt_ins[-1]._t - rvt_outs[0]._t) + time_limit = self._max_mission_time + # wrong reso timing in seconds + reso_penalty = np.sum(reso_dts) + # permit start dv _max_dv0 + dvs[0] = max(0, dvs[0] - self._max_dv0) + + if not self._use_constraints: + # use customized 'unconstrain' to obtain a single objective + dv_val = np.sum(dvs) + time_val = time_all + if time_val > time_limit: + time_val += 10 * (time_val - time_limit) + distance_val = min_sun_distance / AU + # penalty for perhelion < 0.28 + if distance_val < 0.28: + distance_val += 10 * (0.28 - distance_val); + # wrong minimal / maximal distance + distance_penalty = max(0, 0.28 - min_sun_distance / AU) + distance_penalty += max(0, max_sun_distance / AU - 1.2) + value = (100 * dv_val + + 100 * corrected_inclination + + 5000 * (max(0, distance_val - 0.28)) + + 5000 * (max(0, emp_perhelion / AU - 0.28)) + + 0.5 * time_val + + reso_penalty + + 50000 * distance_penalty + ) + else: + corrected_inclination *= DEG2RAD + value = corrected_inclination if self._multi_objective else \ + corrected_inclination + emp_perhelion / AU + return ( + # objectives + [value] + +[emp_perhelion / AU] * self._multi_objective + +[time_all] * self._multi_objective + # constraints + +[np.sum(dvs) - 0.1] * self._use_constraints + +[reso_penalty - 0.1] * self._use_constraints + +[0.28 - min_sun_distance / AU] * self._use_constraints + +[max_sun_distance / AU - 1.2] * self._use_constraints + ) + + def get_nobj(self): + return self._multi_objective * 2 + 1 + + def get_nic(self): + return self._use_constraints * 4 + + def get_bounds(self): + t0 = self._t0 + tof = self._tof + n_legs = self._n_legs + + lb = [t0[0]] + ub = [t0[1]] + lb += [t[0] for t in tof] + ub += [t[1] for t in tof] + + nbetas = n_legs - len(tof) + lb += [self._min_beta] * nbetas + ub += [self._max_beta] * nbetas + return (lb, ub) + + def pretty(self, x): + lambert_legs = [] + resos = [] + rvt_outs, rvt_ins, rvt_pls, _, dvs = self._compute_dvs(x, lambert_legs, resos) + rvt_outs = [rvt.rotate(self._rotation_axis, self._theta) for rvt in rvt_outs] + rvt_ins[1:] = [rvt.rotate(self._rotation_axis, self._theta) for rvt in rvt_ins[1:]] + rvt_pls = [rvt.rotate(self._rotation_axis, self._theta) for rvt in rvt_pls] + + date = [epoch(rvt_pl._t * SEC2DAY) for rvt_pl in rvt_pls] + ep = [rvt_pl._t * SEC2DAY for rvt_pl in rvt_pls] + b_legs = [[rvt_out._r, rvt_out._v] for rvt_out in rvt_outs] + Vinfx, Vinfy, Vinfz = [ + a - b for a, b in zip(b_legs[0][1], self._seq[0].eph(ep[0])[1]) + ] + common_mu = rvt_outs[0]._mu + + lambert_indices = [lam.best_i for lam in lambert_legs] + + transfer_ang = _angle(rvt_outs[0]._r, rvt_outs[1]._r) + + print("Multiple Gravity Assist (MGA) + Resonance problem: ") + print("Planet sequence: ", [pl.name for pl in self._seq]) + + print("Departure: ", self._seq[0].name) + print("\tDate: ", date[0]) + print("\tEpoch: ", ep[0], " [mjd2000]") + print("\tSpacecraft velocity: ", b_legs[0][1], "[m/s]") + + print("\tr:", str(rvt_outs[0]._r)) + print("\tOutgoing V:", str(rvt_outs[0]._v)) + + print("\tLaunch velocity: ", [Vinfx, Vinfy, Vinfz], "[m/s]") + _, _, transfer_i, _, _, _ = ic2par(*(b_legs[0]), common_mu) + print("\tTransfer Angle: ", np.degrees(transfer_ang), "deg") + print("\tOutgoing Inclination:", transfer_i * RAD2DEG, "[deg]") + print("\tNumber of Revolutions:", int((lambert_indices[0] + 1) / 2)) + print("\tLambert Index:", int(lambert_indices[0])) + + lambert_i = 0 + reso_i = 0 + for i in range(1, len(self._seq) - 1): + pl = self._seq[i] + e = ep[i] + d = date[i] + dv = dvs[i] + leg = b_legs[i] + rtv_in = rvt_ins[i] + rtv_out = rvt_outs[i] + rtv_pl = rvt_pls[i] + vr_in = [a - b for a, b in zip(rtv_in._v, rtv_pl._v)] + vr_out = [a - b for a, b in zip(rtv_out._v, rtv_pl._v)] + v_inf = np.linalg.norm(vr_out) + deflection = _angle(vr_in, vr_out) + transfer_ang = _angle(rtv_out._r, rvt_outs[i + 1]._r) if i < len(self._seq) - 2 else 0 + print("Fly-by: ", pl.name) + print("\tDate: ", d) + print("\tEpoch: ", e, " [mjd2000]") + print("\tDV: ", dv, "[m/s]") + print("\tV_inf: ", v_inf, "[m/s]") + print("\tr:", str(rtv_out._r)) + print("\tOutgoing V:", str(rtv_out._v)) + print("\tOutgoing Vr:", str(vr_out)) + print("\tTransfer Angle: ", np.degrees(transfer_ang), "deg") + print("\tGA deflection: ", np.degrees(deflection), "deg") + eph = [rotate_vector(v, self._rotation_axis, self._theta) for v in pl.eph(e)] + if i < len(self._seq) - 1: + assert np.linalg.norm([a - b for a, b in zip(leg[0], eph[0])]) < 0.01 + _, _, transfer_i, _, _, _ = ic2par(eph[0], leg[1], common_mu) + print("\tOutgoing Inclination:", transfer_i * RAD2DEG, "[deg]") + if pl != self._seq[i - 1]: # no lamberts for resonances + print("\tLambert Index:", str(lambert_indices[lambert_i])) + lambert_i += 1 + else: # resonance at Venus + print("\tResonance:", str(resos[reso_i]._resonance)) + print("\tResonance time error:", str(resos[reso_i]._timing_error) + " sec") + reso_i += 1 + + print("Final Fly-by: ", self._seq[-1].name) + print("\tEpoch: ", ep[-1], " [mjd2000]") + print("\tSpacecraft velocity: ", rvt_outs[-1]._v, "[m/s]") + print("\tBeta: ", x[-1]) + print("\tr_p: ", self._seq[-1].radius + self._safe_distance) + + print("Resulting Solar orbit:") + a, e, i, _, _, _ = rvt_outs[-1].kepler() + print("Perihelion: ", (a * (1 - e)) / AU, " AU") + print("Aphelion: ", (a * (1 + e)) / AU, " AU") + print("Inclination: ", i * RAD2DEG, " degrees") + T = [SEC2DAY * (rvt_outs[i + 1]._t - rvt_outs[i]._t) for i in range(len(rvt_outs) - 1)] + print("Time of flights: ", T, "[days]") + + def plot(self, x, axes=None, units=AU, N=60): + import matplotlib as mpl + import matplotlib.pyplot as plt + from pykep.orbit_plots import plot_planet + + rvt_outs, rvt_ins, rvt_pls, _, _ = self._compute_dvs(x) + rvt_outs = [rvt.rotate(self._rotation_axis, self._theta) for rvt in rvt_outs] + rvt_ins[1:] = [rvt.rotate(self._rotation_axis, self._theta) for rvt in rvt_ins[1:]] + rvt_pls = [rvt.rotate(self._rotation_axis, self._theta) for rvt in rvt_pls] + + ep = [epoch(rvt_pl._t * SEC2DAY) for rvt_pl in rvt_pls] + + # Creating the axes if necessary + if axes is None: + mpl.rcParams["legend.fontsize"] = 10 + fig = plt.figure() + axes = fig.gca(projection="3d") + + plt.xlim([-1, 1]) + # planets + for pl, e in zip(self._seq, ep): + plot_planet( + pl, e, units=units, legend=True, color=(0.7, 0.7, 1), axes=axes + ) + + # lamberts and resonances + for i in range(0, len(self._seq) - 1): + pl = self._seq[i] + # stay at planet: it is a resonance colored black + is_reso = pl == self._seq[i + 1] + rvt_out = rvt_outs[i] + tof = rvt_ins[i + 1]._t - rvt_out._t + rvt_out.plot(tof, + units=units, + N=4 * N, + color="k" if is_reso else "r", + axes=axes) + return axes + + def eph(self, rvts, t): + for i in range(0, len(rvts)): + orb = rvts[i] + if i == len(rvts) - 1 or rvts[i + 1]._t > t: + tof = t - orb._t + orb = orb.propagate_lagrangian(tof) + return orb._r, orb._v + + def plot_distance_and_flybys(self, x, axes=None, N=1200, extension=0): + import matplotlib.pyplot as plt + rvt_outs, rvt_ins, rvt_pls, _, _ = self._compute_dvs(x) + ep = [rvt_pl._t * SEC2DAY for rvt_pl in rvt_pls] + T = [SEC2DAY * (rvt_ins[i + 1]._t - rvt_outs[i]._t) for i in range(len(rvt_outs))] + timeframe = np.linspace(0, sum(T) + extension, N) + earth = self._seq[0] + venus = self._seq[-1] + + distances = [] + edistances = [] + vdistances = [] + + for day in timeframe: + t = x[0] + day + pos, _ = self.eph(rvt_outs, t * DAY2SEC) + epos, _ = earth.eph(t) + vpos, _ = venus.eph(t) + distances.append(np.linalg.norm(pos) / AU) + edistances.append(np.linalg.norm(epos) / AU) + vdistances.append(np.linalg.norm(vpos) / AU) + + fl_times = list() + fl_distances = list() + for pl, t in zip(self._seq, ep): + fl_times.append(t - x[0]) + pos, _ = pl.eph(t) + fl_distances.append(np.linalg.norm(pos) / AU) + + if axes is None: + _, axes = plt.subplots() + plt.ylim([0, 1.3]) + + axes.plot(list(timeframe), distances, label="Solar Orbiter") + axes.plot(list(timeframe), edistances, label="Earth") + axes.plot(list(timeframe), vdistances, label="Venus") + plt.scatter(fl_times, fl_distances, marker="o", color="r") + axes.set_xlabel("Days") + axes.set_ylabel("AU") + axes.set_title("Distance to Sun") + axes.legend() + return axes + + +# propagate rvt_outs, rvt_ins, rvt_pls, dvs using MGA / Lambert +def _dv_mga(pl1, pl2, tof, max_revs, rvt_outs, rvt_ins, rvt_pls, dvs, lps=None): + rvt_pl = rvt_pls[-1] # current planet + v_in = rvt_pl._v if rvt_ins[-1] is None else rvt_ins[-1]._v + rvt_pl2 = rvt_planet(pl2, rvt_pl._t + tof) + rvt_pls.append(rvt_pl2) + r = rvt_pl._r + vpl = rvt_pl._v + r2 = rvt_pl2._r + lp = lambert_problem(r, r2, tof, rvt_pl._mu, False, max_revs) + lp = lambert_problem_multirev_ga(v_in, lp, pl1, vpl) + if not lps is None: + lps.append(lp) + v_out = lp.get_v1()[0] + rvt_out = rvt(r, v_out, rvt_pl._t, rvt_pl._mu) + rvt_outs.append(rvt_out) + rvt_in = rvt(r2, lp.get_v2()[0], rvt_pl._t + tof, rvt_pl._mu) + rvt_ins.append(rvt_in) + vr_in = [a - b for a, b in zip(v_in, vpl)] + vr_out = [a - b for a, b in zip(v_out, vpl)] + dv = fb_vel(vr_in, vr_out, pl1) + dvs.append(dv) + + +# propagate rvt_outs, rvt_ins, rvt_pls, dvs using resonances +def compute_resonance(pl, resonances, beta, safe_distance, used_resos, reso_dts, rvt_outs, + rvt_ins, rvt_pls, dvs, resos=None): + rvt_in = rvt_ins[-1] # current spaceship + rvt_pl = rvt_pls[-1] # current planet + reso = resonance(pl, rvt_in, rvt_pl, resonances) + reso_dt, used_reso = reso.select_resonance(beta, safe_distance) + if not resos is None: + resos.append(reso) + used_resos.append(used_reso) + reso_dts.append(reso_dt) + rvt_outs.append(reso._rvt_out) + tof = reso.tof() + time2 = reso._time + tof + rvt_pl2 = rvt_planet(pl, time2) + rvt_pls.append(rvt_pl2) + rvt_in2 = rvt(rvt_pl2._r, reso._rvt_out._v, time2, rvt_pl2._mu) # its a resonance, we arrive with same v as we started + rvt_ins.append(rvt_in2) + dvs.append(0) # # its a resonance, we don't need an impulse + + +def _angle(v1, v2): + ca = np.dot(v1, v2) / np.linalg.norm(v1) / np.linalg.norm(v2) # -> cosine of the angle + return np.arccos(np.clip(ca, -1, 1)) + + +solo_mgar = _solo_mgar_udp(max_revs=2) +solo_mgar_unconstrained = _solo_mgar_udp(max_revs=2, use_constraints=False) +solo_mgar_multi_objective = _solo_mgar_udp(max_revs=2, multi_objective=True) +