1from pykep.core import epoch, DAY2SEC, MU_SUN, lambert_problem, propagate_lagrangian, fb_prop, AU, epoch 2from pykep.planet import jpl_lp 3from pykep.trajopt._lambert import lambert_problem_multirev 4from math import pi, cos, sin, acos, log, sqrt 5import numpy as np 6 7# Avoiding scipy dependency 8def norm(x): 9 return sqrt(sum([it * it for it in x])) 10 11class mga_1dsm: 12 r""" 13 This class transcribes a Multiple Gravity Assist trajectory with one deep space maneuver per leg (MGA-1DSM) into an optimisation problem. 14 It may be used as a User Defined Problem (UDP) for the pygmo (http://esa.github.io/pygmo/) optimisation suite. 15 16 - Izzo, Dario. "Global optimization and space pruning for spacecraft trajectory design." Spacecraft Trajectory Optimization 1 (2010): 178-200. 17 18 The decision vector (chromosome) is:: 19 20 direct encoding: [t0] + [u, v, Vinf, eta1, T1] + [beta, rp/rV, eta2, T2] + ... 21 alpha encoding: [t0] + [u, v, Vinf, eta1, a1] + [beta, rp/rV, eta2, a2] + ... + [T] 22 eta encoding: [t0] + [u, v, Vinf, eta1, n1] + [beta, rp/rV, eta2, n2] + ... 23 24 where t0 is a mjd2000, Vinf is in km/s, T in days, beta in radians and the rest non dimensional. 25 26 .. note:: 27 28 The time of flights of a MGA-1DSM trajectory (and in general) can be encoded in different ways. 29 When they are directly present in the decision vector, we have a *direct* encoding. This is the most 'evolvable' encoding 30 but also the one that requires the most problem knowledge (e.g. to define the bounds on each leg) and is not 31 very flexible in dealing with constraints on the total time of flight. The *alpha* and *eta* encodings, instead, allow 32 to only specify bounds on the time of flight of the entire trajectory, and not on the single legs: a property that is attractive 33 for multi-objective optimization, for example. 34 35 In the *alpha* encoding each leg time-of-flight is decoded as follows, T_i = T log(alpha_i) / \sum_n(log(alpha_n)). 36 In the *eta* encoding each leg time-of-flight is decoded as follows, T_i = (tof_max - \sum_0^(i-1)(T_j)) * eta_i 37 38 The chromosome dimension for the direct and eta encoding is the same, while the alpha encoding requires one more gene. 39 40 .. note:: 41 42 The resulting problem is box-bounded (unconstrained). 43 """ 44 45 def __init__(self, 46 seq = [jpl_lp('earth'), jpl_lp('venus'), jpl_lp('earth')], 47 t0 = [epoch(0), epoch(1000)], 48 tof = [[10, 300], [10, 300]], 49 vinf = [0.5, 2.5], 50 add_vinf_dep = False, 51 add_vinf_arr = True, 52 tof_encoding = 'direct', 53 multi_objective = False, 54 orbit_insertion = False, 55 e_target = None, 56 rp_target = None, 57 eta_lb = 0.1, 58 eta_ub = 0.9, 59 rp_ub = 30, 60 max_revs = 0 61 ): 62 """ 63 pykep.trajopt.mga_1dsm(seq = [jpl_lp('earth'), jpl_lp('venus'), jpl_lp('earth')], t0 = [epoch(0),epoch(1000)], tof = [1.0,5.0], vinf = [0.5, 2.5], multi_objective = False, add_vinf_dep = False, add_vinf_arr=True) 64 65 - seq (``list`` of ``pykep.planet``): the encounter sequence (including the starting launch) 66 - t0 (``list`` of ``pykep.epoch`` or ``floats``): the launch window (in mjd2000 if floats) 67 - tof (``list`` or ``float``): bounds on the time of flight (days). If *tof_encoding* is 'direct', this contains a list 68 of 2D lists defining the upper and lower bounds on each leg. If *tof_encoding* is 'alpha', 69 this contains a list of two floats containing the lower and upper bounds on the time-of-flight. If *tof_encoding* 70 is 'eta' tof is a float defining the upper bound on the time-of-flight 71 - vinf (``list``): the minimum and maximum allowed initial hyperbolic velocity (at launch), in km/sec 72 - add_vinf_dep (``bool``): when True the computed Dv includes the initial hyperbolic velocity (at launch) 73 - add_vinf_arr (``bool``): when True the computed Dv includes the final hyperbolic velocity (at the last planet) 74 - tof_encoding (``str``): one of 'direct', 'alpha' or 'eta'. Selects the encoding for the time of flights 75 - multi_objective (``bool``): when True constructs a multiobjective problem (dv, T) 76 - orbit_insertion (``bool``): when True the arrival dv is computed as that required to acquire a target orbit defined by e_target and rp_target 77 - e_target (``float``): if orbit_insertion is True this defines the target orbit eccentricity around the final planet 78 - rp_target (``float``): if orbit_insertion is True this defines the target orbit pericenter around the final planet (in m) 79 - max_revs (``int``): maximal number of revolutions for lambert transfer 80 """ 81 82 # Sanity checks 83 # 1 - Planets need to have the same mu_central_body 84 if ([r.mu_central_body for r in seq].count(seq[0].mu_central_body) != len(seq)): 85 raise ValueError( 86 'All planets in the sequence need to have identical mu_central_body') 87 # 2 - tof encoding needs to be one of 'alpha', 'eta', 'direct' 88 if tof_encoding not in ['alpha', 'eta', 'direct']: 89 raise TypeError( 90 'tof encoding must be one of \'alpha\', \'eta\', \'direct\'') 91 # 3 - tof is expected to have different content depending on the tof_encoding 92 if tof_encoding == 'direct': 93 if np.shape(np.array(tof)) != (len(seq) - 1, 2): 94 raise TypeError( 95 'tof_encoding is ' + tof_encoding + ' and tof must be a list of two dimensional lists and with length equal to the number of legs') 96 if tof_encoding == 'alpha': 97 if np.shape(np.array(tof)) != (2,): 98 raise TypeError( 99 'tof_encoding is ' + tof_encoding + ' and tof must be a list of two floats') 100 if tof_encoding == 'eta': 101 if np.shape(np.array(tof)) != (): 102 raise TypeError( 103 'tof_encoding is ' + tof_encoding + ' and tof must be a float') 104 # 4 - Check launch window t0. If defined in terms of floats transform into epochs 105 if len(t0) != 2: 106 raise TypeError( 107 't0 is ' + t0 + ' while should be a list of two floats or epochs') 108 if type(t0[0]) is not epoch: 109 t0[0] = epoch(t0[0]) 110 if type(t0[1]) is not epoch: 111 t0[1] = epoch(t0[1]) 112 # 5 - Check that if orbit insertion is selected e_target and r_p are 113 # defined 114 if orbit_insertion: 115 if rp_target is None: 116 raise ValueError( 117 'The rp_target needs to be specified when orbit insertion is selected') 118 if e_target is None: 119 raise ValueError( 120 'The e_target needs to be specified when orbit insertion is selected') 121 if add_vinf_arr is False: 122 raise ValueError( 123 'When orbit insertion is selected, the add_vinf_arr must be True') 124 125 self._seq = seq 126 self._t0 = t0 127 self._tof = tof 128 self._vinf = vinf 129 self._add_vinf_dep = add_vinf_dep 130 self._add_vinf_arr = add_vinf_arr 131 self._tof_encoding = tof_encoding 132 self._multi_objective = multi_objective 133 self._orbit_insertion = orbit_insertion 134 self._e_target = e_target 135 self._rp_target = rp_target 136 self._eta_lb = eta_lb 137 self._eta_ub = eta_ub 138 self._rp_ub = rp_ub 139 140 self.n_legs = len(seq) - 1 141 self.common_mu = seq[0].mu_central_body 142 self.max_revs = max_revs 143 144 def get_nobj(self): 145 return self._multi_objective + 1 146 147 def get_bounds(self): 148 t0 = self._t0 149 tof = self._tof 150 vinf = self._vinf 151 seq = self._seq 152 # Base for all possiblities (eta encoding) 153 lb = [t0[0].mjd2000] + [0.0, 0.0, vinf[0] * 1000, self._eta_lb, 154 1e-3] + [-2 * pi, np.nan, self._eta_lb, 1e-3] * (self.n_legs - 1) 155 ub = [t0[1].mjd2000] + [1.0, 1.0, vinf[1] * 1000, self._eta_ub, 156 1.0 - 1e-3] + [2 * pi, self._rp_ub, self._eta_ub, 1.0 - 1e-3] * (self.n_legs - 1) 157 # Distinguishing among cases (only direct and alpha) 158 if self._tof_encoding == 'alpha': 159 lb = lb + [tof[0]] 160 ub = ub + [tof[1]] 161 elif self._tof_encoding == 'direct': 162 for i in range(self.n_legs): 163 lb[5 + 4 * i] = tof[i][0] 164 ub[5 + 4 * i] = tof[i][1] 165 166 # Setting the minimum rp/rP using the planet safe radius 167 for i, pl in enumerate(seq[1:-1]): 168 lb[7 + 4 * i] = pl.safe_radius / pl.radius 169 return (lb, ub) 170 171 def _decode_times_and_vinf(self, x): 172 # 1 - we decode the times of flight 173 if self._tof_encoding == 'alpha': 174 # decision vector is [t0] + [u, v, Vinf, eta1, a1] + [beta, rp/rV, eta2, a2] + ... + [T] 175 T = list([0] * (self.n_legs)) 176 for i in range(len(T)): 177 T[i] = -log(x[5 + 4 * i]) 178 alpha_sum = sum(T) 179 retval_T = [x[-1] * time / alpha_sum for time in T] 180 elif self._tof_encoding == 'direct': 181 # decision vector is [t0] + [u, v, Vinf, eta1, T1] + [beta, rp/rV, eta2, T2] + ... 182 retval_T = x[5::4] 183 elif self._tof_encoding == 'eta': 184 # decision vector is [t0] + [u, v, Vinf, eta1, n1] + [beta, rp/rV, eta2, n2] + ... 185 dt = self._tof 186 T = [0] * self.n_legs 187 for i in range(self.n_legs): 188 T[i] = (dt - sum(T[:i])) * x[5 + 4*i] 189 retval_T = T 190 191 # 2 - we decode the hyperbolic velocity at departure 192 theta = 2 * pi * x[1] 193 phi = acos(2 * x[2] - 1) - pi / 2 194 195 Vinfx = x[3] * cos(phi) * cos(theta) 196 Vinfy = x[3] * cos(phi) * sin(theta) 197 Vinfz = x[3] * sin(phi) 198 199 return (retval_T, Vinfx, Vinfy, Vinfz) 200 201 # Objective function 202 def fitness(self, x): 203 # 1 - we 'decode' the chromosome recording the various times of flight 204 # (days) in the list T and the cartesian components of vinf 205 T, Vinfx, Vinfy, Vinfz = self._decode_times_and_vinf(x) 206 207 # 2 - We compute the epochs and ephemerides of the planetary encounters 208 t_P = list([None] * (self.n_legs + 1)) 209 r_P = list([None] * (self.n_legs + 1)) 210 v_P = list([None] * (self.n_legs + 1)) 211 DV = list([0.0] * (self.n_legs + 1)) 212 for i in range(len(self._seq)): 213 t_P[i] = epoch(x[0] + sum(T[0:i])) 214 r_P[i], v_P[i] = self._seq[i].eph(t_P[i]) 215 216 # 3 - We start with the first leg 217 v0 = [a + b for a, b in zip(v_P[0], [Vinfx, Vinfy, Vinfz])] 218 r, v = propagate_lagrangian( 219 r_P[0], v0, x[4] * T[0] * DAY2SEC, self.common_mu) 220 221 # Lambert arc to reach seq[1] 222 dt = (1 - x[4]) * T[0] * DAY2SEC 223 l = lambert_problem_multirev(v, lambert_problem( 224 r, r_P[1], dt, self.common_mu, cw=False, max_revs=self.max_revs)) 225 v_end_l = l.get_v2()[0] 226 v_beg_l = l.get_v1()[0] 227 228 # First DSM occuring at time nu1*T1 229 DV[0] = norm([a - b for a, b in zip(v_beg_l, v)]) 230 231 # 4 - And we proceed with each successive leg 232 for i in range(1, self.n_legs): 233 # Fly-by 234 v_out = fb_prop(v_end_l, v_P[i], x[ 235 7 + (i - 1) * 4] * self._seq[i].radius, x[6 + (i - 1) * 4], self._seq[i].mu_self) 236 # s/c propagation before the DSM 237 r, v = propagate_lagrangian( 238 r_P[i], v_out, x[8 + (i - 1) * 4] * T[i] * DAY2SEC, self.common_mu) 239 # Lambert arc to reach Earth during (1-nu2)*T2 (second segment) 240 dt = (1 - x[8 + (i - 1) * 4]) * T[i] * DAY2SEC 241 l = lambert_problem_multirev(v, lambert_problem(r, r_P[i + 1], dt, 242 self.common_mu, cw=False, max_revs=self.max_revs)) 243 v_end_l = l.get_v2()[0] 244 v_beg_l = l.get_v1()[0] 245 # DSM occuring at time nu2*T2 246 DV[i] = norm([a - b for a, b in zip(v_beg_l, v)]) 247 248 # Last Delta-v 249 if self._add_vinf_arr: 250 DV[-1] = norm([a - b for a, b in zip(v_end_l, v_P[-1])]) 251 if self._orbit_insertion: 252 # In this case we compute the insertion DV as a single pericenter 253 # burn 254 DVper = np.sqrt(DV[-1] * DV[-1] + 2 * 255 self._seq[-1].mu_self / self._rp_target) 256 DVper2 = np.sqrt(2 * self._seq[-1].mu_self / self._rp_target - 257 self._seq[-1].mu_self / self._rp_target * (1. - self._e_target)) 258 DV[-1] = np.abs(DVper - DVper2) 259 260 if self._add_vinf_dep: 261 DV[0] += x[3] 262 263 if not self._multi_objective: 264 return (sum(DV),) 265 else: 266 return (sum(DV), sum(T)) 267 268 def pretty(self, x): 269 """ 270 prob.plot(x) 271 272 - x: encoded trajectory 273 274 Prints human readable information on the trajectory represented by the decision vector x 275 276 Example:: 277 278 print(prob.pretty(x)) 279 """ 280 # 1 - we 'decode' the chromosome recording the various times of flight 281 # (days) in the list T and the cartesian components of vinf 282 T, Vinfx, Vinfy, Vinfz = self._decode_times_and_vinf(x) 283 284 # 2 - We compute the epochs and ephemerides of the planetary encounters 285 t_P = list([None] * (self.n_legs + 1)) 286 r_P = list([None] * (self.n_legs + 1)) 287 v_P = list([None] * (self.n_legs + 1)) 288 DV = list([0.0] * (self.n_legs + 1)) 289 for i in range(len(self._seq)): 290 t_P[i] = epoch(x[0] + sum(T[0:i])) 291 r_P[i], v_P[i] = self._seq[i].eph(t_P[i]) 292 293 # 3 - We start with the first leg 294 print("First Leg: " + self._seq[0].name + " to " + self._seq[1].name) 295 print("Departure: " + str(t_P[0]) + 296 " (" + str(t_P[0].mjd2000) + " mjd2000) ") 297 print("Duration: " + str(T[0]) + "days") 298 print("VINF: " + str(x[3] / 1000) + " km/sec") 299 300 v0 = [a + b for a, b in zip(v_P[0], [Vinfx, Vinfy, Vinfz])] 301 r, v = propagate_lagrangian( 302 r_P[0], v0, x[4] * T[0] * DAY2SEC, self.common_mu) 303 304 print("DSM after " + str(x[4] * T[0]) + " days") 305 306 # Lambert arc to reach seq[1] 307 dt = (1 - x[4]) * T[0] * DAY2SEC 308 l = lambert_problem_multirev(v, lambert_problem( 309 r, r_P[1], dt, self.common_mu, cw=False, max_revs=self.max_revs)) 310 v_end_l = l.get_v2()[0] 311 v_beg_l = l.get_v1()[0] 312 313 # First DSM occuring at time nu1*T1 314 DV[0] = norm([a - b for a, b in zip(v_beg_l, v)]) 315 print("DSM magnitude: " + str(DV[0]) + "m/s") 316 317 # 4 - And we proceed with each successive leg 318 for i in range(1, self.n_legs): 319 print("\nleg no. " + str(i + 1) + ": " + 320 self._seq[i].name + " to " + self._seq[i + 1].name) 321 print("Duration: " + str(T[i]) + "days") 322 # Fly-by 323 v_out = fb_prop(v_end_l, v_P[i], x[ 324 7 + (i - 1) * 4] * self._seq[i].radius, x[6 + (i - 1) * 4], self._seq[i].mu_self) 325 print( 326 "Fly-by epoch: " + str(t_P[i]) + " (" + str(t_P[i].mjd2000) + " mjd2000) ") 327 print( 328 "Fly-by radius: " + str(x[7 + (i - 1) * 4]) + " planetary radii") 329 # s/c propagation before the DSM 330 r, v = propagate_lagrangian( 331 r_P[i], v_out, x[8 + (i - 1) * 4] * T[i] * DAY2SEC, self.common_mu) 332 print("DSM after " + str(x[8 + (i - 1) * 4] * T[i]) + " days") 333 # Lambert arc to reach Earth during (1-nu2)*T2 (second segment) 334 dt = (1 - x[8 + (i - 1) * 4]) * T[i] * DAY2SEC 335 l = lambert_problem_multirev(v, lambert_problem(r, r_P[i + 1], dt, 336 self.common_mu, cw=False, max_revs=self.max_revs)) 337 v_end_l = l.get_v2()[0] 338 v_beg_l = l.get_v1()[0] 339 # DSM occuring at time nu2*T2 340 DV[i] = norm([a - b for a, b in zip(v_beg_l, v)]) 341 print("DSM magnitude: " + str(DV[i]) + "m/s") 342 343 # Last Delta-v 344 print("\nArrival at " + self._seq[-1].name) 345 DV[-1] = norm([a - b for a, b in zip(v_end_l, v_P[-1])]) 346 print( 347 "Arrival epoch: " + str(t_P[-1]) + " (" + str(t_P[-1].mjd2000) + " mjd2000) ") 348 print("Arrival Vinf: " + str(DV[-1]) + "m/s") 349 if self._orbit_insertion: 350 # In this case we compute the insertion DV as a single pericenter 351 # burn 352 DVper = np.sqrt(DV[-1] * DV[-1] + 2 * 353 self._seq[-1].mu_self / self._rp_target) 354 DVper2 = np.sqrt(2 * self._seq[-1].mu_self / self._rp_target - 355 self._seq[-1].mu_self / self._rp_target * (1. - self._e_target)) 356 DVinsertion = np.abs(DVper - DVper2) 357 print("Insertion DV: " + str(DVinsertion) + "m/s") 358 359 print("Total mission time: " + str(sum(T) / 365.25) + " years (" + str(sum(T)) + " days)") 360 361 # Plot of the trajectory 362 def plot(self, x, ax = None): 363 """ 364 ax = prob.plot(x, ax=None) 365 366 - x: encoded trajectory 367 - ax: matplotlib axis where to plot. If None figure and axis will be created 368 - [out] ax: matplotlib axis where to plot 369 370 Plots the trajectory represented by a decision vector x on the 3d axis ax 371 372 Example:: 373 374 ax = prob.plot(x) 375 """ 376 import matplotlib as mpl 377 from mpl_toolkits.mplot3d import Axes3D 378 import matplotlib.pyplot as plt 379 from pykep.orbit_plots import plot_planet, plot_lambert, plot_kepler 380 381 if ax is None: 382 mpl.rcParams['legend.fontsize'] = 10 383 fig = plt.figure() 384 axis = fig.gca(projection='3d') 385 else: 386 axis = ax 387 388 axis.scatter(0, 0, 0, color='y') 389 390 # 1 - we 'decode' the chromosome recording the various times of flight 391 # (days) in the list T and the cartesian components of vinf 392 T, Vinfx, Vinfy, Vinfz = self._decode_times_and_vinf(x) 393 394 # 2 - We compute the epochs and ephemerides of the planetary encounters 395 t_P = list([None] * (self.n_legs + 1)) 396 r_P = list([None] * (self.n_legs + 1)) 397 v_P = list([None] * (self.n_legs + 1)) 398 DV = list([None] * (self.n_legs + 1)) 399 400 for i, planet in enumerate(self._seq): 401 t_P[i] = epoch(x[0] + sum(T[0:i])) 402 r_P[i], v_P[i] = planet.eph(t_P[i]) 403 plot_planet(planet, t0=t_P[i], color=( 404 0.8, 0.6, 0.8), legend=True, units=AU, axes=axis, N=150) 405 406 # 3 - We start with the first leg 407 v0 = [a + b for a, b in zip(v_P[0], [Vinfx, Vinfy, Vinfz])] 408 r, v = propagate_lagrangian( 409 r_P[0], v0, x[4] * T[0] * DAY2SEC, self.common_mu) 410 411 plot_kepler(r_P[0], v0, x[4] * T[0] * DAY2SEC, self.common_mu, 412 N=100, color='b', units=AU, axes=axis) 413 414 # Lambert arc to reach seq[1] 415 dt = (1 - x[4]) * T[0] * DAY2SEC 416 417 l = lambert_problem_multirev(v, lambert_problem( 418 r, r_P[1], dt, self.common_mu, cw=False, max_revs=self.max_revs)) 419 420 plot_lambert(l, sol=0, color='r', units=AU, axes=axis) 421 v_end_l = l.get_v2()[0] 422 v_beg_l = l.get_v1()[0] 423 424 # First DSM occuring at time nu1*T1 425 DV[0] = norm([a - b for a, b in zip(v_beg_l, v)]) 426 427 # 4 - And we proceed with each successive leg 428 for i in range(1, self.n_legs): 429 # Fly-by 430 v_out = fb_prop(v_end_l, v_P[i], x[ 431 7 + (i - 1) * 4] * self._seq[i].radius, x[6 + (i - 1) * 4], self._seq[i].mu_self) 432 # s/c propagation before the DSM 433 r, v = propagate_lagrangian( 434 r_P[i], v_out, x[8 + (i - 1) * 4] * T[i] * DAY2SEC, self.common_mu) 435 plot_kepler(r_P[i], v_out, x[8 + (i - 1) * 4] * T[i] * DAY2SEC, 436 self.common_mu, N=100, color='b', units=AU, axes=axis) 437 # Lambert arc to reach Earth during (1-nu2)*T2 (second segment) 438 dt = (1 - x[8 + (i - 1) * 4]) * T[i] * DAY2SEC 439 440 l = lambert_problem_multirev(v, lambert_problem(r, r_P[i + 1], dt, 441 self.common_mu, cw=False, max_revs=self.max_revs)) 442 443 plot_lambert(l, sol=0, color='r', legend=False, 444 units=AU, N=1000, axes=axis) 445 446 v_end_l = l.get_v2()[0] 447 v_beg_l = l.get_v1()[0] 448 # DSM occuring at time nu2*T2 449 DV[i] = norm([a - b for a, b in zip(v_beg_l, v)]) 450 plt.show() 451 return axis 452 453 def get_extra_info(self): 454 return ("\n\t Sequence: " + [pl.name for pl in self._seq].__repr__() + 455 "\n\t Add launcher vinf to the objective?: " + self._add_vinf_dep.__repr__() + 456 "\n\t Add final vinf to the objective?: " + self._add_vinf_arr.__repr__()) 457