Modeling Hohmann and Bielliptic maneuvers

Orekit Python Wrapper usage
1

Hi all,

I recently started using the Orekit Python wrapper and wanted to model an impulsive Hohmann trasfer. After going trough different examples in the official wrapper repository and checking different code snippets in this forum, I was able to understand the logic behind how ImpulsiveManeuver instances are applied.

However, I am not sure how can I define the deltaV vectors aligned with local velocity in both periapsis and apoapsis of the orbit. The following code shows my current implementation of the Hohmann impulsive transfer, which is expected to output two instances of the ImpulsiveManeuver class, being the first one the apogee raising and the second one for orbit circularization:

def hohmann(k, r0_vec, v0_vec, r_f, ISP=300):
    """ Computes required impulses for a Hohmann transfer """

    # Initial position and velocity magnitude
    r0_norm, v0_vec = [norm(vec) for vec in (r0_vec, v0_vec)]
    a_trans = (r0_norm + rf_norm) / 2

    # Compute the magnitude of Hohmann's deltaV vectors
    dv_a = np.sqrt(2 * k / r0_norm - k / a_trans) - np.sqrt(k / r0_norm)
    dv_b = np.sqrt(k / rf_norm) - np.sqrt(2 * k / rf_norm - k / a_trans)

    # Transform previous magnitudes into vectors aligned with local velocity
    # in the local orbital frame
    ?? ?? ?? ??

    # Define the events which trigger the impulses
    at_periapsis = PositionAngleDetector(orekit_ss, PositionAngle.TRUE, 0)
    at_apoapsis = PositionAngleDetector(orekit_ss, PositionAngle.TRUE, np.pi)

    # Build the impulses
    imp_A = ImpulseManeuver(at_periapsis, dVa_vec, ISP)
    imp_B = ImpulseManeuver(at_apoapsis, dVb_vec, ISP)

    return imp_A, imp_B

I know there is a LOFType class, but I am not sure how to properly add it to previous algorithm… On the other hand, a similar procedure could be applied for Bi-elliptic transfers.

Thanks for this amazing software and its open-source nature!

2

Hi @jorgepiloto,

Welcome to the Orekit forum!

You indeed need to use a LOFOffset attitude provider.
Given the inertialFrame you’re using.

attitudeProvider = LofOffset(inertialFrame, LOFType.TNW)
imp_A = ImpulseManeuver(at_periapsis, attitudeProvider , dVa_vec, ISP)

With dVa_vec = dv_a * [1.0, 0.0, 0.0]

That way your DV will be along velocity; since TNW local orbital frame has its X axis along velocity (the “T” stands for “Tangential”).

Beware that I’m not an expert in Python and I’ve written the code without testing it; so you may need to tweak that example to make it work.
But this should give you a general idea of what you need to do.

Maxime

3

Sorry for this late response @MaximeJ! After trying you suggestion I was able to finally apply an impulsive maneuver to a KeplerianOrbit instance.

However, after checking against GMAT2020a and poliastro (both open-source orbital mechanics software), I found final results not to match the ones claimed by previous software. Let me attach my full current implementation. Only numpy and orekit packages are required for running it.

Even after forcing vectors to be in the GCRF frame, final position and velocity vectors are not the same. Not sure if I am missing something…

4

Hi @jorgepiloto,

I ran your script and you do indeed get quite some difference in the end status (for others to see):

Assert final position
rf_vec = [18234803.003842704, -1089464.5055657763, -0.0]
rf_expected = [-18447181.34824541, -30659530.26513587, 2.116148304903275e-10]

Assert final velocity
vf_vec = [332.64926224347874, 3899.2434203912685, 0.0]
vf_expected = [2859.891379982459, -1720.738617222231, 1.600645301237219e-14]

Cannot spot anything directly in the code (nice code!), have you double checked the magnitudes of the maneuvers with other systems? Are the coordinate systems right?

Regards

5

Hi again,

You could also double check that your orbit definition is same in all systems, your PV is converted to a keplerian orbit:
<KeplerianOrbit: Keplerian parameters: {a: 8531260.567751326; e: 0.1560450014601089; i: 0.0; pa: 0.0; raan: 0.0; v: 0.0;}>

6

Hi @petrus.hyvonen,

thanks for your reply. After checking units and frames the divergence of the results is still present. Let me attach the GMAT script just in case anyone wants to give it a try:

7

Hi,

I haven’t got GMAT up and running, but I see that it is using a numerical propagator and your orekit the analytical keplerian - however I do not expect that to give as large differences as you get.

My gut feeling is that it is something in the orbit definition/time/coordinate systems, often is… :slight_smile: Note that the both the time and coordinate system definition differs in the two codes - did you also note that GMAT is using some own definition of MJD:
http://gmat.sourceforge.net/docs/nightly/html/SpacecraftEpoch.html

8

Oh! I had no idea about GMAT’s own definition of MDJ… Thanks for the advice!

I am also suspecting about the fact that after the last impulse, the orbit becomes a circular one. For example, this is warned in the ApsideDetector events, see:

Beware that apside detection will fail for almost circular orbits. If for example an apside detector is used to trigger an ImpulseManeuver and the maneuver change the orbit shape to circular, then the detector may completely fail just after the maneuver has been performed!

9

Hi Jorge,

I had a look at both your Orekit and GMAT scripts. I tried several things but I haven’t managed to obtain the same results than GMAT either. Here are possibilities of explanation:

  1. Epoch definition

As Petrus mentioned, the epoch definition differs between both. To match GMAT, you could define the epoch this way in Orekit

utc = TimeScalesFactory.getUTC()
epoch_00 = AbsoluteDate(2000, 1, 1, 11, 59, 28.0, utc)

I checked that epoch_ff then matches the last timestamp in GMAT’s report: 01 Jan 2000 17:29:18.489 UTC.
But this had no effect on the results.

  1. Delta-V values

The Delta-V values computed manually in the Python script don’t match the converged values in GMAT.
I wrote the GMAT values in the Orekit script. This does not bring the same results as GMAT though.

  1. Equatorial orbit

The orbit you are using is nearly equatorial, it can bring numerical issues when using KeplerianOrbit. You should use EquinoctialOrbit instead. The constructor has the same syntax: EquinoctialOrbit (ORbit Extrapolation KIT 10.3 API)

This also requires a change in PositionAngleDetector, the orbit type argument must be changed to OrbitType.EQUINOCTIAL.

However, this does not have any effect on the results.

  1. Maneuvers not happening at the right time

I am far from being an expert in maneuvers in Orekit, so I am not sure about this behaviour:

The first maneuver does not happen at the first perigee (like 12 minutes after the epoch) as it does in GMAT, but at the second perigee instead (around 8800 seconds).

It might be because the true anomaly of the initial KeplerianOrbit (but it is also the case with an equinoctial orbit) is negative? I thought the true anomaly is normally in the 0-360° range?
I tried forcing the initial true anomaly to be positive, but it did not change anything, the first maneuver still happens at the second perigee.

You could take a look at the PositionAngleDetector::g method: src/main/java/org/orekit/propagation/events/PositionAngleDetector.java · develop · Orekit / Orekit · GitLab

Cheers
Clément

1 Like
10

Hi!
I also had a look and ran the orekit and gmat scripts and by following some suggestions here, I finally managed to get the same GMAT results.

However, as yzokras just mentioned, I used the deltaV values converged from GMAT, (as they differed by ~100m/s).
I also used ApsideDetectors instead.

Here are the modifications I did on your script Jorge, orekit_hohmann_mod.py · GitHub

1 Like
11

Hi @yzokras and @marcpm,

thanks for your feedback on the issue, it was very useful. Let me present some results in the following lines:

Delta-V computed values were wrong

As pointed out by @yzokras, dVs were not matching expected ones. After a deep research, I was able to find a bug in my current implementation. Notice the magnitude of the initial position and velocity vectors are computed as:

r0_norm, v0_norm = r0_vec.getNorm(), v0_vec.getNorm()

however, for Hohmann we need as initial state vectors the ones for peripasis! Therefore, previous line needs to be as follows:

a, ecc = ss_0.getA(), ss_0.getE() 
r0_norm = a * (1 - ecc)
v0_norm = np.sqrt(2 * k / r0_norm - k / a)

With this fix, the dVs were properly computed. However, even with this change, the maneuver still did not reach expected final state vector. In fact, the PositionAngleDetector was being triggered at wrong epoch, this still needs further investigation…

Using ApsideDetector brings right result

As @marcpm claimed, using the ApsideDetector for triggering the maneuver (once previous dVs bug has been fixed) made the final result to match expected one.

1 Like