Thrust and rotation strategy to circularize a standard GTO orbit using ion propulsion?

Question

The previous question How much time does it take to circularize a GTO orbit using ion propulsion? has a good, quantitative answer based on a known delta-v.

But I'm wondering how would you actually aim a continuous ion thrust as you move around in a highly elliptical orbit at the beginning? Is there a way to do this without wasting a lot of propellant? Or do you use the thrust only during a fraction of the orbit (in the beginning) when near periapsis, and therefore it takes longer to accumulate the delta-v because the thruster is actually turned off a lot of the time?

Here's a quick, silly calculation showing that for 0.2N and 2200kg, starting at 300km LEO, it takes about 600 days to slowly spiral out to GEO. The satellite rotates once per revolution around the earth, so I have to slow down the rotation rate a teeny tiny bit once in a while to keep them synchronized so the thrust is always tangent to the circular orbit.

---

Now suppose I want to start from a standard Hohman transfer ellipse (GTO), and run my program again. Periapsis altitude 300km, apoapsis at GEO. Let's assume there's nothing else at GEO to collide with to make it simple.

Q1: In order to use the least amount of ion propellant, should the spacecraft be slowly rotating, or have a fixed attitude, perhaps tangent to apoapsis? If rotating, at a constant, or ramped rate?

Q2: How would one modulate the thrust and attitude to circularize a standard elliptical GTO using the 0.2N ion propulsion in the least amount of time?

def deriv(X, t):

    r, v  = X.reshape(2, -1)
    vhat  = v / np.sqrt((v**2).sum())
    rhat  = r / np.sqrt((r**2).sum())
    rsq   = (r**2).sum()

    acc_grav   = -GMe * rhat  / rsq
    acc_thrust = vhat * acc   # right now, thrust is just in vhat direction

    return np.hstack((v, acc_grav + acc_thrust))


import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import odeint as ODEint

GMe = 3.986004418E+14  # m^3/s^2

re  = 6371.E+03   # m
alt = 300.E+03    # m
a   = re + alt    # m
v0  = np.sqrt(GMe / a)  # vis-viva equation for circular

T = 2. * np.pi * a / v0
print "period (min) ", T/60., " at ", alt, " km altitude."

acc  = 0.2/2200.  # m/s^2  0.2N and 2200kg

X0   = np.array([a, 0, 0, 0, v0, 0])  # initial state vector

days = np.linspace(0, 600, 100000)
t    = days * 24. * 3600.

tol  = 1E-10  # good enough for rough answer but not an accurate trajectory
answer, blob = ODEint(deriv, X0, t, rtol=tol, atol=tol, full_output=True)

xp, yp  = 0.001*answer.T[:2]
alti_km = np.sqrt(xp**2 + yp**2) - 0.001*re

if 1 == 1:
    plt.figure()
    plt.subplot(3, 1, 1)
    plt.plot(days, xp)
    plt.plot(days, yp)
    plt.title('x, y (km)')
    plt.subplot(3, 1, 2)
    plt.plot(xp, yp)
    plt.title('xy (km)')
    plt.subplot(3, 1, 3)
    plt.plot(days, alti_km)
    plt.title('altitude (km)')
    plt.show()

Answer 1

TL; DR: Trajectory optimization for continuous thrust is difficult and this field is very active in research.

2021 clarifications:

Methodology

1. For the least amount of fuel, the best is the thrust the least amount of time as possible and only when it's extremely efficient ($\eta \geq 0.98$). But that also implies that it will take an incredibly long amount of time to reach the desired orbit.
2. Conversely, to reach an orbit in the least amount of time, one should thrust as much as possible, even when it isn't that efficient.
3. Typically, in trajectory optimization, we set a maximum transfer time and optimize for the lowest total thrust over time (i.e. lowest overall $\Delta v$) in that maximum time. This allows a fine balance between "least amount of fuel" and "get this spacecraft to its mission orbit and start the science/making money".

In this older answer, I mention the Ruggiero control laws: these provide the most efficient thrusting angles based on true anomaly to change an orbital element (or set thereof). The paper also provides the partial derivatives of the equations that one may use these to determine the efficiency ($\eta$) of modifying an orbital element at a given true anomaly.

Answers

These answers are specific to raising the semi-major axis and reaching an eccentricity of zero. It does not answer the solution for a real GEO orbit, which has additional constraints on inclination (near zero) and other orbital elements specific to a given mission.

Q1: In order to use the least amount of ion propellant, should the spacecraft be slowly rotating, or have a fixed attitude, perhaps tangent to apoapsis? If rotating, at a constant, or ramped rate?

If your thrust is always in the direction of velocity vector, you're continuously increasing the velocity of the spacecraft and therefore you'll be increasing your SMA (along with modifying other elements). To circularize an orbit, you need to raise only your periapse until it matches your apoapse. To raise your SMA, you need to thrust all the time... but that will affect your eccentricity. Therefore, the best approach is probably to raise your SMA until GEO and then circularize your orbit. To consume the least amount of fuel, you should only turn on your engine when you're within a very small region of your periapse such that your apoapse is very slightly increased one half-orbit later. Repeat that until you've reached your target SMA. Then circularize.

Q2: How would one modulate the thrust and attitude to circularize a standard elliptical GTO using the 0.2N ion propulsion in the least amount of time?

One would be continuously thrusting and varying the in-plane and out-of-plane angles such that both the inclination and the SMA are modified at any given time throughout the orbit, but by weighting one of the orbital elements more then the other based on their efficiency. The high-fidelity Nyx astrodynamics toolkit provide the Ruggiero control laws out of the box.

Full details (2017)

Edit: Concerning the rotation of your spacecraft, you'll want to plot the in-plane and out-of-plane thrusting angles (with respect to the RNC frame of the spacecraft). That will give you an idea of how much the engines needs to gimbal by before thrusting. That angle depends entirely on what orbital parameters you're changing (cf. Naasz, Ruggiero and Petropoulos control laws). I expect that angle to not vary by much throughout an orbit correction, but I may be wrong (hopefully will have an answer soon). This expectation is based on the equations which perform the instantaneously optimal correction of each parameter, cf. figures 4 and 5 at the bottom.

If using continuous thrust, you'll generally want to solve the optimal control problem. When solving an optimal control problem, you're minimizing a functional, i.e. you're searching for a function (the control) which will minimize your cost function at each time step (cf. screenshot of my group's presentation below).

There are different ways to solve the optimal control problem for low thrust trajectory design, and (as often, sadly) each group of researchers is somewhat silos.

From the (unpublished) research that two colleagues and I did, there are a few strategies to solve this, which I summarize below. All but the first one solve an optimal control problem.

1. The simplest: use a sub-optimal closed-loop control law (like the Petropoulos Q-Law) which is fed to an optimizer like a genetic algorithm. This leads to a sub-optimal solution for the trajectory but is often very close to the optimal solution. All you need is to code up the Q-Law and the GA and define an initial and final orbit. Other controls worth checking out include Naasz and Ruggiero (but note that the initial formulation for the Naasz control was only validated on cases when you want to increase the value of an orbital parameter, and you need to perform a slight change in the equations to make it work the other way (I can provide the info if needed, I validated the changes in my own research)).

2. The most theoretical and math-heavy method: the indirect methods, which solve the two point boundary value problem by finding the Hessian. Requires the initial and final conditions desired. There are several such methods, mostly started by Betts in 1998, including Direct Shooting, Indirect Shooting, and used more recently, indirect collocation. One of the issues you'll encounter is the selection of the initial Lagrange multipliers, knowing that a poor choice may prevent convergence. In addition, these methods lead to open loop solutions.

3. Direct methods: time is discretized and at each step you're trying to find the best control. This is also an open loop method, and because of how the problem is posed, you encounter the curse of dimensionality (cf. Lantoine 2012). Yet, this is considered by some as the state of the art of trajectory optimization (cf. "Spacecraft trajectory optimization", Conway 2010).

4. A better way: dynamic programming: the idea here is the split up the problem into simpler independent problems (a "divide and conquer"-like approach) and apply the principle of optimality (Bellman's graph traversing principle where the shortest path from a to c corresponds to the shortest path from a to b and from b to c) (cf. figure 2). The main drawback here is that splitting up the state space into independent sub problems means that you're going to be using a ton of memory to store the ongoing solutions... and then you'll need to traverse this graph. The main advantage though is that you'll always find the global optimal solution, and the returned functional is a closed loop control.

5. Differential dynamic programming: this isn't a very new method per se, but it's definitely gaining tracking again in recent history. Specifically, it is used in NASA Mystic which does the trajectory optimization for the Dawn mission (which has three Safran PPS-1350 but only one turned on at each time, it's a pretty cool mission). The idea of this method is to pick an initial guess for the control, propagate the (real) dynamics forward, compute the cost, perform a correction of the control by propagating backwards (but using a quadratic approximation of the dynamics (i.e. you'll need the state transition matrix and the Jacobian of the STM too, which is the Hessian of the dynamics)), and finally check if you're minimizing the cost function by propagating forward again (cf. figure 3). The references you'll want to check for this method include Liao 1991, Pantoja 1988, Lantoine 2012, Ozaki 2015, and eventually Aziz 2016 (for an explanation of Hybrid DDP).

I hope this helps!

Figure 1

Figure 2

Figure 3

Figure 4

Figure 5

Answer 2

I originally posted this answer here. This is a figure that I have from a class assignment from a few years back. While definitely not a practical trajectory, it shows the characteristics of how to transfer from an elliptical to a circular orbit. This solution was computed using indirect optimization. This problem assumed constant thrust magnitude (so the thruster is always firing), with the control variable being the thrust angle. The objective function was min(t_f).

As the figure shows, the thrust vectors (yellow arrows) are approximately in the velocity direction near periapsis and apoapsis, and they are approximately perpendicular between the apses. In general, when the thrust vectors are in line with the velocity vector (parallel or anti-parallel), they are adjusting the orbital energy (and therefore increasing/decreasing semi-major axis). When the thrust vectors are perpendicular to velocity, they are adjusting the eccentricity (shape) and argument of periapsis (orientation) of the orbit.

This problem is 2D, but to achieve an inclination change you will want to thrust out of plane near the apses. You could visualize this as pointing "into the page" at periapsis, and slowly rotating up until it's in-plane at the mid-point, and finally "out of the page" at apoapsis.

A more realistic transfer from GTO to GEO would use MANY more spirals than in this problem (which is a goofy transfer from some crazy elliptical orbit to Mars' orbit). For a CubeSat or smallsat with very limited propulsion, a transfer around 6 months could take over 250 spirals!

Another notable distinction to make is that this problem solved a minimum-time as opposed to a minimum-fuel transfer. In general, minimum-fuel problems are notoriously much more difficult to solve using indirect optimization, and it is common to use homotopy to first solve the problem min(u^2) (where u is the control effort, in this case the thrust magnitude), and then with fancy math step closer and closer to min(u). The min(u^2) problem has a much wider radius of convergence so successful initial guesses of the costates are much easier to find. These numerical difficulties are part of the reason why direct methods are almost always preferred in practice. For very-low-thrust transfers, a minimum-time transfer can often be close to a minimum-propellant transfer.

When you have high thrust available, a GTO to GEO transfer would be completed with impulsive maneuvers at the apses. As your max thrust level decreases, you will have to either thrust over larger sections of the orbit, you will have to thrust for more orbits, or both. Available propellant and time constraints will drive the final trajectory selection.