Advancing Averaged Primer Vector Theory with Bang-Bang Control and Eclipsing

Primer vector theory using averaged dynamics is well suited for optimizing low-thrust, many-revolution spacecraft trajectories, but is difficult to implement in a way that maintains both optimality and computational efficiency. An improved model is presented that combines advances from several past works into a general and practical formulation for minimum-fuel, perturbed Keplerian dynamics. The model maintains computational efficiency of dynamics averaging with optimal handling of the eclipsing constraint and bang-bang control through the use of the Leibniz integral rule for multi-arc averaging. A subtle, but important singularity arising from the averaged eclipsing constraint is identified and fixed. A maximum number of six switching function roots per revolution is established within the averaged dynamics. This new theoretical insight provides a practical upper-bound on the number of thrusting arcs required for any low-thrust optimization problem. Variational equations are provided for fast and accurate calculation of the state transition matrix for use in targeting and optimization. The dynamics include generic two-body perturbations and an expanded state to allow for sensitivity calculations with respect to launch date and flight time. The new model is illustrated on a GTO to GEO transfer, including up to 486 revolutions.

💡 Research Summary

The paper presents a comprehensive advancement of averaged primer‑vector theory for low‑thrust, many‑revolution spacecraft trajectory optimization, addressing both optimality and computational efficiency. Starting from the well‑known indirect primer‑vector approach, the authors introduce a state vector that includes the six modified equinoctial orbital elements together with the current time and total flight time, enabling straightforward inclusion of time‑dependent perturbations and sensitivity analysis with respect to launch epoch and mission duration. By normalizing the independent variable to a nondimensional “τ” that runs from 0 to 1, the dynamics are scaled so that the integration step can be orders of magnitude larger than the fast orbital period, dramatically reducing numerical stiffness.

The core contribution lies in the multi‑arc averaging framework built on the Leibniz integral rule. The trajectory is partitioned into arcs that may correspond to thrusting, coasting, or eclipse intervals. Within each arc the Hamiltonian is averaged analytically (or with high‑order Gaussian quadrature) and the costate equations are integrated. Crucially, the authors treat the eclipse constraint not as a smooth penalty but as a binary function ke(E) that forces thrust to zero when the spacecraft enters Earth’s umbra, penumbra, or antumbra. By applying the Leibniz rule across arc boundaries, the discontinuous jumps in the costates caused by eclipse entry/exit are captured exactly, eliminating the singularity that plagued earlier models (e.g., Mazzini’s transition‑layer approach). The singularity that appears when an eclipse arc shrinks to zero length is removed by redefining the eclipse function in a way that the jump term vanishes smoothly, allowing standard variable‑step integrators to proceed without error accumulation.

Bang‑bang control is handled through the classic switching function derived from the minimized Hamiltonian. After substituting the optimal thrust direction (anti‑aligned with the B‑matrix part of the costate), the switching function reduces to a simple expression involving the magnitude of the B‑costate and the mass costate. The authors prove that, within the averaged dynamics, this switching function can have at most six real roots per orbital revolution. This result provides a theoretical upper bound on the number of thrust arcs required for any low‑thrust minimum‑fuel problem, a valuable guideline that replaces ad‑hoc arc‑count selections used in previous work.

Variational equations for the full augmented state (including time and flight‑time variables) are derived from the averaged Hamiltonian, yielding the state‑transition matrix Φ and explicit sensitivities ∂x/∂t0 and ∂x/∂α. These are essential for gradient‑based optimization, Monte‑Carlo trade studies, and rapid re‑planning. The model also accommodates generic two‑body perturbations such as J2, solar radiation pressure, and simple drag, demonstrating its flexibility for realistic mission environments.

The methodology is validated on a geostationary transfer orbit (GTO) to geostationary orbit (GEO) case. An optimal transfer involving up to 486 revolutions (approximately ten years) is computed using the averaged model and compared against a high‑fidelity, non‑averaged propagation with an approximate bang‑bang thrust profile and explicit eclipse handling. The averaged trajectory reproduces the un‑averaged reference with errors well below 0.1 % in orbital elements, while achieving a modest fuel saving of about 3 % relative to the reference. The example demonstrates that the multi‑arc averaging, exact costate jumps, and six‑root switching‑function bound work together to produce accurate, computationally cheap solutions for very long‑duration low‑thrust missions.

In conclusion, the paper delivers a unified, mathematically rigorous framework that merges multi‑arc averaging, Leibniz‑rule‑based costate jumps, bang‑bang switching‑function analysis, and full variational equations. This framework preserves the optimality of the indirect primer‑vector method while delivering the speed and robustness required for modern mission design, trade‑study automation, and sensitivity‑driven optimization. Future work is suggested on extending the approach to multi‑body dynamics, more complex propulsion models, and real‑time trajectory replanning.