Infinite loop in PID stepsize controller due to float precision error #703
Description
There is a special case that leads to the same step in PID being recalculated infinitely.
This arises in PIDController.adapt_step_size() when t0 is high, dt is clipped to dtmin, and force_dtmin=True
dt sets next_t1 as: next_t1 = next_t0 + dt
this is then read back next loop and prev_dt is recalculated as: prev_dt = t1 - t0
Due to float precision error, sometimes prev_dt != dt, and (prev_dt <= self.dtmin) ends up False
So, keep_step stays False, but dt never drops below dtmin.
See below example:
import jax
import jax.numpy as jnp
jax.config.update('jax_enable_x64', True)
DTYPE = jnp.float64
# Case:
# dt is clipped to dtmin, and we wish to continue solver (force_dtmin=True)
# let's say dtmin = 0.2
dt = jnp.asarray(0.2, dtype=DTYPE)
# Let's say we are at an arbitrary time point>>dtmin in the solver
next_t0 = jnp.asarray(967.1378944377053, dtype=DTYPE)
next_t1 = next_t0 + dt
# next_t0, next_t1 are returned
# an ODE step is taken
# adapt_step_size() is called again with t0=next_t0, t1=next_t1
# We calculate the actual dt taken
# prev_dt = t1 - t0
prev_dt = next_t1 - next_t0
print(f"desired dt: {dt}; actual dt: {prev_dt}")
# Because t0 >> dt, there is a chance for float error such that prev_dt > dt (and dtmin)
# in keep_step = keep_step | (prev_dt <= self.dtmin)
print(prev_dt <= dt) # is FALSE!
# then, keep_step=False
# factor < 1
# dt is clipped to dtmin again
# ...
# Infinite loop!My solution was to add a flag in the controller_state (called keep_next_step)
if self.dtmin is not None:
keep_step = keep_step | keep_next_stepif self.dtmin is not None:
if not self.force_dtmin:
result = RESULTS.where(dt < self.dtmin, RESULTS.dt_min_reached, result)
# flag next step to be kept if dtmin is reached
# or if it was reached previously and dt is unchanged
keep_next_step = (dt <= self.dtmin) | (keep_next_step & (factor == 1))
dt = jnp.maximum(dt, self.dtmin)Which works for my stuff, but unsure if it breaks things down the line..