State transition matrix and mass

Hi folks,

I can’t find a native way of including the mass in the computation of the STM (with NumericalPropagator).
If confirmed, I think it’s a feature we need.

Cheers,
Romain.

Hello,

One would need to field everything so it should be simpler than that :sweat_smile:.

Cheers,
Vincent

Hi
It would be useful to include also the “longitude” (time dependent?) in the STM matrix; in the GTO transfer orbit, in my case, a cost function is defined that depends on the satellite’s longitude in the fixed terrestrial reference system (the same state vector with a different date at the end of the boost sequence generates different costs); in my experiments I added to the Jacobian matrix the derivative of the mass (calculated through the flow rates of the maneuvers) and of the time (calculated based on the stop time of the last apogee maneuver); for the mass everything works and the derivative of the cost function is well represented, while I cannot determine the gradient of the part of the cost that depends on the longitude.

1 Like

Yeah, that’s what I mean, it’s possible with FieldNumericalPropagator, but not the embedded STM in NumericalPropagator. It’s a bit weird actually. I’ll look how we could do this without breaking APIs.

@roccafrancesco I’m not sure to follow what you mean with the longitude. It’s a constraint on your final state, so it’s neither a state nor a parameter of the propagation model, therefore I don’t really see how it could be included in the variational equations. However with a full Field computation, there’s more applications (but it’s also slower).

Cheers,
Romain.

by longitude i mean the longitude of the spacecraft a the end of the last boost; the cost function include inclination cost,eccentricity cost and a “randevou cost”; the last component depend diereclty by longitude of spacecraft; By the way what i see is adding only derivative of the final time in the jacobian is not enough; the complete gradient is a combination of a cost function gradient respect to final state and the jacobian calculated with STM matrix; i wonder if using gast function directly can help; in this way maybe i dont need any transformation in Earth Reference frame and express the longidude derivative directly using cartesian element and time at the end of boorst;

I’ll merge today this feature.
It’s optional, by default the mass is not included in the STM.

As of now, there is a few limitations:

  • it is not compatible with all native parameter drivers. For example the TriggerDate
  • only Maneuver inheritors are considered. If you have a custom ForceModel modifying mass, it won’t see it because at the moment there is no direct access to the mass equation. To change that we would need to expose some sort of getMassDerivative in ForceModel. I think it could make sense. By default it would return 0

I saw also you added SphericalConstantThrustPropulsionModel that is the classes that i needed to implement for my GTO transfer Orbit optimezer; now i i included this class instead to use mine; the result in my case doesnt’ change for the gradiente calculation of cost function respect to the control parameter; i use azimuth,declination,duration,median as controls and im trying to calculate this cost function using the rule chain; what i do is calculate the gradient of cost function respect to orbital parameter ( keplerian,cartesian or equitoctial) and assemble them with the jacobian calculated by orekit; while the finite difference gradient is calculated whell and bring me to reach a solution for my optimizer; the calculation of this gradient in mixed mode doenst work; there is something i miss?

At the moment the STM with mass is not compatible with drivers like Duration and MedianDate of maneuvers. This is something i’ll try fixing in the coming days or weeks.

Cheers,
Romain.

it is ok; i mean i added the derivate of the mass manually and works; problem seem is composing analytic and numerical derivate; i tried also with direction derivate composition but doenst works; here the code;

/**

  • Calcola ∇h(x) con h(x)=f(g(x)) usando differenze finite direzionali

  • D_{v_i}f con v_i = colonna i di J_g(x).

  • @param extendedJ Jacobiana J_g(x) in R^(m×n)

  • @param stateVector vettore u = g(x) in R^m

  • @param costF funzione f: R^m → R (black-box)

  • @return gradiente ∇h(x) in R^n
    */
    private RealVector computeGradientViaDirectionalFD(
    RealMatrix extendedJ,
    RealVector stateVector,
    TwiceDifferentiableFunction costF
    ) {
    int n = extendedJ.getColumnDimension();
    RealVector grad = new ArrayRealVector(n);

    // 1) passo ottimo per differenze finite centrali:
    double h = FastMath.cbrt(Precision.EPSILON); // ≃ ε^(1/3)

    for (int i = 0; i < n; i++) {
    RealVector v = extendedJ.getColumnVector(i);

     // se v=0, ∂h/∂x_i = 0
     if (v.getNorm() == 0.0) {
         grad.setEntry(i, 0.0);
         continue;
     }
    
     // 2) spostamento Δu = h * v
     RealVector step = v.mapMultiply(h);
    
     // 3) valutazioni in u ± Δu
     double fPlus  = costF.value(stateVector.add(step));
     double fMinus = costF.value(stateVector.subtract(step));
    
     // 4) differenza centrale / (2h) → ∂h/∂x_i
     double derivative = (fPlus - fMinus) / (2.0 * h);
     grad.setEntry(i, derivative);
    

    }

    return grad;
    }

From the beginning the problem is always the same: the gradient with respect to the MEDIAN date; below the calculation of the gradient with respect to the parameters in 4 maneuvers ordered by Azimuth, Declination, Duration, Median; I calculated the gradient with finite differences and using Orekit via directional derivative; the only columns that have completely different values ​​are related to the MEDIAN date;

//FINITE GRADIENT
SQP grandient J{-0,3053850512; -32,5952346797; -422,0171079513; 160,8337792221; 0,5295239346; -17,2739674968; -261,1395807918; 176,4341718293; 19,4689659562; -4,9075625861; -31,0909050196; 120,7768057806; 2,5813631142; -0,8661122562; -3,148477395; 11,863828561}

//OREKIT
SQP grandient J{-0,3055421189; -32,5936577621; -418,6772621055; 945,3938704497; 0,5294236282; -17,2753028409; -295,1458682752; 955,7209276277; 19,4689878688; -4,9064925672; -49,2705031938; 486,7539394584; 2,5813626028; -0,8660034394; 0,3192549127; 81,2042592137}

Beware that adaptive step integrators increase noise in finite differences.

What triggers your maneuvers? Dates only or some state-based conditions?

i just added to my propagator builder a Maneuver class composed by DataBasedManeuverTriggrs and SphericalConstantThrustPropulsionModel

SphericalConstantThrustPropulsionModel propModel1=new SphericalConstantThrustPropulsionModel(isp, thrust,
new Vector3D(alfa,beta), “ORM1”) ;
DateBasedManeuverTriggers dataBase1 = new DateBasedManeuverTriggers(“ORM1”, ORM1_START, dur1);
Maneuver ORM1= new Maneuver(att,dataBase1 ,propModel1);

Ok thanks.
And you say that Median doesn’t work but Duration does?
It’s really weird because the two parameter drivers use pretty much the same code internally.

this is what i see comparing the gradients calcutated in the two ways i explaned above ;by the way also the gradient respect to the duration seem little bit different

i also tried with fixed step integrator and i removed the handler for retriving intermediate spacecraftState but the result doesn’t change; continue to have discrepancy between values of finite gradiente and gradient obtained by orekit jacobian;

I’m actually a bit puzzled too. I’m trying to enable TriggerDate with 7x7 STM and comparing with finite differences (see here commented out assertion). I’m also getting large differences for MedianDate. What are we missing here @luc?