Simulation stability

[sanguinariojoe]

Right now the simulation stability completely rests on the user's selected time step. For instance, let's reduce the time step in example_005.ipynb just replacing
N = tfin * 100 + 1
by
N = tfin * 10 + 1
and enjoy the arising instabilities at the end of the simulation (you may also increase tfin to 200 seconds if you wanna crash your plane :-) )

I don't know much about flight dynamics, but probably a Courant condition can be optionally imposed. Right?

[AlexS12]

You are right, right now the user has to select the right time-step.

I think, this instabilities that you showed us, are more related to stability regions of the RK than to Courant condition, but anyway, it's something we must deal with carefully.

Moreover, at some point we would like to try other schemes (ie. Adams-Bashforth) that I have seen they use in other simulators. It would be great to be able to do some numerical research on the integration of this and future systems.

[sanguinariojoe]

I'm afraid you are probably wrong, basically because you have a lot of explicit integrations in your model.

The first explicit scheme results from the decoupling of the Momentum and kinematic equations, at src/pyfme/models/euler_flat_earth.py, where you are separately integrating the speed, without considering the acceleration changes coming from momentum equation, in a explicit explicit (and therefore eventually unstable) way.

But that's not the very only place where you are actually using a explicit scheme. In fact, the forces are computed (sampled) just once by time step, regardless the internal changes of speed and position (as far as I know).

Hence, you can improve your internal integration scheme as many as you want... Right now, if you are not able to introduce a Courant condition your system would become unstable.

Anyway, since your integration is explicit, I also would choice an Adams-Basforth instead of a Runge-Kutta.

[AlexS12]

We still have some problems with numerical stability even if every object is updated at each internal time step. Maybe this is not well implemented or maybe we need to tune some of the internal parameters of the method.

Moreover, it seems that there is a problem with the conversion from z_earth to height that may be affecting the integration.

This is a working example:

from pyfme.aircrafts import Cessna172

from pyfme.environment.atmosphere import ISA1976
from pyfme.environment.wind import NoWind
from pyfme.environment.gravity import VerticalConstant
from pyfme.environment import Environment

from pyfme.utils.trimmer import steady_state_trim

from pyfme.models.state.position import EarthPosition

from pyfme.models import EulerFlatEarth

from pyfme.utils.input_generator import Constant, Doublet
from pyfme.simulator import Simulation


aircraft = Cessna172()
atmosphere = ISA1976()
gravity = VerticalConstant()
wind = NoWind()

environment = Environment(atmosphere, gravity, wind)

pos = EarthPosition(x=0, y=0, height=1000)
psi = 0.5  # rad
TAS = 45  # m/s
controls0 = {'delta_elevator': 0, 'delta_aileron': 0, 'delta_rudder': 0, 'delta_t': 0.5}

trimmed_state, trimmed_controls = steady_state_trim(
    aircraft,
    environment,
    pos,
    psi,
    TAS,
    controls0
)    

system = EulerFlatEarth(t0=0, full_state=trimmed_state)

de0 = trimmed_controls['delta_elevator']

controls = controls = {
    'delta_elevator': Doublet(t_init=2, T=1, A=0.1, offset=de0),
    'delta_aileron': Constant(trimmed_controls['delta_aileron']),
    'delta_rudder': Constant(trimmed_controls['delta_rudder']),
    'delta_t': Constant(trimmed_controls['delta_t'])
}

sim = Simulation(aircraft, system, environment, controls, dt=0.3)
results_03 = sim.propagate(25)

sim = Simulation(aircraft, system, environment, controls, dt=0.05)
results_005 = sim.propagate(25)

kwargs = {'subplots': True,
          'sharex': True,
          'figsize': (12, 100)}

ax = results_005.plot(marker='.', color='r', **kwargs)
ax = results_03.plot(ax=ax, marker='x', color='k', ls='', **kwargs)

plt.show()

These are the outputs:

[sanguinariojoe]

I certainly don't know nothing about the z_earth to height conversion, and its implications... Would it be related with the initial sock (~2.5 seconds)?

Anyway, these plots show the relevance of the time step. I thought a bit about that, and I think you'll pass a really hard time to get an actual Courant condition. However, you can simplify a bit the things featuring your airplane natural frequencies (by a decay test, maybe?), and then setting the time step as a fraction of the minimum resulting natural period.

Think about that and let's start discussing

[AlexS12]

The height z_earthmight be having influence in the atmosphere properties calculations and thus in the forces used to integrate. Anyway, it seem to be wrong for both simulations... I will have a deeper look at it.

I like your approach, @sanguinariojoe. I was thinking about setting the max_step option in the solver. As the order of natural frequencies for regular aircrafts is known, it should be easy to set a sensible time step for the integration, then the one chosen by the user would be related to the number of results that are saved.

On the other hand there is still something that intrigues me:
Assuming there aren't any errors on how different elements are updated (which probably is not true as this is only "the first attempt" **), the solver should choose an adequate adaptive time step and as every implied variable is updated: the full state vector, the forces and the moments (which wasn't the case before this refactor) the result of the integration should be independent of the time step chosen by the user. Do you think the same @Juanlu001?

---

** This is done in: pyfme/models/euler_flat_earth.py --> fun(self, t, x) (function to integrate), where self._update_full_system_state_from_state(x, self.state_vector_dot) and updated_simulation = self.update_simulation(t, self.full_state) are called.

What update_simulation does is:

    def update(self, time, state):
        self.environment.update(state)

        controls = self._get_current_controls(time)

        self.aircraft.calculate_forces_and_moments(
            state,
            self.environment,
            controls
        )
        return self

That besides that ugly return self, does in principle what is expected to do.

[astrojuanlu]

the solver should choose an adequate adaptive time step and as every implied variable is updated

The flow is a bit difficult to follow and I would need some time and pen and paper to figure out an answer to this :) Give me some days and I'll try to think something, if you haven't found the problem yet.

[AlexS12]

The height z_earthmight be having influence in the atmosphere properties calculations and thus in the forces used to integrate. Anyway, it seem to be wrong for both simulations... I will have a deeper look at it.

It was related to the position array dtype, solved in: 8dd95a8

Now we can focus on objects update and numerical stability