Repository navigation
In formationBarycenter classical orbital-element averaging is singular near circular orbits #1656
Description
Describe the bug
With useOrbitalElements = True, FormationBarycenter::computeBaricenter() averages a, e and i linearly, and the argument of periapsis omega and the true anomaly f separately with the weighted circular mean atan2(sum m_i sin(alpha_i), sum m_i cos(alpha_i)).
In a near-circular orbit, omega and f of each spacecraft are not well defined individually. Only their sum, the argument of latitude u = omega + f, is. In a formation, the small eccentricity vectors of the members point in different directions, so the members' omega (and with them f) are spread over the whole circle while their u are almost equal. The circular means of omega and of f are then sums of vectors that nearly cancel, and their directions are arbitrary. The sum of the two means is not the mean argument of latitude. The scalar mean of e also ignores the directions of the eccentricity vectors.
The module documentation only lists the singularity of the ascending node for near-equatorial orbits as a limitation.
To reproduce
Three or four spacecraft of equal mass in a projected circular orbit formation of radius 100 m about a circular 7000 km orbit (i = 0.9 rad), without a spacecraft at the center. The members have e = 7.1e-6 and omega of, for example, 237.3, 117.3 and 357.3 deg:
import numpy as np
from Basilisk.architecture import messaging
from Basilisk.fswAlgorithms import formationBarycenter
from Basilisk.utilities import SimulationBaseClass, macros, orbitalMotion as om
mu = om.MU_EARTH * 1e9 # [m^3/s^2]
oe = om.ClassicElements()
oe.a, oe.e, oe.i, oe.Omega, oe.omega, oe.f = 7000e3, 0.0, 0.9, 0.3, 0.0, 1.0 # [m, -, rad, rad, rad, rad]
rc, vc = (np.array(x) for x in om.elem2rv(mu, oe))
n = np.sqrt(mu / oe.a**3) # [rad/s]
def barycenter(phasesDeg, useOrbitalElements, rho=100.0): # rho [m]
sim = SimulationBaseClass.SimBaseClass()
sim.CreateNewProcess("p").addTask(sim.CreateNewTask("t", macros.sec2nano(1.0)))
module = formationBarycenter.FormationBarycenter()
module.useOrbitalElements = useOrbitalElements
module.mu = mu
keep = []
for ph in np.radians(phasesDeg):
rho_H = [rho / 2 * np.cos(ph), -rho * np.sin(ph), rho * np.cos(ph)] # [m]
rhoP_H = [-rho / 2 * n * np.sin(ph), -rho * n * np.cos(ph), -rho * n * np.sin(ph)] # [m/s]
rd, vd = om.hill2rv(rc, vc, rho_H, rhoP_H)
keep.append(messaging.NavTransMsg().write(messaging.NavTransMsgPayload(r_BN_N=list(rd), v_BN_N=list(vd))))
keep.append(messaging.VehicleConfigMsg().write(messaging.VehicleConfigMsgPayload(massSC=100.0)))
module.addSpacecraftToModel(keep[-2], keep[-1])
rec = module.transOutMsg.recorder()
sim.AddModelToTask("t", module)
sim.AddModelToTask("t", rec)
sim.InitializeSimulation()
sim.ConfigureStopTime(0)
sim.ExecuteSimulation()
return np.array(rec.r_BN_N[-1])
for phases in ([0, 120, 240], [0, 120.1, 240], [0, 121, 240], [0, 125, 240], [0, 140, 240], [0, 90, 180, 270]):
rCart = barycenter(phases, False)
rOe = barycenter(phases, True)
print(f"phases {phases}: |r_cartesian - r_ref| = {np.linalg.norm(rCart - rc):.2f} m, "
f"|r_OE - r_cartesian| = {np.linalg.norm(rOe - rCart):.1f} m")It prints
phases [0, 120, 240]: |r_cartesian - r_ref| = 0.00 m, |r_OE - r_cartesian| = 1098382.5 m
phases [0, 120.1, 240]: |r_cartesian - r_ref| = 0.06 m, |r_OE - r_cartesian| = 17703.2 m
phases [0, 121, 240]: |r_cartesian - r_ref| = 0.63 m, |r_OE - r_cartesian| = 1777.4 m
phases [0, 125, 240]: |r_cartesian - r_ref| = 3.16 m, |r_OE - r_cartesian| = 367.6 m
phases [0, 140, 240]: |r_cartesian - r_ref| = 12.40 m, |r_OE - r_cartesian| = 106.1 m
phases [0, 90, 180, 270]: |r_cartesian - r_ref| = 0.00 m, |r_OE - r_cartesian| = 895148.1 m
When a spacecraft is at the center of the formation (chief plus three deputies), its exactly circular orbit dominates the means and the result is close to the Cartesian barycenter. The same holds for two spacecraft of equal mass, where the two circular means happen to add up to the mean u.
Expected behavior
For a formation whose members are close to each other, the orbital-element barycenter is close to the members, for any eccentricity, including near-circular orbits. If a configuration is singular for the chosen elements, the module documentation says so.
Screenshots
Na
Desktop (please complete the following information):
- OS: Linux Ubuntu
- Version: Basilisk 2.13.0b0
- Python version: 3.11.10
Additional context
Possible fix: average elements that are not singular for circular (and equatorial) orbits, for example the equinoctial elements of orbitalMotion (clElem2eqElem / eqElem2clElem): a, P1 = e sin(Omega + omega), P2 = e cos(Omega + omega), Q1 = tan(i / 2) sin(Omega) and Q2 = tan(i / 2) cos(Omega) linearly, and the mean longitude l (or the true longitude L) with the circular mean. This also removes the documented equatorial singularity. Averaging the mean anomaly (or mean longitude) instead of the true anomaly would keep the barycenter on the mean motion of the formation for eccentric orbits.
Metadata
Metadata
Assignees
Labels
Type
Projects
- Status
Activity
carlo98 commentedon Nov 1, 2026
Hi @zakariyaghazzol, I am not. Thank you 💪
zakariyaghazzol commentedon Nov 1, 2026
Hey @carlo98, are you already working on a fix for this? If not, I can investigate the equinoctial approach you suggested, starting with reproducing your example. Could we also clarify how the averaged orbit should behave for noncircular orbits, so the fix covers those cases as well?
zakariyaghazzol commentedon Nov 1, 2026
@carlo98, in that case. I’m leaning toward the mean-longitude option you mentioned. For spacecraft with the same semimajor axis, it keeps their shared mean motion in the two-body case. Is that the behavior you had in mind for eccentric orbits?
carlo98 commentedon Nov 3, 2026
Yes, that’s the behavior I had in mind.
Two additional notes to consider:
Consider also the use of the retrograde-factor variant of the equinoctial elements (to deal with i=pi).
In the documentation add a note on how e>=1 is dealt with.
schaubh commentedon Nov 3, 2026
We don't want to change the existing behavior of this module. If we just switch to doing mass averaged equinoctial elements, then we get a different OE-barycenter than if we do this using classical elements. Averaging classical elements just doesn't make sense near their singularities.
We could introduce a flag to do the orbit element averaging using a range of orbit element types. This keeps the existing classical element averaging, but introduces alternate equinoctial averaging if operating near equatorial or circular orbits.
schaubh commentedon Nov 4, 2026
@zakariyaghazzol, let me know if you are already working on this. Otherwise I can knock this out readily.
zakariyaghazzol commentedon Nov 4, 2026
Yes @schaubh, I’ve got it working locally with the optional setting you suggested. Just finishing up the review before opening a PR.
[AVSLab#1656] Clarify formation barycenter equinoctial option names
[AVSLab#1656] ensure the unit tests is robust to work across the comp…
[AVSLab#1656] Add multi-spacecraft equinoctial singularity regressions
[AVSLab#1656] Strengthen classical formation barycenter compatibility…
[AVSLab#1656] Verify runtime orbital selector transitions with distin…
[#1656] Add optional equinoctial formation averaging
[#1656] Clarify formation barycenter equinoctial option names
[#1656] ensure the unit tests is robust to work across the computers.
[#1656] Add multi-spacecraft equinoctial singularity regressions
[#1656] Strengthen classical formation barycenter compatibility tests
[#1656] Verify runtime orbital selector transitions with distinct for…