1import pykep as pk 2import numpy as np 3import math 4 5mpcorbline = "K14Y00D 24.3 0.15 K1794 105.44160 34.12337 117.64264 1.73560 0.0865962 0.88781021 1.0721510 2 MPO369254 104 1 194 days 0.24 M-v 3Eh MPCALB 2803 2014 YD 20150618" 6 7 8class lt_margo: 9 """ 10 This class can be used as a User Defined Problem (UDP) in the pygmo2 software and if successfully solved, 11 represents a low-thrust interplanetary trajectory from the Earth (or from the 12 Sun-Earth L1 or L2 Lagrangian point) to a target NEO. The trajectory is modeled using the Sims-Flanagan model, 13 extended to include the Earth's gravity (assumed constant along each segment). 14 The propulsion model can be both nuclear (NEP) or solar (SEP). 15 16 This problem was developed and used during the European Space Agency interplanetary cubesat M-ARGO study. 17 18 The decision vector (chromosome) is:: 19 20 [t0, tof, mf] + [throttles1] + [throttles2] + ... 21 """ 22 23 def __init__(self, 24 target=pk.planet.mpcorb(mpcorbline), 25 n_seg=30, 26 grid_type="uniform", 27 t0=[pk.epoch(8000), pk.epoch(9000)], 28 tof=[200, 365.25 * 3], 29 m0=20.0, 30 Tmax=0.0017, 31 Isp=3000.0, 32 earth_gravity=False, 33 sep=False, 34 start="earth"): 35 """ 36 prob = lt_margo(target = pk.planet.mpcorb(mpcorbline), n_seg = 30, grid_type = "uniform", t0 = [epoch(8000), epoch(9000)], tof = [200, 365.25*3], m0 = 20.0, Tmax = 0.0017, Isp = 3000.0, earth_gravity = False, sep = False, start = "earth") 37 38 Args: 39 - target (``pykep.planet``): target planet 40 - n_seg (``int``): number of segments to use in the problem transcription (time grid) 41 - grid_type (``string``): "uniform" for uniform segments, "nonuniform" toi use a denser grid in the first part of the trajectory 42 - t0 (``list``): list of two pykep.epoch defining the bounds on the launch epoch 43 - tof (``list``): list of two floats defining the bounds on the time of flight (days) 44 - m0 (``float``): initial mass of the spacecraft (kg) 45 - Tmax (``float``): maximum thrust at 1 AU (N) 46 - Isp (``float``): engine specific impulse (s) 47 - earth_gravity (``bool``): activates the Earth gravity effect in the dynamics 48 - sep (``bool``): Activates a Solar Electric Propulsion model for the thrust - distance dependency. 49 - start(``string``): Starting point ("earth", "l1", or "l2"). 50 51 .. note:: 52 53 L1 and L2 are approximated as the points on the line connecting the Sun and the Earth at a distance of, respectively, 0.99 and 1.01 AU from the Sun. 54 55 .. note:: 56 57 If the Earth's gravity is enabled, the starting point cannot be the Earth 58 """ 59 60 # 1) Various checks 61 if start not in ["earth", "l1", "l2"]: 62 raise ValueError("start must be either 'earth', 'l1' or 'l2'") 63 if grid_type not in ["uniform", "nonuniform"]: 64 raise ValueError( 65 "grid_type must be either 'uniform' or 'nonuniform'") 66 if earth_gravity and start == "earth": 67 raise ValueError( 68 "If Earth gravity is enabled the starting point cannot be the Earth") 69 70 # 2) Class data members 71 # public: 72 self.target = target 73 # private: 74 self.__n_seg = n_seg 75 self.__grid_type = grid_type 76 self.__sc = pk.sims_flanagan.spacecraft(m0, Tmax, Isp) 77 self.__earth = pk.planet.jpl_lp('earth') 78 self.__earth_gravity = earth_gravity 79 self.__sep = sep 80 self.__start = start 81 # grid construction 82 if grid_type == "uniform": 83 grid = np.array([i / n_seg for i in range(n_seg + 1)]) 84 elif grid_type == "nonuniform": 85 grid_f = lambda x: x**2 if x < 0.5 else 0.25 + 1.5 * \ 86 (x - 0.5) # quadratic in [0,0.5], linear in [0.5,1] 87 grid = np.array([grid_f(i / n_seg) for i in range(n_seg + 1)]) 88 # index corresponding to the middle of the transfer 89 fwd_seg = int(np.searchsorted(grid, 0.5, side='right')) 90 bwd_seg = n_seg - fwd_seg 91 fwd_grid = grid[:fwd_seg + 1] 92 bwd_grid = grid[fwd_seg:] 93 self.__fwd_seg = fwd_seg 94 self.__fwd_grid = fwd_grid 95 self.__fwd_dt = np.array([(fwd_grid[i + 1] - fwd_grid[i]) 96 for i in range(fwd_seg)]) * pk.DAY2SEC 97 self.__bwd_seg = bwd_seg 98 self.__bwd_grid = bwd_grid 99 self.__bwd_dt = np.array([(bwd_grid[i + 1] - bwd_grid[i]) 100 for i in range(bwd_seg)]) * pk.DAY2SEC 101 102 # 3) Bounds 103 lb = [t0[0].mjd2000] + [tof[0]] + [0] + [-1, -1, -1] * n_seg 104 ub = [t0[1].mjd2000] + [tof[1]] + [m0] + [1, 1, 1] * n_seg 105 self.__lb = lb 106 self.__ub = ub 107 108 # Fitness function 109 def fitness(self, x): 110 obj = -x[2] 111 112 ceq = list() 113 cineq = list() 114 throttles_constraints = [] 115 mismatch_constraints = [] 116 117 throttles = [x[3 + 3 * i: 6 + 3 * i] for i in range(self.__n_seg)] 118 for t in throttles: 119 throttles_constraints.append(t[0]**2 + t[1]**2 + t[2]**2 - 1.) 120 cineq.extend(throttles_constraints) 121 122 rfwd, rbwd, vfwd, vbwd, mfwd, mbwd, _, _, _, _, dfwd, dbwd = self._propagate( 123 x) 124 125 if self.__earth_gravity: 126 flyby_constraints = [] 127 for d in dfwd + dbwd: 128 d2 = d[0]**2 + d[1]**2 + d[2]**2 129 l22 = (0.01 * pk.AU)**2 130 flyby_constraints.append(1. - d2 / l22) 131 cineq.extend(flyby_constraints) 132 133 mismatch_constraints.extend([a - b for a, b in zip(rfwd[-1], rbwd[0])]) 134 mismatch_constraints.extend([a - b for a, b in zip(vfwd[-1], vbwd[0])]) 135 mismatch_constraints.append(mfwd[-1] - mbwd[0]) 136 ceq.extend(mismatch_constraints) 137 138 # Making the mismatches non dimensional 139 ceq[0] /= pk.AU 140 ceq[1] /= pk.AU 141 ceq[2] /= pk.AU 142 ceq[3] /= pk.EARTH_VELOCITY 143 ceq[4] /= pk.EARTH_VELOCITY 144 ceq[5] /= pk.EARTH_VELOCITY 145 ceq[6] /= self.__sc.mass 146 147 # We assemble the constraint vector 148 retval = [obj] 149 retval.extend(ceq) 150 retval.extend(cineq) 151 152 return retval 153 154 # Get bounds 155 def get_bounds(self): 156 return (self.__lb, self.__ub) 157 158 # Get number of inequality contraints 159 def get_nic(self): 160 if self.__earth_gravity: 161 return 2 * self.__n_seg 162 else: 163 return self.__n_seg 164 165 # Get number of equality contraints 166 def get_nec(self): 167 return 7 168 169 # Gradient sparsity 170 def gradient_sparsity(self): 171 # Objective 172 retval = [[0, 2]] 173 # Mismatches 174 retval += [[i, j] for i in range(1, 7) 175 for j in range(3 * (1 + self.__n_seg))] 176 retval += [[7, j] for j in range(1, 3 * (1 + self.__n_seg))] 177 # Throttles 178 retval += [[8 + i, 3 + 3 * i + j] 179 for i in range(self.__n_seg) for j in range(3)] 180 # Prevent flyby 181 if self.__earth_gravity: 182 retval += [[8 + self.__n_seg + i, j] 183 for i in range(self.__fwd_seg) for j in [0, 1] + list(range(3, 6 + 3 * i))] 184 retval += [[8 + self.__n_seg + self.__fwd_seg + i, j] for i in range(self.__bwd_seg) for j in [ 185 0, 1, 2] + list(range(3 + 3 * (self.__n_seg - self.__bwd_seg + i), 3 + 3 * self.__n_seg))] 186 return retval 187 188 def get_name(self): 189 return "MARGO cubesat transfer to " + self.target.name 190 191 def get_extra_info(self): 192 retval = "\tTarget planet: " + self.target.name 193 retval += "\n\tEarth gravity: " + str(self.__earth_gravity) 194 retval += "\n\tSolar Electric Propulsion: " + str(self.__sep) 195 retval += "\n\tStart mass: " + str(self.__sc.mass) + " kg" 196 retval += "\n\tMaximum thrust as 1AU: " + str(self.__sc.thrust) + " N" 197 retval += "\n\tSpecific impulse: " + str(self.__sc.isp) + " s" 198 retval += "\n\n\tLaunch window: [" + \ 199 str(self.__lb[0]) + ", " + str(self.__ub[0]) + "] - MJD2000" 200 retval += "\n\tBounds on time of flight: [" + str( 201 self.__lb[1]) + ", " + str(self.__ub[1]) + "] - days" 202 retval += "\n\n\tNumber of segments: " + str(self.__n_seg) 203 retval += "\n\tGrid type: " + self.__grid_type 204 205 return retval 206 207 # SEP model 208 def _sep_model(self, r): 209 SAA = 0 210 Psupply = 13.75 211 eff = 0.92 212 213 Pbmp = (-40.558 * r**3 + 173.49 * r**2 - 214 259.19 * r + 141.86) * math.cos(SAA) 215 P = -146.26 * r**3 + 658.52 * r**2 - 1059.2 * r + 648.24 # 6 panels 216 # P = -195.02 * r**3 + 878.03 * r**2 - 1412.3 * r + 864.32 # 8 panels 217 if Pbmp < Psupply: 218 P -= (Psupply - Pbmp) 219 Pin = eff * P 220 if Pin > 120: 221 Pin = 120 # thermal max 120W 222 223 Tmax = (26.27127 * Pin - 708.973) / 1000000 224 if Tmax < 0: 225 Tmax = 0 226 227 Isp = -0.0011 * Pin**3 + 0.175971 * Pin**2 + 4.193797 * Pin + 2037.213 228 229 return Tmax, Isp 230 231 # Propagates the trajectory 232 def _propagate(self, x): 233 # 1 - We decode the chromosome 234 t0 = x[0] 235 T = x[1] 236 m_f = x[2] 237 # We extract the number of segments for forward and backward 238 # propagation 239 n_seg = self.__n_seg 240 fwd_seg = self.__fwd_seg 241 bwd_seg = self.__bwd_seg 242 # We extract information on the spacecraft 243 m_i = self.__sc.mass 244 max_thrust = self.__sc.thrust 245 isp = self.__sc.isp 246 veff = isp * pk.G0 247 # And on the leg 248 throttles = [x[3 + 3 * i: 6 + 3 * i] for i in range(n_seg)] 249 # Return lists 250 n_points_fwd = fwd_seg + 1 251 n_points_bwd = bwd_seg + 1 252 rfwd = [[0.0] * 3] * (n_points_fwd) 253 vfwd = [[0.0] * 3] * (n_points_fwd) 254 mfwd = [0.0] * (n_points_fwd) 255 ufwd = [[0.0] * 3] * (fwd_seg) 256 dfwd = [[0.0] * 3] * (fwd_seg) # distances E/S 257 rbwd = [[0.0] * 3] * (n_points_bwd) 258 vbwd = [[0.0] * 3] * (n_points_bwd) 259 mbwd = [0.0] * (n_points_bwd) 260 ubwd = [[0.0] * 3] * (bwd_seg) 261 dbwd = [[0.0] * 3] * (bwd_seg) 262 263 # 2 - We compute the initial and final epochs and ephemerides 264 t_i = pk.epoch(t0) 265 r_i, v_i = self.__earth.eph(t_i) 266 if self.__start == 'l1': 267 r_i = [r * 0.99 for r in r_i] 268 v_i = [v * 0.99 for v in v_i] 269 elif self.__start == 'l2': 270 r_i = [r * 1.01 for r in r_i] 271 v_i = [v * 1.01 for v in v_i] 272 t_f = pk.epoch(t0 + T) 273 r_f, v_f = self.target.eph(t_f) 274 275 # 3 - Forward propagation 276 fwd_grid = t0 + T * self.__fwd_grid # days 277 fwd_dt = T * self.__fwd_dt # seconds 278 # Initial conditions 279 rfwd[0] = r_i 280 vfwd[0] = v_i 281 mfwd[0] = m_i 282 # Propagate 283 for i, t in enumerate(throttles[:fwd_seg]): 284 if self.__sep: 285 r = math.sqrt(rfwd[i][0]**2 + rfwd[i][1] 286 ** 2 + rfwd[i][2]**2) / pk.AU 287 max_thrust, isp = self._sep_model(r) 288 veff = isp * pk.G0 289 ufwd[i] = [max_thrust * thr for thr in t] 290 if self.__earth_gravity: 291 r_E, v_E = self.__earth.eph(pk.epoch(fwd_grid[i])) 292 dfwd[i] = [a - b for a, b in zip(r_E, rfwd[i])] 293 r3 = sum([r**2 for r in dfwd[i]])**(3 / 2) 294 disturbance = [mfwd[i] * pk.MU_EARTH / 295 r3 * ri for ri in dfwd[i]] 296 rfwd[i + 1], vfwd[i + 1], mfwd[i + 1] = pk.propagate_taylor_disturbance( 297 rfwd[i], vfwd[i], mfwd[i], ufwd[i], disturbance, fwd_dt[i], pk.MU_SUN, veff, -10, -10) 298 else: 299 rfwd[i + 1], vfwd[i + 1], mfwd[i + 1] = pk.propagate_taylor( 300 rfwd[i], vfwd[i], mfwd[i], ufwd[i], fwd_dt[i], pk.MU_SUN, veff, -10, -10) 301 302 # 4 - Backward propagation 303 bwd_grid = t0 + T * self.__bwd_grid # days 304 bwd_dt = T * self.__bwd_dt # seconds 305 # Final conditions 306 rbwd[-1] = r_f 307 vbwd[-1] = v_f 308 mbwd[-1] = m_f 309 # Propagate 310 for i, t in enumerate(throttles[-1:-bwd_seg - 1:-1]): 311 if self.__sep: 312 r = math.sqrt(rbwd[-i - 1][0]**2 + rbwd[-i - 1] 313 [1]**2 + rbwd[-i - 1][2]**2) / pk.AU 314 max_thrust, isp = self._sep_model(r) 315 veff = isp * pk.G0 316 ubwd[-i - 1] = [max_thrust * thr for thr in t] 317 if self.__earth_gravity: 318 r_E, v_E = self.__earth.eph(pk.epoch(bwd_grid[-i - 1])) 319 dbwd[-i - 1] = [a - b for a, b in zip(r_E, rbwd[-i - 1])] 320 r3 = sum([r**2 for r in dbwd[-i - 1]])**(3 / 2) 321 disturbance = [mfwd[i] * pk.MU_EARTH / 322 r3 * ri for ri in dbwd[-i - 1]] 323 rbwd[-i - 2], vbwd[-i - 2], mbwd[-i - 2] = pk.propagate_taylor_disturbance( 324 rbwd[-i - 1], vbwd[-i - 1], mbwd[-i - 1], ubwd[-i - 1], disturbance, -bwd_dt[-i - 1], pk.MU_SUN, veff, -10, -10) 325 else: 326 rbwd[-i - 2], vbwd[-i - 2], mbwd[-i - 2] = pk.propagate_taylor( 327 rbwd[-i - 1], vbwd[-i - 1], mbwd[-i - 1], ubwd[-i - 1], -bwd_dt[-i - 1], pk.MU_SUN, veff, -10, -10) 328 329 return rfwd, rbwd, vfwd, vbwd, mfwd, mbwd, ufwd, ubwd, fwd_dt, bwd_dt, dfwd, dbwd 330 331 # Visualizes the trajectory 332 def plot_traj(self, x, units=pk.AU, plot_segments=True, plot_thrusts=False, axes=None): 333 """ 334 ax = prob.plot_traj(self, x, units=AU, plot_segments=True, plot_thrusts=False, axes=None) 335 336 Args: 337 - x (``list``, ``tuple``, ``numpy.ndarray``): Decision chromosome, e.g. (``pygmo.population.champion_x``). 338 - units (``float``): the length unit to be used in the plot 339 - plot_segments (``bool``): when True plots also the segments boundaries 340 - plot_thrusts (``bool``): when True plots also the thrust vectors 341 342 Returns: 343 matplotlib.axes: axes where to plot 344 345 Visualizes the the trajectory in a 3D plot 346 """ 347 348 if not len(x) == len(self.get_bounds()[0]): 349 raise ValueError("Invalid length of the decision vector x") 350 351 import matplotlib as mpl 352 import matplotlib.pyplot as plt 353 from mpl_toolkits.mplot3d import Axes3D 354 355 # Creating the axes if necessary 356 if axes is None: 357 mpl.rcParams['legend.fontsize'] = 10 358 fig = plt.figure() 359 axes = fig.gca(projection='3d') 360 361 n_seg = self.__n_seg 362 fwd_seg = self.__fwd_seg 363 bwd_seg = self.__bwd_seg 364 t0 = x[0] 365 T = x[1] 366 isp = self.__sc.isp 367 veff = isp * pk.G0 368 fwd_grid = t0 + T * self.__fwd_grid # days 369 bwd_grid = t0 + T * self.__bwd_grid # days 370 371 throttles = [x[3 + 3 * i: 6 + 3 * i] for i in range(n_seg)] 372 alphas = [min(1., np.linalg.norm(t)) for t in throttles] 373 374 times = np.concatenate((fwd_grid, bwd_grid)) 375 376 rfwd, rbwd, vfwd, vbwd, mfwd, mbwd, ufwd, ubwd, fwd_dt, bwd_dt, dfwd, dbwd = self._propagate( 377 x) 378 379 # Plotting the Sun, the Earth and the target 380 axes.scatter([0], [0], [0], color='y') 381 pk.orbit_plots.plot_planet(self.__earth, pk.epoch( 382 t0), units=units, legend=True, color=(0.7, 0.7, 1), axes=axes) 383 pk.orbit_plots.plot_planet(self.target, pk.epoch( 384 t0 + T), units=units, legend=True, color=(0.7, 0.7, 1), axes=axes) 385 386 # Forward propagation 387 xfwd = [0.0] * (fwd_seg + 1) 388 yfwd = [0.0] * (fwd_seg + 1) 389 zfwd = [0.0] * (fwd_seg + 1) 390 xfwd[0] = rfwd[0][0] / units 391 yfwd[0] = rfwd[0][1] / units 392 zfwd[0] = rfwd[0][2] / units 393 394 for i in range(fwd_seg): 395 if self.__sep: 396 r = math.sqrt(rfwd[i][0]**2 + rfwd[i][1] 397 ** 2 + rfwd[i][2]**2) / pk.AU 398 _, isp = self._sep_model(r) 399 veff = isp * pk.G0 400 if self.__earth_gravity: 401 r3 = sum([r**2 for r in dfwd[i]])**(3 / 2) 402 disturbance = [mfwd[i] * pk.MU_EARTH / 403 r3 * ri for ri in dfwd[i]] 404 pk.orbit_plots.plot_taylor_disturbance(rfwd[i], vfwd[i], mfwd[i], ufwd[i], disturbance, fwd_dt[ 405 i], pk.MU_SUN, veff, N=10, units=units, color=(alphas[i], 0, 1 - alphas[i]), axes=axes) 406 else: 407 pk.orbit_plots.plot_taylor(rfwd[i], vfwd[i], mfwd[i], ufwd[i], fwd_dt[ 408 i], pk.MU_SUN, veff, N=10, units=units, color=(alphas[i], 0, 1 - alphas[i]), axes=axes) 409 xfwd[i + 1] = rfwd[i + 1][0] / units 410 yfwd[i + 1] = rfwd[i + 1][1] / units 411 zfwd[i + 1] = rfwd[i + 1][2] / units 412 if plot_segments: 413 axes.scatter(xfwd[:-1], yfwd[:-1], zfwd[:-1], 414 label='nodes', marker='o', s=5, c='k') 415 416 # Backward propagation 417 xbwd = [0.0] * (bwd_seg + 1) 418 ybwd = [0.0] * (bwd_seg + 1) 419 zbwd = [0.0] * (bwd_seg + 1) 420 xbwd[-1] = rbwd[-1][0] / units 421 ybwd[-1] = rbwd[-1][1] / units 422 zbwd[-1] = rbwd[-1][2] / units 423 424 for i in range(bwd_seg): 425 if self.__sep: 426 r = math.sqrt(rbwd[-i - 1][0]**2 + rbwd[-i - 1] 427 [1]**2 + rbwd[-i - 1][2]**2) / pk.AU 428 _, isp = self._sep_model(r) 429 veff = isp * pk.G0 430 if self.__earth_gravity: 431 r3 = sum([r**2 for r in dbwd[-i - 1]])**(3 / 2) 432 disturbance = [mfwd[i] * pk.MU_EARTH / 433 r3 * ri for ri in dbwd[-i - 1]] 434 pk.orbit_plots.plot_taylor_disturbance(rbwd[-i - 1], vbwd[-i - 1], mbwd[-i - 1], ubwd[-i - 1], disturbance, -bwd_dt[ 435 -i - 1], pk.MU_SUN, veff, N=10, units=units, color=(alphas[-i - 1], 0, 1 - alphas[-i - 1]), axes=axes) 436 else: 437 pk.orbit_plots.plot_taylor(rbwd[-i - 1], vbwd[-i - 1], mbwd[-i - 1], ubwd[-i - 1], -bwd_dt[-i - 1], 438 pk.MU_SUN, veff, N=10, units=units, color=(alphas[-i - 1], 0, 1 - alphas[-i - 1]), axes=axes) 439 xbwd[-i - 2] = rbwd[-i - 2][0] / units 440 ybwd[-i - 2] = rbwd[-i - 2][1] / units 441 zbwd[-i - 2] = rbwd[-i - 2][2] / units 442 if plot_segments: 443 axes.scatter(xbwd[1:], ybwd[1:], zbwd[1:], marker='o', s=5, c='k') 444 445 # Plotting the thrust vectors 446 if plot_thrusts: 447 throttles = np.array(throttles) 448 xlim = axes.get_xlim() 449 xrange = xlim[1] - xlim[0] 450 ylim = axes.get_ylim() 451 yrange = ylim[1] - ylim[0] 452 zlim = axes.get_zlim() 453 zrange = zlim[1] - zlim[0] 454 455 scale = 0.1 456 457 throttles[:, 0] *= xrange 458 throttles[:, 1] *= yrange 459 throttles[:, 2] *= zrange 460 461 throttles *= scale 462 463 for (x, y, z, t) in zip(xfwd[:-1] + xbwd[:-1], yfwd[:-1] + ybwd[:-1], zfwd[:-1] + zbwd[:-1], throttles): 464 axes.plot([x, x + t[0]], [y, y + t[1]], [z, z + t[2]], c='g') 465 466 return axes 467 468 def plot_dists_thrust(self, x, axes=None): 469 """ 470 axes = prob.plot_dists_thrust(self, x, axes=None) 471 472 Args: 473 - x (``list``, ``tuple``, ``numpy.ndarray``): Decision chromosome, e.g. (``pygmo.population.champion_x``). 474 475 Returns: 476 matplotlib.axes: axes where to plot 477 478 Plots the distance of the spacecraft from the Earth/Sun and the thrust profile 479 """ 480 481 if not len(x) == len(self.get_bounds()[0]): 482 raise ValueError("Invalid length of the decision vector x") 483 484 import matplotlib as mpl 485 import matplotlib.pyplot as plt 486 487 # Creating the axes if necessary 488 if axes is None: 489 mpl.rcParams['legend.fontsize'] = 10 490 fig = plt.figure() 491 axes = fig.add_subplot(111) 492 493 n_seg = self.__n_seg 494 fwd_seg = self.__fwd_seg 495 bwd_seg = self.__bwd_seg 496 t0 = x[0] 497 T = x[1] 498 fwd_grid = t0 + T * self.__fwd_grid # days 499 bwd_grid = t0 + T * self.__bwd_grid # days 500 501 throttles = [np.linalg.norm(x[3 + 3 * i: 6 + 3 * i]) 502 for i in range(n_seg)] 503 504 dist_earth = [0.0] * (n_seg + 2) # distances spacecraft - Earth 505 dist_sun = [0.0] * (n_seg + 2) # distances spacecraft - Sun 506 times = np.concatenate((fwd_grid, bwd_grid)) 507 508 rfwd, rbwd, vfwd, vbwd, mfwd, mbwd, ufwd, ubwd, fwd_dt, bwd_dt, _, _ = self._propagate( 509 x) 510 511 # Forward propagation 512 xfwd = [0.0] * (fwd_seg + 1) 513 yfwd = [0.0] * (fwd_seg + 1) 514 zfwd = [0.0] * (fwd_seg + 1) 515 xfwd[0] = rfwd[0][0] / pk.AU 516 yfwd[0] = rfwd[0][1] / pk.AU 517 zfwd[0] = rfwd[0][2] / pk.AU 518 r_E = [ri / pk.AU for ri in self.__earth.eph(pk.epoch(fwd_grid[0]))[0]] 519 dist_earth[0] = np.linalg.norm( 520 [r_E[0] - xfwd[0], r_E[1] - yfwd[0], r_E[2] - zfwd[0]]) 521 dist_sun[0] = np.linalg.norm([xfwd[0], yfwd[0], zfwd[0]]) 522 523 for i in range(fwd_seg): 524 xfwd[i + 1] = rfwd[i + 1][0] / pk.AU 525 yfwd[i + 1] = rfwd[i + 1][1] / pk.AU 526 zfwd[i + 1] = rfwd[i + 1][2] / pk.AU 527 r_E = [ 528 ri / pk.AU for ri in self.__earth.eph(pk.epoch(fwd_grid[i + 1]))[0]] 529 dist_earth[ 530 i + 1] = np.linalg.norm([r_E[0] - xfwd[i + 1], r_E[1] - yfwd[i + 1], r_E[2] - zfwd[i + 1]]) 531 dist_sun[ 532 i + 1] = np.linalg.norm([xfwd[i + 1], yfwd[i + 1], zfwd[i + 1]]) 533 534 # Backward propagation 535 xbwd = [0.0] * (bwd_seg + 1) 536 ybwd = [0.0] * (bwd_seg + 1) 537 zbwd = [0.0] * (bwd_seg + 1) 538 xbwd[-1] = rbwd[-1][0] / pk.AU 539 ybwd[-1] = rbwd[-1][1] / pk.AU 540 zbwd[-1] = rbwd[-1][2] / pk.AU 541 r_E = [ 542 ri / pk.AU for ri in self.__earth.eph(pk.epoch(bwd_grid[-1]))[0]] 543 dist_earth[-1] = np.linalg.norm([r_E[0] - xbwd[-1], 544 r_E[1] - ybwd[-1], r_E[2] - zbwd[-1]]) 545 dist_sun[-1] = np.linalg.norm([xbwd[-1], ybwd[-1], zbwd[-1]]) 546 547 for i in range(bwd_seg): 548 xbwd[-i - 2] = rbwd[-i - 2][0] / pk.AU 549 ybwd[-i - 2] = rbwd[-i - 2][1] / pk.AU 550 zbwd[-i - 2] = rbwd[-i - 2][2] / pk.AU 551 r_E = [ 552 ri / pk.AU for ri in self.__earth.eph(pk.epoch(bwd_grid[-i - 2]))[0]] 553 dist_earth[-i - 2] = np.linalg.norm( 554 [r_E[0] - xbwd[-i - 2], r_E[1] - ybwd[-i - 2], r_E[2] - zbwd[-i - 2]]) 555 dist_sun[-i - 556 2] = np.linalg.norm([xbwd[-i - 2], ybwd[-i - 2], zbwd[-i - 2]]) 557 558 axes.set_xlabel("t [mjd2000]") 559 # Plot Earth distance 560 axes.plot(times, dist_earth, c='b', label="sc-Earth") 561 # Plot Sun distance 562 axes.plot(times, dist_sun, c='y', label="sc-Sun") 563 axes.set_ylabel("distance [AU]", color='k') 564 axes.set_ylim(bottom=0.) 565 axes.tick_params('y', colors='k') 566 axes.legend(loc=2) 567 # draw threshold where Earth gravity equals 0.1*Tmax 568 if self.__earth_gravity: 569 axes.axhline(y=np.sqrt(pk.MU_EARTH * self.__sc.mass / (self.__sc.thrust * 0.1) 570 ) / pk.AU, c='g', ls=":", lw=1, label="earth_g = 0.1*Tmax") 571 # Plot thrust profile 572 axes = axes.twinx() 573 if self.__sep: 574 max_thrust = self.__sc.thrust 575 thrusts = np.linalg.norm( 576 np.array(ufwd + ubwd), axis=1) / max_thrust 577 # plot maximum thrust achievable at that distance from the Sun 578 distsSun = dist_sun[:fwd_seg] + \ 579 dist_sun[-bwd_seg:] + [dist_sun[-1]] 580 Tmaxs = [self._sep_model(d)[0] / max_thrust for d in distsSun] 581 axes.step(np.concatenate( 582 (fwd_grid, bwd_grid[1:])), Tmaxs, where="post", c='lightgray', linestyle=':') 583 else: 584 thrusts = throttles.copy() 585 # duplicate the last for plotting 586 thrusts = np.append(thrusts, thrusts[-1]) 587 axes.step(np.concatenate( 588 (fwd_grid, bwd_grid[1:])), thrusts, where="post", c='r', linestyle='--') 589 axes.set_ylabel("T/Tmax$_{1AU}$", color='r') 590 axes.tick_params('y', colors='r') 591 axes.set_xlim([times[0], times[-1]]) 592 axes.set_ylim([0, max(thrusts) + 0.2]) 593 594 return axes 595 596 def double_segments(self, x): 597 """ 598 new_prob, new_x = prob.double_segments(self,x) 599 600 Args: 601 - x (``list``, ``tuple``, ``numpy.ndarray``): Decision chromosome, e.g. (``pygmo.population.champion_x``). 602 603 Returns: 604 - the new udp having twice the segments 605 - list: the new chromosome to be used as initial guess 606 607 Returns the decision vector encoding a low trust trajectory having double the number of segments with respect to x 608 and a 'similar' throttle history. In case high fidelity is True, and x is a feasible trajectory, the returned decision vector 609 also encodes a feasible trajectory that can be further optimized 610 Returns the problem and the decision vector encoding a low-thrust trajectory having double the number of 611 segments with respect to x and the same throttle history. If x is a feasible trajectory, the new chromosome is "slightly 612 unfeasible", due to the new refined Earth's gravity and/or SEP thrust that are now computed in the 2 halves of each segment. 613 """ 614 615 if not len(x) == len(self.get_bounds()[0]): 616 raise ValueError("Invalid length of the decision vector x") 617 618 new_x = np.append(x[:3], np.repeat(x[3:].reshape((-1, 3)), 2, axis=0)) 619 620 new_prob = lt_margo( 621 target=self.target, 622 n_seg=2 * self.__n_seg, 623 grid_type=self.__grid_type, 624 t0=[pk.epoch(self.__lb[0]), pk.epoch(self.__ub[0])], 625 tof=[self.__lb[1], self.__ub[1]], 626 m0=self.__sc.mass, 627 Tmax=self.__sc.thrust, 628 Isp=self.__sc.isp, 629 earth_gravity=self.__earth_gravity, 630 sep=self.__sep, 631 start=self.__start 632 ) 633 634 return new_prob, new_x 635 636 def pretty(self, x): 637 """ 638 prob.pretty(x) 639 640 Args: 641 - x (``list``, ``tuple``, ``numpy.ndarray``): Decision chromosome, e.g. (``pygmo.population.champion_x``). 642 643 Prints human readable information on the trajectory represented by the decision vector x 644 """ 645 646 if not len(x) == len(self.get_bounds()[0]): 647 raise ValueError("Invalid length of the decision vector x") 648 649 n_seg = self.__n_seg 650 m_i = self.__sc.mass 651 t0 = x[0] 652 T = x[1] 653 m_f = x[2] 654 thrusts = [np.linalg.norm(x[3 + 3 * i: 6 + 3 * i]) 655 for i in range(n_seg)] 656 657 tf = t0 + T 658 mP = m_i - m_f 659 deltaV = self.__sc.isp * pk.G0 * np.log(m_i / m_f) 660 661 dt = np.append(self.__fwd_dt, self.__bwd_dt) * T / pk.DAY2SEC 662 time_thrusts_on = sum(dt[i] for i in range( 663 len(thrusts)) if thrusts[i] > 0.1) 664 665 print("Departure:", pk.epoch(t0), "(", t0, "mjd2000 )") 666 print("Time of flight:", T, "days") 667 print("Arrival:", pk.epoch(tf), "(", tf, "mjd2000 )") 668 print("Delta-v:", deltaV, "m/s") 669 print("Propellant consumption:", mP, "kg") 670 print("Thrust-on time:", time_thrusts_on, "days") 671