1from pykep.trajopt import mga_1dsm, launchers 2from pykep.planet import jpl_lp 3from pykep import epoch_from_string 4 5import numpy as np 6from numpy.linalg import norm 7from math import log, acos, cos, sin, asin, exp 8from copy import deepcopy 9 10class _tandem_udp(mga_1dsm): 11 """ 12 This class represents a rendezvous mission to Saturn modelled as an MGA-1DSM transfer. Mission parameters are 13 inspired to the TandEM mission. A launcher model (i.e. Atlas 501) is also used, so that the final mass delivered 14 at Saturn is the main objective of this optimization problem. 15 16 The problem draws inspiration from the work performed in April 2008 by the 17 European Space Agency working group on mission analysis on the mission named TandEM. TandEM is an interplanetary 18 mission aimed at reaching Titan and Enceladus (two moons of Saturn). 19 20 .. note:: 21 22 The Titan and Enceladus Mission (TandEM), an ambitious scientific mission to study the Saturnian system 23 with particular emphasis on the moons Titan and Enceladus, was selected in October 2007 as a candidate mission 24 within the ESA Cosmic Vision plan. In February 2009, TandEM exited the Cosmic Vision programme when ESA 25 and NASA chose EJSM-Laplace as the L-class outer Solar System mission candidate. 26 27 .. note:: 28 29 A significantly similar version of this problem was part of the no longer maintained GTOP database, 30 https://www.esa.int/gsp/ACT/projects/gtop/gtop.html. The exact definition is, though, different and results 31 cannot thus not be compared to those posted in GTOP. 32 """ 33 def __init__(self, prob_id = 1, constrained = True): 34 """ 35 The TandEM problem of the trajectory gym consists in 48 different instances varying in fly-by sequence and 36 the presence of a time constraint. 37 38 Args: 39 - prob_id (``int``): The problem id defines the fly-by sequence. 40 - constrained (``bool``): Activates the constraint on the time of flight 41 (fitness will thus return two numbers, the objective function and the inequality constraint violation) 42 """ 43 # Redefining the planets as to change their safe radius 44 earth = jpl_lp('earth') 45 earth.safe_radius = 1.05 46 # We need the Earth eph in the fitness 47 venus = jpl_lp('venus') 48 venus.safe_radius = 1.05 49 mars = jpl_lp('mars') 50 mars.safe_radius = 1.05 51 jupiter = jpl_lp('jupiter') 52 jupiter.safe_radius = 1.7 53 saturn = jpl_lp('saturn') 54 55 # Defining the different sequences 56 seq_tandem = [] 57 seq_tandem.append([earth, venus, venus, venus, saturn]) 58 seq_tandem.append([earth, venus, venus, earth, saturn]) 59 seq_tandem.append([earth, venus, venus, mars, saturn]) 60 seq_tandem.append([earth, venus, venus, jupiter, saturn]) 61 62 seq_tandem.append([earth, venus, earth, venus, saturn]) 63 seq_tandem.append([earth, venus, earth, earth, saturn]) 64 seq_tandem.append([earth, venus, earth, mars, saturn]) 65 seq_tandem.append([earth, venus, earth, jupiter, saturn]) 66 67 seq_tandem.append([earth, venus, mars, venus, saturn]) 68 seq_tandem.append([earth, venus, mars, earth, saturn]) 69 seq_tandem.append([earth, venus, mars, mars, saturn]) 70 seq_tandem.append([earth, venus, mars, jupiter, saturn]) 71 72 seq_tandem.append([earth, earth, venus, venus, saturn]) 73 seq_tandem.append([earth, earth, venus, earth, saturn]) 74 seq_tandem.append([earth, earth, venus, mars, saturn]) 75 seq_tandem.append([earth, earth, venus, jupiter, saturn]) 76 77 seq_tandem.append([earth, earth, earth, venus, saturn]) 78 seq_tandem.append([earth, earth, earth, earth, saturn]) 79 seq_tandem.append([earth, earth, earth, mars, saturn]) 80 seq_tandem.append([earth, earth, earth, jupiter, saturn]) 81 82 seq_tandem.append([earth, earth, mars, venus, saturn]) 83 seq_tandem.append([earth, earth, mars, earth, saturn]) 84 seq_tandem.append([earth, earth, mars, mars, saturn]) 85 seq_tandem.append([earth, earth, mars, jupiter, saturn]) 86 87 if prob_id > 24 or type(prob_id) != int or prob_id < 1: 88 raise ValueError("TandEM problem id must be an integer in [1, 24]") 89 90 super().__init__( 91 seq = seq_tandem[prob_id - 1], 92 t0 = [5475, 9132], 93 tof = [[20, 2500], [20, 2500], [20, 2500], [20, 2500]], 94 vinf = [2.5, 4.9], 95 add_vinf_dep = False, 96 add_vinf_arr = True, 97 tof_encoding = 'direct', 98 multi_objective = False, 99 orbit_insertion = True, 100 e_target = 0.98531407996358, 101 rp_target = 80330000, 102 eta_lb = 0.01, 103 eta_ub = 0.99, 104 rp_ub = 10 105 ) 106 107 self.prob_id = prob_id 108 self.constrained = constrained 109 110 def fitness(self, x): 111 T, Vinfx, Vinfy, Vinfz = self._decode_times_and_vinf(x) 112 # We transform it (only the needed component) to an equatorial system rotating along x 113 # (this is an approximation, assuming vernal equinox is roughly x and the ecliptic plane is roughly xy) 114 earth_axis_inclination = 0.409072975 115 # This is different from the GTOP tanmEM problem, I think it was bugged there as the rotation was in the wrong direction. 116 Vinfz = - Vinfy * sin(earth_axis_inclination) + Vinfz * cos(earth_axis_inclination) 117 # And we find the vinf declination (in degrees) 118 sindelta = Vinfz / x[3] 119 declination = asin(sindelta) / np.pi * 180. 120 # We now have the initial mass of the spacecraft 121 m_initial = launchers.atlas501(x[3] / 1000., declination) 122 # And we can evaluate the final mass via Tsiolkowsky 123 Isp = 312. 124 g0 = 9.80665 125 DV = super().fitness(x)[0] 126 DV = DV + 165. # losses for 3 swgbys + insertion 127 m_final = m_initial * exp(-DV / Isp / g0) 128 # Numerical guard for the exponential 129 if m_final == 0: 130 m_final = 1e-320 131 if self.constrained: 132 retval = [-log(m_final), sum(T) - 3652.5] 133 else: 134 retval = [-log(m_final), ] 135 return retval 136 137 def get_nic(self): 138 return int(self.constrained) 139 140 def get_name(self): 141 return "TandEM problem, id: " + str(self.prob_id) 142 143 def get_extra_info(self): 144 retval = "\t Sequence: " + \ 145 [pl.name for pl in self._seq].__repr__() + "\n\t Constrained: " + \ 146 str(self.constrained) 147 return retval 148 149 def pretty(self, x): 150 """ 151 prob.plot(x) 152 153 - x: encoded trajectory 154 155 Prints human readable information on the trajectory represented by the decision vector x 156 157 Example:: 158 159 print(prob.pretty(x)) 160 """ 161 super().pretty(x) 162 T, Vinfx, Vinfy, Vinfz = self._decode_times_and_vinf(x) 163 # We transform it (only the needed component) to an equatorial system rotating along x 164 # (this is an approximation, assuming vernal equinox is roughly x and the ecliptic plane is roughly xy) 165 earth_axis_inclination = 0.409072975 166 # This is different from the GTOP tanmEM problem, I think it was bugged there as the rotation was in the wrong direction. 167 Vinfz = - Vinfy * sin(earth_axis_inclination) + Vinfz * cos(earth_axis_inclination) 168 # And we find the vinf declination (in degrees) 169 sindelta = Vinfz / x[3] 170 declination = asin(sindelta) / np.pi * 180. 171 m_initial = launchers.atlas501(x[3] / 1000., declination) 172 # And we can evaluate the final mass via Tsiolkowsky 173 Isp = 312. 174 g0 = 9.80665 175 DV = super().fitness(x)[0] 176 DV = DV + 165. # losses for 3 swgbys + insertion 177 m_final = m_initial * exp(-DV / Isp / g0) 178 print("\nInitial mass:", m_initial) 179 print("Final mass:", m_final) 180 print("Declination:", declination) 181 182 def __repr__(self): 183 return "TandEM (Trajectory Optimisation Gym P12, multiple instances)" 184 185# Problem P12: TandEM mission MGA1DSM, single objective, direct encoding, time constrained 186tandem = _tandem_udp 187