Re: [Orekit Users] Fwd: Problem with additional equations

[MAISONOBE Luc <luc.maisonobe@c-s.fr>]

Christophe Le Bris <chris.lebris@gmail.com> a écrit :

> Hello,

Hi Christophe,

> I tried to reproduce the behavior exposed in my previous email with a
> minimal test.
> I integrate numerically a trajectory with only central force and an
> additional equation that increases the value of an additional state of
> 1/sec. The "Additional equation" should also increase the SMA of 2m/sec
> but, after integration for the duration of one keplerian orbit, my
> additional state changed to the value of the initial orbital period but the
> orbit did not changed => the increase of 2m/sec has not been taken into
> account.

I confirm the behaviour you get, and also confirm your thorough analysis.
The latest version on Hipparchus indeed does not allow the secondary
equation to modify the main equation.

However, I am on the fence with it. Should we really allow an additional
equation to modify the main equation? If I remember well, the fact that the
derivatives of the main state are provided as an argument is mainly because
they may be needed to compute the secondary derivatives, it was not because
users could modify the primary derivatives.

This was however not obvious in the Apache Commons Math era because the
signature of the method was anyway to provide arrays and have the user
fill them. So the interdace SecondaryEquations did specify the
computeDerivatives method as:

  /** Compute the derivatives related to the secondary state parameters.
   * @param t current value of the independent <I>time</I> variable
   * @param primary array containing the current value of the primary  
state vector
   * @param primaryDot array containing the derivative of the primary  
state vector
   * @param secondary array containing the current value of the  
secondary state vector
   * @param secondaryDot placeholder array where to put the  
derivative of the secondary state vector
   * @exception MaxCountExceededException if the number of functions  
evaluations is exceeded
   * @exception DimensionMismatchException if arrays dimensions do  
not match equations settings
   */
  void computeDerivatives(double t, double[] primary, double[] primaryDot,
                          double[] secondary, double[] secondaryDot)
      throws MaxCountExceededException, DimensionMismatchException;

As you see, primaryDot and secondaryDot are passed the same way,  
despite only the
secondaryDot is really intended to be filled by user.

Since Hipparchus, the design is that arrays the user should fill must be
return values (i.e. allocated within the method), so the equivalent interface,
which is now called SecondaryODE defines:

    /** Compute the derivatives related to the secondary state parameters.
     * @param t current value of the independent <I>time</I> variable
     * @param primary array containing the current value of the  
primary state vector
     * @param primaryDot array containing the derivative of the  
primary state vector
     * @param secondary array containing the current value of the  
secondary state vector
     * @return derivative of the secondary state vector
     * @exception MathIllegalStateException if the number of  
functions evaluations is exceeded
     * @exception MathIllegalArgumentException if arrays dimensions  
do not match equations settings
     */
    double[] computeDerivatives(double t, double[] primary, double[]  
primaryDot, double[] secondary)
        throws MathIllegalArgumentException, MathIllegalStateException;

It makes more clear than only the derivatives of the secondary state  
vector are expected
from this method. The primaryDot array is still passed because it may  
be needed in
to compute the other derivatives, but is was not intended that users  
change it.

For the sake of separation of concerns, I would prefer to have realy  
the primary
equations completely computed by one part of the code, and secondary equations
by another part and not messing with the primary. So if you remove the  
secondary
equations, you only don't get the secondary state they compute, but  
you do not change
at all the premary equations. When the secondary equations do change  
the primary
derivatives, they are not secondary anymore: they are needed to compute the
primary state.

I understand that moving the mapper.insertEquationData call  
corresponding to the main
state from before the computation of the secondary equations to after  
the computation
allows as a side effect changes made by the secondary equations to be  
taken into
account, but it is not really clean. It implies that the primary  
equations were
indeed incomplete, and it implies the secondary equation returns two  
different results,
on in its regular return value, and one by changing a value in an  
array passed by the caller,
and expected this array to be used only after the call, not before.

In your use case, would it be possible to split your computation in  
two parts, one in
the secondary equation (to compute the return value as in your return  
new double[]
{ 2, 0, 0, 0, 0, 0, 0 } example) and another one as a force model (to  
compute the
main state derivative as in your pDot[0] = 1 example)?

best regards,
Luc

> *Initial orbit = Keplerian parameters: {a: 7800000.0; e: 0.001; i: 28.0;
> pa: 0.0; raan: 0.0; v: 0.0;}Initial keplerian period = 6855.717044697212
> secFinal value of added state = 6855.717044697215Final orbit = Keplerian
> parameters: {a: 7800000.000000003; e: 0.001000000000000334; i:
> 28.000000000000007; pa: -1.1173044601630712E-25; raan:
> 2.925205344836198E-29; v: -5.234453515214931E-13;}*
> I expect the problem to be in Hipparchus 1.2 (at least) and not Orekit. In
> *org.hipparchus.ode.ExpandableODE*, I think the line 137
> "mapper.insertEquationData(index, primaryStateDot, yDot);" should be placed
> just before the "return".
>
>
> If I execute the same program with Orekit 7.0/Commons Math 3.4.1, the
> behavior is different and is ok:
>
>
>
>
>
> *Initial orbit = keplerian parameters: {a: 7800000.0; e: 0.001; i: 28.0;
> pa: 0.0; raan: 0.0; v: 0.0;}Initial keplerian period = 6855.717044697212
> secFinal value of added state = 6855.717044697215Final orbit = keplerian
> parameters: {a: 7813711.434089395; e: 0.0010000000000002728; i:
> 27.999999999999993; pa: -1.3080031127518713E-13; raan:
> 3.7263393091775305E-17; v: -0.4748813935419053;}*
>
> Thank you!
>
> Christophe
>
>
> ---------- Forwarded message ----------
> From: Christophe Le Bris <chris.lebris@gmail.com>
> Date: 2017-11-30 18:21 GMT+01:00
> Subject: Problem with additional equations
> To: "orekit-users@orekit.org" <orekit-users@orekit.org>
>
>
> Hello,
>
> I have a problem when integrating "additional equations" that impact the
> derivatives of the primary state (ie: value returned by "computeDerivative"
> method is not null).
> I have the feeling that these values are not taken into account during the
> integration.
>
> I use Orekit 9.1 with Hipparchus 1.2.
> It seems that the origin of my problem is located in Hipparchus library (in
> the "computeDerivative" method of  the class "ExpandableODE") and not in
> Orekit but I would like to know if someone already used this functionality
> with the last versions of Orekit/Hipparchus ?
>
> Thank you.
>
> Christophe