ChrisRackauckas
commented
4 years ago Is deciding whether to update the Jacobian up to Sundials.jl or Sundials IDA? I found from online at http://www.scholarpedia.org/article/Sundials_equation_solvers#IDA.2C_for_DAE_Systems, saying that The nonlinear iteration is Modified or Inexact Newton, and Jacobian information is updated only when necessary for convergence, not at every step.
It's up to Sundials IDA. However, you can probably trigger it by doing a reinit!. That should happen naturally if you add a discrete callback at the time and then just do u_modified!(integrator,true). Then the reinit! should require a new Jacobian to be drawn. A reinit is required anyways since the history values wouldn't be correct and it would need to hit back to first order.
cuihantao
jd-lara
TL;DR: IDA won't converge after an algebraic equation switching, due to the fact that Jacobian is only updated once for each step.
Hi,
First of all, thanks for developing and supporting the DifferentialEquations.jl eco-system. I found this work impressive and have started trying to use it for power system simulation problems, which are formulated as index-1 DAEs with a large set of nonlinear algebraic equations. At predefined times, algebraic equations will be modified, such as removing part of the equation to simulate an out-of-service of a device.
My simulation setup is a hybrid of Python and Julia for now. The Python part implements the models, i.e., the functions that compute the residuals of
G(u, du, t) = resand the Jacobians given bygamma * dG/d(du) + dG/du. The changes of algebraic equations at predefined times are also implemented in Python because for now they are part of the model.The Julia part constructs the DAEFunction, DAEProblem and uses IDA for solving. The Julia part looks like
The switching event happens at the end of
t=2.0.switch_times = [1.9999, 2.0, 2.0001]to force the solver to stop right before, at and after the switching. I have tried to switch out difference devices and simulate the DAE problem. Results of the solved ones can match that from my Python solver implementing the Trapezoidal rule. But there are cases that IDA won't converge after the switching.I tried to print out the maximum mismatch (residual) from Julia and compare with my Python code. The debugging printouts and my interpretations are as follows.
h=0.000025and updates the Jacobian once. Since the algebraic switching is drastic, the approximated Jacobian failed to make it converge ath=2.000025. Notice the residuals are6.784984056352011(iter. 1),0.3345103040021776(iter. 2), and0.17705125103928232(iter. 3).h=0.00000625and updates the Jacobian once. No convergence either. Notice the residuals are6.784984056352011(iter. 1),0.33448411910463083(iter. 2), and0.17703102372314716(iter. 3). The residuals are almost the same as that withh=h=2.000025(with the tiny differences caused by the differential variables).I compared the residuals with my Python program. Residuals from iter. 1 and iter. 2 match identically with my Trapezoidal rule solver. Starting from iter. 3, IDA's residuals are much larger than that from Trapezoidal solver. Since the Jacobian only accurate to for iter. 1 and hardly usable for the subsequent iterations, it explains why IDA won't converge.
Questions:
Is deciding whether to update the Jacobian up to
Sundials.jlor Sundials IDA? I found from online at http://www.scholarpedia.org/article/Sundials_equation_solvers#IDA.2C_for_DAE_Systems, saying thatThe nonlinear iteration is Modified or Inexact Newton, and Jacobian information is updated only when necessary for convergence, not at every step.Is there any workaround? It seems that I cannot force a call to the
jac!at each iteration, unless I use some messy global variables to store theJmatrix and forcefully write to it in the residual function.Ideally, an additional logic for updating Jacobian can be implemented - either one Jacobian update per step, or one update per iteration. In the proximity of
tstopsor after a non-convergence, update the Jacobian per iteration.Thanks in advance for your time.