Problems with abstol / abstolv of integrator

[Joost Lauwers]

Hi,

I'm working on a parameter estimation problem that involves a fairly large stiff DAE system in which some of the states are small ( in the order 1e-8). I'm having several problems and I would like to point out something very peculiar imo.

In my parameter estimation script I get errors of the following kind:

 " Newton/Linesearch algorithm failed to converge.eval_h failed:  on line 1155 of file ".../casadi/casadi/interfaces/sundials/idas_internal.cpp. Module "IDACalcICB" returned flag -4 ("IDA_CONV_FAIL"). Consult Idas documentation." 

I've played around a bit with the number of parameters I estimate, the objective function and the tolerance errors and I believe the failure is, at least partially, due to the stiffness of the problem. In order to have some more control on this, I would like to set the error tolerance specifically per state I am integrating:

Thus, for instance by setting 

>>> abs_tol = [1e-4]*2 + [1e-5]*2

>>> integrator.setOption("abstolv",abs_tol)

instead of 

>>> integrator.setOption("abstol",1e-4)

Now, I believe that the first command doesn't really does what I should do. I've noticed that  for instance setting

>>> abs_tol = [1e-4]*4 # 4 is the number of states

>>> integrator.setOption("abstolv",abs_tol)

or

>>> integrator.setOption("abstol",1e-4)

does NOT give the same result (while I would expect otherwise). That is: in the first case the solver crashes with the error "Module "IDACalcIC" returned flag -22 ("IDA_ILL_INPUT"). Consult Idas documentation.".

My code is a bit lenghty and complex to put here so I made a 'simple' artificial parameter estimation problem script (see attachment "Ridiculous example2_estimation_single_shooting_forum") in which the parameters of a stiff DAE system are estimated. I hope the structure of the code is a bit clear: I generate data using the analytical solution of a system of DAEs + some noise. Based on this data I estimate the parameters, repeat the integration with the found parameters and plot the results. At line 73-80 you'll see the following

>>> abs_tol = [1e-13]*4

>>> integrator.setOption("abstolv",abs_tol) # tolerance

>>> #abs_tol = 1e-13

>>> #integrator.setOption("abstol",abs_tol) # tolerance

If you use the first tolerance setting, I cannot calculate the parameters. Using the second set, the parameter estimation works.

So, is there a problem with integrator.setOption("abstolv",...) , why does the first not work and how can I set the error for each state independently? 

Curiously, this problem does NOT occur when I use my parameter esimation script as a simulator, as in the file of attachment " Ridiculous example2_estimation_simulation_single_shooting_forum". This is the same problem, but I fixed the parameters and the initial conditions on the right values, I removed the constraint that the initial conditions have to be consistent with the algebraic equations and I optimize my objective function f = 0.0 for a parameter q which does not appear in the equations. (this is just a simulation of course). Here both settings, i.e.,

>>> abs_tol = [1e-13]*4

>>> integrator.setOption("abstolv",abs_tol) # tolerance

or 

>>> #abs_tol = 1e-13

>>> #integrator.setOption("abstol",abs_tol) # tolerance

 lead to a succesful simulation. I really do not understand what is different here than in the normal estimation procedure of 'Ridiculous_example2_estimation_single_shooting_forum.py'.

I hope there's some explanation so I can work on it in my original problem!

Kind regards and many thanks!

Joost

[Joel Andersson]

Hi Joost!

A failure in IDACalcIC or IDACalcICB suggests that IDAS cannot find consistent initial conditions. This can be resolved by giving a better initial guess for the initial conditions (input "z0" for IDACalcIC and "rz0" for IDACalcICB). Do you manually set "z0" and/or "rz0" or do you rely on the default values?

Note that in CasADi, it's easy to use a tool such as IPOPT to get the consistent "z0".

Joel

[Joost Lauwers]

I forgot about the rz0 !  I've implemented it and at least the solver seems to work now (although I seem to get stuck in an infeasible solution, but maybe I've made an other error somewhere)

Still, I think it's curious that 

>>> abs_tol = [1e-4]*4 # 4 is the number of states

>>> integrator.setOption("abstolv",abs_tol)

and

>>> integrator.setOption("abstol",1e-4)

do not seem to give the same results in my 'example files'

In any case. Thanks for the help!

Joost

[Joel Andersson]

Hi Joost,

CasADi just passes the information to SUNDIALS/IDAS. The length of the vector should be the number of states + number of algebraic variables. Maybe this is the problem?

Otherwise, check the SUNDIALS/IDAS user guide for IDASVtolerances and IDASStolerances respectively.

Best regards,
Joel