MEE-SMA Q-law: GTO to GEO

One of the most common low-thrust many revolution transfer is to go from GTO to GEO. In this example, we make use of pyqlaw to construct this transfer.

import numpy as np
from numpy.random import rand
import matplotlib.pyplot as plt
import time

import sys
sys.path.append("../../../")    # path to pyqlaw
import pyqlaw

We begin by defining our canonical scales and initial conditions.

# initial and final elements: [a,e,i,RAAN,omega,ta]
LU = 42164.0
GM_EARTH = 398600.44
VU = np.sqrt(GM_EARTH/LU)
TU = LU/VU

rp_gto = 200 + 6378
ra_gto = 35786 + 6378
sma_gto = (rp_gto + ra_gto)/(2*LU)
ecc_gto = (ra_gto - rp_gto)/(ra_gto + rp_gto)
KEP0 = [sma_gto,ecc_gto,np.deg2rad(23),0,0,0]
KEPF = [1,0,np.deg2rad(3),0,0,0]

# spacecraft parameters
MU = 1000        # spacecraft wet mass, kg
tmax_si = 0.45   # spacecraft thrust, Newton
isp_si  = 1500   # spacecraft specific impulse, seconds
mdot_si = tmax_si/(isp_si*9.81)  # kg/s

# non-dimensional quantities
mass0 = 1.0
tmax = tmax_si * (1/MU)*(TU**2/(1e3*LU))
mdot = np.abs(mdot_si) *(TU/MU)
tmax, mdot

(0.0020070507277914697, 0.0004193687522487825)

Because GEO is circular and planar, Keplerian elements are not good representations for the Lyapunov controller due to the associated singularities. We will instead use the state representation based on modified equinoctial elements, but with the semilatus rectum replaced by semi-major axis:

\[\boldsymbol{x} = [a,f,g,h,k,L]^T\]

This set of orbital elements is called mee_with_a within pyqlaw. We will convert our initial Keplerian elements to this representation.

oe0 = pyqlaw.kep2mee_with_a(np.array(KEP0))
oeT = pyqlaw.kep2mee_with_a(np.array(KEPF))

We can now construct the Q-law problem object, making sure that we choose the appropriate elements_type. We also have a choice of integrator, where we can either choose:

• "rk4": fixed-step 4th order Runge-Kutta, or

• "rkf45": 4th order Runge-Kutta with 5th order step-size correction

In addition, by setting use_sundman = True, we are taking steps in terms of fixed angle (i.e. fixed anomaly) rather than time. This is particularly useful to decrease computational time (by reducing number of integration steps) without compromising the accuracy too much for highly elliptical orbits, where small time-steps are required near periapsis, while large time-steps can be used near apoapsis.

For now, we will assume Keplerian dynamics, so we do not pass any perturbations (we set it to None - which is the default anyway).

prob = pyqlaw.QLaw(
    integrator="rk4",
    elements_type="mee_with_a",
    verbosity=2,
    print_frequency=2000,
    use_sundman=True,
    perturbations=None,
)

We now “set” the problem, i.e. we assign it its initial and final conditions, initial mass, spacecraft engine parameters, etc. One important parameter to keep in mind is woe - this is the weights associated to each state (i.e. orbital elements) component within the Lyapunov function. A larger weight means that the controller will put higher emphasis on correcting this particular component.

# setup problem
tf_max = 365.25 * 86400 / TU       # max time, in canonical scales
t_step = np.deg2rad(5)             # integration step, in angles
woe = [3.0, 1.0, 1.0, 1.0, 1.0]
prob.set_problem(oe0, oeT, mass0, tmax, mdot, tf_max, t_step,
    woe = woe)
prob.pretty()

Transfer:
  a  : 5.7800e-01 -> 1.0000e+00 (weight: 3.00)
  f  : 7.3009e-01 -> 0.0000e+00 (weight: 1.00)
  g  : 0.0000e+00 -> 0.0000e+00 (weight: 1.00)
  h  : 2.0345e-01 -> 2.6186e-02 (weight: 1.00)
  k  : 0.0000e+00 -> 0.0000e+00 (weight: 1.00)

Time-optimal transfer

We can now solve for the transfer. We will begin by constructing a solution that will always be thrusting. In optimal control problem (OCP) jargon, this is a bang-bang (as opposed to bang-off-bang) solution, and is the control history found for time-optimal problems. Of course, we must keep in mind that we are not solving an OCP here, so the solution we will obtain is sub-optimal.

# solve
prob.solve()
prob.pretty_results()

 iter   |  time      |  del1       |  del2       |  del3       |  del4       |  del5       |  el6        |
      0 |  1.039e-02 | -4.2198e-01 |  7.3010e-01 |  0.0000e+00 |  1.7726e-01 |  0.0000e+00 |  6.2435e-02 |
   2000 |  7.486e+01 | -3.3556e-01 |  6.7217e-01 |  5.8932e-02 |  1.2068e-01 |  1.1470e-02 |  1.6575e+02 |
   4000 |  1.730e+02 | -1.9396e-01 |  5.7053e-01 |  1.1938e-01 |  4.8560e-02 |  1.6363e-02 |  3.2133e+02 |
   6000 |  3.110e+02 | -4.5635e-02 |  2.5424e-01 |  1.0588e-01 |  4.8795e-03 |  6.4765e-03 |  4.7994e+02 |
Target elements successfully reached!
Exit code : 2
Converge  : True
Final state:
  a  : 1.0018e+00 (error: 1.8276e-03)
  f  : -2.4312e-03 (error: 2.4312e-03)
  g  : 5.1535e-03 (error: 5.1535e-03)
  h  : 2.6192e-02 (error: 5.9591e-06)
  k  : 3.1508e-06 (error: 3.1508e-06)
Transfer time : 394.14835992080316
Final mass    : 0.8347064940991077

We can now plot and visualize the results:

# plots
fig1, ax1 = prob.plot_elements_history(to_keplerian=True, TU=TU/86400, time_unit_name="day")
fig2, ax2 = prob.plot_trajectory_3d(interpolate=False, sphere_radius=6378/LU)
fig3, ax3 = prob.plot_controls()
plt.show()

Introducing throttles

We can see in the third plot from the previous cell that the throttle is always at 100% - as expected. We can now introduce efficiency thresholds when we call the solve() method, (potentially) resulting in mass savings. It’s important to remember that introducing efficiency thresholds does not necessarily mean mass savings in all cases - because we thrust less often, it will always take longer to reach the final targeted elements - and the fact that it takes longer also means we will be thrusting for longer duration in total, thus resulting in more mass expenditure.

# solve
prob.solve(eta_r=0.1)
prob.pretty_results()

 iter   |  time      |  del1       |  del2       |  del3       |  del4       |  del5       |  el6        |
      0 |  1.039e-02 | -4.2198e-01 |  7.3010e-01 |  0.0000e+00 |  1.7726e-01 |  0.0000e+00 |  6.2435e-02 |
   2000 |  7.132e+01 | -3.7904e-01 |  7.2973e-01 |  3.1390e-02 |  1.6966e-01 |  5.4212e-03 |  1.7083e+02 |
   4000 |  1.509e+02 | -3.3565e-01 |  7.2378e-01 |  6.7669e-02 |  1.6053e-01 |  1.1689e-02 |  3.3571e+02 |
   6000 |  2.386e+02 | -2.9457e-01 |  7.1253e-01 |  1.0511e-01 |  1.5043e-01 |  1.8040e-02 |  4.9833e+02 |
   8000 |  3.394e+02 | -2.4926e-01 |  6.9347e-01 |  1.4888e-01 |  1.3837e-01 |  2.5773e-02 |  6.6352e+02 |
  10000 |  4.457e+02 | -2.0973e-01 |  6.6511e-01 |  1.9186e-01 |  1.2434e-01 |  3.2888e-02 |  8.2371e+02 |
  12000 |  5.572e+02 | -1.7922e-01 |  6.2490e-01 |  2.3206e-01 |  1.0796e-01 |  3.8371e-02 |  9.7282e+02 |
  14000 |  6.790e+02 | -1.4735e-01 |  5.6812e-01 |  2.7032e-01 |  8.8035e-02 |  4.3422e-02 |  1.1310e+03 |
  16000 |  8.097e+02 | -1.2048e-01 |  4.8337e-01 |  2.9958e-01 |  6.4930e-02 |  4.5387e-02 |  1.2906e+03 |
  18000 |  9.468e+02 | -9.9431e-02 |  3.6218e-01 |  2.9628e-01 |  4.0623e-02 |  3.8171e-02 |  1.4456e+03 |
  20000 |  1.095e+03 | -5.9155e-02 |  1.9291e-01 |  2.3226e-01 |  1.9019e-02 |  2.0033e-02 |  1.6066e+03 |
Target elements successfully reached!
Exit code : 2
Converge  : True
Final state:
  a  : 9.9766e-01 (error: 2.3445e-03)
  f  : -6.9620e-03 (error: 6.9620e-03)
  g  : 3.9621e-03 (error: 3.9621e-03)
  h  : 2.6210e-02 (error: 2.4209e-05)
  k  : 4.3703e-05 (error: 4.3703e-05)
Transfer time : 1249.148452090831
Final mass    : 0.7428518861140109

# plots
fig1, ax1 = prob.plot_elements_history(to_keplerian=True, TU=TU/86400, time_unit_name="day")
fig2, ax2 = prob.plot_trajectory_3d(interpolate=False, sphere_radius=6378/LU)
fig3, ax3 = prob.plot_controls()
fig4, ax4 = prob.plot_efficiency(TU=TU/86400, time_unit_name="day")
plt.show()

Introducing Non-Keplerian perturbations

We will now consider non-Keplerian perturbations into our problem. This requires us to go back to the initialization of the pyqlaw.QLaw object.

To compute perturbations, we require some SPICE kernels (to compute the position of the planets/Sun, as well as to compute the orientation of the principal axes frame with respect to the inertial frame.

All kernels furnished (i.e. loaded) in the next cell can be downloaded on NAIF’s website: https://naif.jpl.nasa.gov/pub/naif/generic_kernels/

import os
import spiceypy as spice
spice.furnsh(os.path.join(os.getenv("SPICE"), "lsk", "naif0012.tls"))
spice.furnsh(os.path.join(os.getenv("SPICE"), "spk", "de440.bsp"))
spice.furnsh(os.path.join(os.getenv("SPICE"), "pck", "gm_de440.tpc"))
spice.furnsh(os.path.join(os.getenv("SPICE"), "pck", "earth_200101_990825_predict.bpc"))
spice.furnsh(os.path.join(os.getenv("SPICE"), "fk", "earth_assoc_itrf93.tf"))

We can now construct the perturbations struct. As input, we need to also provide the initial epoch, i.e. the epoch which corresponds to time t = 0 within the Q-law algorithm.

# initialize object for perturbations
et_ref = spice.str2et("2028-01-01T00:00:00")   # seconds past J2000
perturbations = pyqlaw.SpicePerturbations(
    et_ref, LU, TU,
    frame_qlaw='J2000',
    J2=0.00108263,
)

We can now construct the problem object, making sure we pass perturbations

# construct problem with perturbations
prob = pyqlaw.QLaw(
    integrator="rk4",
    elements_type="mee_with_a",
    verbosity=2,
    print_frequency=3000,
    use_sundman = True,
    perturbations = perturbations,
)

The remaining is the same as before - we set_problem(), then solve()!

prob.set_problem(oe0, oeT, mass0, tmax, mdot, tf_max, t_step,
    woe = woe)
prob.solve()
prob.pretty_results()

 iter   |  time      |  del1       |  del2       |  del3       |  del4       |  del5       |  el6        |
      0 |  1.039e-02 | -4.2198e-01 |  7.3010e-01 | -5.7078e-09 |  1.7726e-01 |  0.0000e+00 |  6.2435e-02 |
   3000 |  1.205e+02 | -2.7214e-01 |  6.3139e-01 |  9.1438e-02 |  8.5709e-02 |  1.5508e-02 |  2.4237e+02 |
   6000 |  3.112e+02 | -4.4788e-02 |  2.5548e-01 |  1.0858e-01 |  4.9118e-03 |  6.7834e-03 |  4.8033e+02 |
Target elements successfully reached!
Exit code : 2
Converge  : True
Final state:
  a  : 1.0021e+00 (error: 2.1131e-03)
  f  : -3.4311e-03 (error: 3.4311e-03)
  g  : 6.5437e-03 (error: 6.5437e-03)
  h  : 2.6208e-02 (error: 2.1648e-05)
  k  : 2.4339e-05 (error: 2.4339e-05)
Transfer time : 394.78700994280655
Final mass    : 0.8344386642362569

# plots
fig1, ax1 = prob.plot_elements_history(to_keplerian=True, TU=TU/86400, time_unit_name="day")
fig2, ax2 = prob.plot_trajectory_3d(interpolate=False, sphere_radius=6378/LU)
fig3, ax3 = prob.plot_controls()
plt.show()