Tudat / tudat

NOTE: This Tudat version is no longer supported. See https://docs.tudat.space/en/stable/ and https://github.com/tudat-team/tudat-bundle for the new version
BSD 3-Clause "New" or "Revised" License
87 stars 142 forks source link

Counter for number of function evaluations during variable step size propagation #375

Closed mfacchinelli closed 6 years ago

mfacchinelli commented 6 years ago

Hi Dominic,

I am starting to evaluate the performance of the new USM propagators, but it would be nice to use the number of function evaluations (so also including discarded step sizes), instead of the computation time. Is this feature implemented in Tudat?

Have a nice day, Michele

mfacchinelli commented 6 years ago

Also, I have been using the Gauss propagator with Keplerian elements in my analysis, and I noticed it that with relatively 'loose' propagation settings, like these:

boost::make_shared< RungeKuttaVariableStepSizeSettings< > >(
                                rungeKuttaVariableStepSize, simulationStartEpoch, 100.0,
                                RungeKuttaCoefficients::rungeKuttaFehlberg56, 1e-1, 1e5, 1e-7, 1e-7 );

the orbit diverges very quickly. Is this the expected behavior? I have attached a figure showing the behavior of the step size over time, for various propagators. time_steps

DominicDirkx commented 6 years ago

Hi Michele,

For your first point: sure, this is not too hard. Add a member variable 'int numberOfFunctionEvauations_' to the DynamicsStateDerivativeModel class, and increment it by one upon every call to the computeStateDerivative member function. Write a get function that retrieves this int, and you whould have what you want. In case you want to reuse the same DynamicsSimulator object several times (or would like to merge this feature into the Tudat master), there should also be some functionality to ensure that the counter is reset to 0 at the start of every propagation.

I hope this helps, let me know if anything is unclear, or if you get stuck with it.

Dominic

DominicDirkx commented 6 years ago

Concerning your second question, at first I was quite mystified by it, but I think I now get what's going on: When propagating Kepler orbtis,5 out of the 6 elements change very slowly, and often very little compared to the their absolute values (depending on your initial elements). It may be that the variations in these 5 elements is all smaller than the integrator settings. In this case, for the purposes of your integrator, these variables 'don't change', and it can take a very large time step. The 6th element is the mean anomaly, and will change linearly with time (plus small variations). The linear behaviour is captured at any time step, and the small variations may have the same effect as small variations in the other 5 elements.

Could you give some details on the simulation settings you are using? Also, the Kepler element history from one of the more accurate propagations should be usable to understand whether the above explanation is indeed the possible issue.

Dominic

mfacchinelli commented 6 years ago

Hi Dominic,

Turns out that what you mention is already in Tudat. You can find the get function at line 488 of the file you mentioned (i.e., here). There is also a reset function, a couple of lines later. So now I should just add the resetFunctionEvaluationCounter( ) at the beginning of integrateEquationsFromIntegrator( ) in the integrateEquations.h file, and add a line ++functionEvaluationCounter_;in the propagators? Also, to retrieve the value, I will simply add a print statement at the end of the same integration function.

mfacchinelli commented 6 years ago

I am using a highly elliptical orbit as initial conditions (to test the propagators in an aerobraking-like scenario). The initial conditions are in Keplerian elements:

[27290000, 0.8375, 93, 158.7, 23.4, 180]^T

for a spacecraft around Mars. Also, the output obtained with a more accurate Cowell propagator, with integrations settings:

boost::make_shared< RungeKuttaVariableStepSizeSettings< > >(
                            rungeKuttaVariableStepSize, simulationStartEpoch, 10,
                            RungeKuttaCoefficients::rungeKuttaFehlberg78, 1e-3, 1e4, 1e-15, 1e-15 );

gives the plot attached (the dashed red line is the average Martian radius). cowell_reference

mfacchinelli commented 6 years ago

Ok, one more thing about the counter. I think that the counter is reset automatically by the dynamicsSimulator.h function (here). Also, I should probably call the dynamicsStateDerivative_->getNumberOfFunctionEvaluations( ) function from within this same file, so probably just after the integrateEquations( ) function. But if I want to add the counter increment, I should add the dynamicsStateDerivativeModel.h file in the header of the propagators?

mfacchinelli commented 6 years ago

Ok, one last comment. I re-read your first reply, and it appears that the counter is already incremented (here).

DominicDirkx commented 6 years ago

Huh, that's interesting. I wonder if I wrote it and forgot, or if Aleix added this when also adding the computation time

mfacchinelli commented 6 years ago

I have also added a new counter, to keep track of the function evaluations for each time step. This gives more insight in the integration process. I have pushed these changes to my repository (91816080434d3a17edf8dc2dad23d11039fca6a5). Thanks for the guidance!

mfacchinelli commented 6 years ago

Hi Dominic, regarding the Gauss propagator, I think it is as you said. I also tested the propagators with a near-circular orbit, and Gauss performs really well, both in terms of error and number of function evaluations.