Open PavelBal opened 2 years ago
Great initiative I would suggest to put in in the docs once finalized
Hi @PavelBal , Regarding the CVODES/CVODE, they are multiple solvers. Two often used ones are bdf (for stiff) and adams (for non-stiff). Stiff case is hard to integrate than non-stiff case. But bdf takes a bit longer integration time than adams. In practice, for a new set of ODEs, I would start with bdf solver first. If the simulation looks fine, check the simulation from adams. Somehow, the choice of solvers would depend on the ultimate goal of the ODE-model: if it is only for one-time simulation. I prefer to choose bdf. If it is related to optimization/inference, some accuracy-efficiency trade-off needs to consider. I believe that can vary in a case-by-case manner and depends on available approaches for optimization/inference.
Besides the concern of numerical solver, I believe the failure of 'ideal' solver can result from purely mathematical expression of the right-hand side (RHS). Please check this. Indeed wrong choice of a RHS can cause later numerically unstable simulation, which is independent of the choice of numerical solvers.
Hi All,
I think it would be good to have some more details on this in the documentation and some suggestions for optimization. Especially if you want to speed up a batch run of simulations, this can also make a big difference for computation time in my experience. I found this explanation helpful. It also points to the importance of the relative and absolute tolerance ratio, as this decides at which point the absolute tolerance becomes dominant. In general, I think reducing the relative tolerance is mainly relevant for models which are highly nonlinear and contain rapidly changing values, e.g., target binding at a high and saturating dose. Reducing the absolute tolerance mainly makes sense if you simulate values below 1e-10 (in base units I assume) which are still relevant and need to be precise, or if you reduce the relative tolerance. Reducing Hmax might be mainly relevant for problems that include tipping points rather than a gradual approach to steady-state. If there is a long phase (relative to 60 minutes) with rather linear behavior, increasing Hmax can speed up the simulations significantly without losing precision.
@msevestre @Yuri05 @Huan-Yang @StephanSchaller @abdelr @VanessaBa
I am not sure where this discussion belongs to, so I start it here and we can transfer it later.
We should have some kind of a wiki for best practices how to parametrize the solver settings available in PK-Sim/MoBi. Furthermore, we can collect ideas which improvements we can implement, and how to make models more robust. Feel free to edit this post.
Solver settings
AbsTol
Default value is
1e-10
What it does - ???
General rule - decreasing the value improves the stability but increases computational time
Experience: for many "simple" models, values between 1e-7 and 1e-9 are sufficient
For complex models (e.g. QSP, DDI), decreasing the value up to 1e-12 may be required
Values lower than 1e-12 often lead to less stability
RelTol
Default value is 1e-5
What it does - ???
General rule - decreasing the value improves the stability but increases computational time
Experience: I sometimes decrease this value to 1e-7, but I do not really have a rule of thumb what to decrease -
AbsTol
orRelTol
I heard that
AbsTol
andRelTol
should not have values too far apart - but I do not really know what "too far apart" means and how to set these valuesH0
Default value is 1e-10
Initial step size (in minutes?)
Experience: Never changed it
HMin
Default value is 0
Minimal allowed step size
Experience: Never changed it
HMax
Default value is 60
Maximal allowed step size (minutes)
Reducing this value might improve the stability of the model but increases the computational time. The value should be selected with regard to the expected kinetics of the model. E.g. for the diabetes model, a value of 5 minutes makes more sense, as we do not expect any simulation scenarios where the RHS remain unchanged for 60 minutes.
MxStep
Default value is 100000
The default value is probably too high and leads sometimes to long computational times before the solver actually breaks and returns an error.
UseJacobian
Default is TRUE
Modeling
Events
TIME
as condition, use>=
resp<=
istead of=
(or was it the other way around?)Time > 1000 ? X : Y
. Create an event atTime = 1000
and chage the value of the parameter fromX
toY
.