Default steady-state-by-simulation tolerances should be higher

[dilpath]

Can include a MWE if desired, not yet prepared.

Using the Blasi problem [1] with a parameter vector supplied by @plakrisenko .

At default tolerances, using simulation to find the steady state, wrms_ [2] seems to fluctuate a bit and never gets below the threshold of 1. Printed here are the wrms_ values at [2] that were printed immediately before the error.

152.175
26.3248
152.175
26.3248
152.175
99.5254
152.175
99.5254
152.175
99.5254
152.175
26.3249
26.3249
152.175
26.3249
26.3249
26.3249
152.175
26.3249
26.3248
26.3247
152.175
26.3248
26.3249
152.175
26.3248
26.3249
26.3249
26.3248
26.3248
26.3248
26.3248
152.175
26.3249
152.175
26.3248
152.175
99.5253

[Warning] AMICI:CVODES:CVode:OTHER: AMICI ERROR: in module CVODES in function CVode : At t = 5.01122e+199, the right-hand side routine failed in an unrecoverable manner. 
[Warning] AMICI:equilibration: AMICI equilibration failed: AMICI failed to integrate the forward problem

[Warning] AMICI:simulation: AMICI simulation failed:
Steady state computation failed. First run of Newton solver failed: No convergence was achieved. Simulation to steady state failed. Second run of Newton solver failed: No convergence was achieved.
Error occurred in:
0       0x7fe37c02f3f6 amici::AmiException::AmiException() + 38
1       0x7fe37c02f4bb amici::AmiException::AmiException(char const*, ...) + 123
2       0x7fe37bf862e0 /home/dilan/Documents/future_annex/steady_state_paper/20220331_debug_conservation_laws/before/packages/AMICI/python/sdist/amici/_amici.cpython-38-x86_64-linux-gnu.so(+0xab2e0) [0x7fe37bf862e0]
3       0x7fe37c0862c9 /home/dilan/Documents/future_annex/steady_state_paper/20220331_debug_conservation_laws/before/packages/AMICI/python

Printed here is the same thing, but with the steady state tolerances set to two orders of magnitude above default simulation tolerances. The wrms_ values decrease monotonically (unlike above) to a value less than the threshold 1,.

1291.25
1049.53
853.222
693.769
564.23
458.968
373.41
303.853
247.294
201.298
163.883
133.443
108.672
78.884
57.2981
41.6431
30.2707
22.0027
15.9956
9.63246
5.86245
3.6014
2.20007
1.32276
# Success
0.794022

A possible solution would be to change steady-state-by-simulation tolerances to be some factor of simulation tolerances, e.g.:

solver->setSteadyStateToleranceFactor(1e+2)

// Replace everywhere `solver->getAbsoluteToleranceSteadyState()` with simulation tolerance increased by the factor:
solver->getAbsoluteTolerance() * solver->getSteadyStateToleranceFactor()

// Replace everywhere `solver->getRelativeToleranceSteadyState()` with simulation tolerance increased by the factor:
solver->getRelativeTolerance() * solver->getSteadyStateToleranceFactor()

[1] https://github.com/Benchmarking-Initiative/Benchmark-Models-PEtab/tree/master/Benchmark-Models/Blasi_CellSystems2016
[2] https://github.com/AMICI-dev/AMICI/blob/master/src/steadystateproblem.cpp#L456=

[FFroehlich]

Sounds reasonable, that's also how I have been setting those tolerances. Recent observations suggest that this problem is also partially resolved by using the newton step as convergence criteria, as implemented in #1737.

I am a bit hesitant to completely remove routines to set absolute and relative steadystate tolerances. I would prefer if the values set by those routines still take precedence over the ratio set by solver->setSteadyStateToleranceFactor(1e+2).

[dweindl]

Thanks for looking into that. I agree that changing default tolerances would be a good idea. 10x or 100x the integration tolerances. I would still prefer handling it via absolute abs/rel tolerances.