ODE solver efficiency

[giovaniceotto]

The Flight class in rocketpyAlpha.py currently employs Scipy's LSODA wrapper to solve to dynamic equations of flight. However, it creates a new solver instance for every micro time step, which is highly inefficient.

The LSODA algorithm uses a variable time step Adams Bashforth implementation for non-stiff problems. Since Adams Bashforth is a multistep method, every time a new solver instance is created, the previous steps are lost and new data needs to be created using other algorithms to start the Adams Bashforth. Therefore, a lot of time is spent using other solvers instead of the Adams Bashforth itself. For very small time steps, it may be the case that none of the advantages of a multi step method are being exploited.

This inefficient implementation should be fixed either by using Scipy's LSODA class and its step method or by moving to a new ODE solver.

[giovaniceotto]

In the _ode_solver_improvement_ branch, I just submitted a commit which does exactly what was suggested.

Instead of creating a new Scipy _solveivp class instance at every macro iteration, a single LSODA solver instance is created every time the derivative function changes (such as when events happen). With this new implementation, every macro iteration runs several steps of the LSODA solver.

One drawback is the loss of the event detection which is natively supported by the _solveivp class. This has been substituted by simple event detection of the apogee event and the impact event after every step of the LSODA solver. To check for apogee, after every step the algorithm checks if its the first time that the rocket's vertical velocity has become negative. To check for impact, the same is done with the rocket's altitude.

This temporary event detection method is not very accurate since it does not employ any root finding. Root finding should be employed in order to get better results, specially for the impact event which can happen at high velocities when no parachute events are triggered.

Let's talk efficiency now. Using parachute trigger functions running at 100 Hz, a typical high powered rocket flight would need 40k derivative function evaluations when the old implementation was used (considering no parachute actually opened). With the new implementation, this number decreased to about 20k function evaluations. In average, every step runs 2 iterations of LSODA's Adams Moulton algorithm. All of these numbers were accessed using absolute and relative tolerances of 1e-6.

Since the derivative function evaluation is still the major bottleneck of the code, the execution time was cut in a little more than half, which is a promising result.

[giovaniceotto]

Currently the rail phase of the flight is being simulated using VODE instead of LSODA. Furtheremore, it is a completely separated implementation from the one which simulates every other flight phase, using LSODA. To make matters worse, the detection algorithm for the event which characterizes when the rocket leaves the rail is very naive and makes no use of root finding.

One alternative would be to use Scipy's _solveivp class to integrate the equations of motion of this flight phase. This has the benefit of using LSODA and it supports a more sophisticated event detection. Since this flight phase does not check for any trigger functions, this method would be ideal.

[giovaniceotto]

The fact that the rail flight phase is separate from the others and is not checking for any event triggers might be a problem. Even though it is very unlikely, events might be triggered during rail phase and later on, time discrete thrust control systems may also need to operate during the rail phase.

Therefore, I would suggest merging the rail phase with the rest of the flight phases. This would invalidate the suggestion of using Scipy's _solveivp class to make use of its native event detection capabilities, but the same root finding event detection employed for the apogee and impact event can be used, forming a general and well established algorithm which is mature enough to deal with all flight phases and both discrete and continuous events throughout them.

[giovaniceotto]

A commit has just been submitted to the _ode_solver_improvement_ branch in which this last suggestion has been implemented.

Considering the current implementation in the master branch, these are the main changes:

• Numba has been removed temporarily, since it was taking too much time to compile. I offers an advantage when running the simulation for the exact same Rocket in the exact same Environment. Since this is rarely done for now, it was decided that Numba was not worth it. For future references, Numba does have a significant speed up (abou x2) after compilation is done. But for now, compilation takes about 10 seconds.

• The rail flight phase has been incorporated to the simulation with the other flight phases. Now, every step of the flight is simulated under a unifying framework which uses LSODA and takes full advantage of its variable step size multi step implementation.

• The unified simulating framework has been completely reformulated to allow for better handling of parachute trigger function and also discrete control functions which are expected to be implemented in the near future. Its works in 4 layers:

  • The top one is characterized by an instance of the new auxiliar class FlightPhases. It stores each phase of the flight (a FlightPhases.FlightPhase object), specially concerning derivative changes. It starts with two phases, one with the derivative for rail phase simulation and another one which just represents the maximum allowed simulation time. During the simulation, new flight phases are added when the rocket leaves the rail and when parachutes are triggered, both of these specifying a new derivative. Furthermore, parachute events also make use of the ability of specifying callback functions. These functions are evaluated before a new flight phase starts and parachute events use them to set the current parachute CdS.
  • The next layer is made out of an instance of the new TimeNodes auxiliary class. For each flight phase, time nodes (TimeNodes.TimeNode objects) by which the integration must pass during the time of the flight phase are saved to this instance. Each TimeNode object which is stored inside this instance has a specific time stamp, a list of parachutes whose triggers need to be feed and also a list of callback functions which must be called (at the time stamp).
  • Next is the actual LSODA layer which steps through time as necessary, going from the time stamp of the current TimeNode to the time stamp of the next TimeNode.
  • Finally, the last layer comprises of overshootable TimeNodes. Essentially, this is exactly the same as a TimeNodes object, however, it spans from the initial to the final time of a single LSODA integration step, instead of spanning the time of a FlightPhase. After the integration step is over, the code looks for all parachute triggers and callbacks which are overshootable, that is, which we can simulate past and then go back to them without too much cost associated with it. For parachutes running at high sampling rates but only raising a flight event once, it can be significantly faster to let the LSODA solver step as needed (instead of constraining its maximum step size to the time between two TimeNode objects of the TimeNodes instance) and then going back and feeding the triggers that should have been fed in the time spanned by the integration step. Discrete controllers would probably change something about the rocket at every call, so using them as overshootables might not be a good idea since every time a overshootable is triggered, the simulation is rolled back. This is the reason that discrete controllers are better suited as TimeNodes.

• Root finding has been implemented for detection of out of rail event, apogee event and also impact event. Out of rail and impact are solved by interpolating the last integration step using a cubic polynomial defined by the points and their derivatives. The roots of the polynomials are found using the analytic formula. In the case of the apogee, linear interpolation is used and the root is also found analytically. All methods provide good results with a accuracy much higher than needed, and they are also faster than using iterative root finding methods.

• Simulation accuracy can now be controlled much better and default settings allow for an accuracy of mm and mm/s in position and speed, respectively.

• Noise feed to parachute triggers have been relocated to the Rocket class. Now, when creating a parachute, its noise properties can be specified. This is a much cleaner and flexible implementation.

In terms of speed, using parachutes as overshootables allow for a simulation which runs in just under 900 ms, considering 100 Hz sampling rate and parachutes being ejected (about 240 s of simulation time). Turning overshooting off (this is controlled by a parameter of the Flight class), the simulation runs in just under 8 seconds, which shows how significant overshootables can be. The execution using the current implementation in the master branch clocks at 10 seconds.

The fact that the new implementation and the old one run at similar speeds when not using overshootables, even though the LSODA solver is performing about half the number of derivative function evaluations is really disturbing. The main reason for it is the use of the TimeNode object which is only used for code organization and readability and could very easily be substituted by a list. The creation of TimeNode objects during simulations takes seconds which is definitely not good. This could easily be fixed, but would be detrimental to code readability. Therefore, it will be left as is for now and overshootables will be preferred. When discrete time controllers are implemented, this may be considered again, perhaps.

Just as a side note, using Numba (second runs, since the first one is slow due to JIT compilation) is significantly faster, clocking in at just over 4 seconds, presenting a significant advantage. Since discrete time controllers might make more use of simulating the same Rocket, this may be a good alternative.

[giovaniceotto]

Seems like this issue has been solved! Marking it as closed.