"[IDAS ERROR] IDAGetDky Illegal value for k." with event callbacks

[JuhaHeiskala]

When solving the "Robertson model" with event callbacks the following error occurs:

[IDAS ERROR] IDAGetDky
Illegal value for k.

┌ Warning: IDAGetDky failed with error code =
│ flag = -25
└ @ Sundials ~/.julia/dev/Sundials/src/simple.jl:20

The event defitions are from Octave ode15i function tests. The events should occur nearly simultaneously, but with the error the second event is not detected.

The error happens when "save_everystep=true" option is given to "solve".
If "save_everystep=false" the error does not happen and the both events are correctly detected.
Also without events the problem is solved successfully.

Julia 1.5.3 and DifferentialEquations 6.15 with Sundials v4.3.0

Below is test code.

import DifferentialEquations  
const DiffEq = DifferentialEquations

function rob(out, du, u, p, t)
    out[:] .= [-(du[1] + 0.04*u[1] - 1e4*u[2]*u[3]),
          -(du[2] - 0.04*u[1] + 1e4*u[2]*u[3] + 3e7*u[2]^2),
          u[1] + u[2] + u[3] - 1]
    
    return out
end

function cond1(u, t, integrator)
    if (t < 1e1)
        val = -1
    else
        val = 1
    end
    return val    
end

function cond2(u, t, integrator)
    if (t < 1e1)
        val = -2
    else
        val = 3
    end
    return val    
end


dae_f = DiffEq.DAEFunction(rob)

tspan = (0.0, 100.0)
u0 = [1.0, 0.0, 0.0]
du0 = [-1e-4, 1e-4, 0.0]
differential_vars = [true, true, false]

dae_prob = DiffEq.DAEProblem(dae_f, du0, u0, tspan; differential_vars)

affect1 = (integrator)-> (println("Event 1 at time: $(integrator.t)"); integrator)
affect2 = (integrator)-> (println("Event 2 at time: $(integrator.t)"); DiffEq.terminate!(integrator))

cb1 = DiffEq.ContinuousCallback(cond1, affect1, nothing)
cb2 = DiffEq.ContinuousCallback(cond2, affect2)

cb_set = DiffEq.CallbackSet(cb1, cb2)

println("Sol 1")
sol1 = DiffEq.solve(dae_prob, DiffEq.IDA(); save_everystep=false,callback=cb_set)

println("\nSol 2")
sol2 = DiffEq.solve(dae_prob, DiffEq.IDA(); save_everystep=true,callback=cb_set)

println("\nSol 3")
sol1 = DiffEq.solve(dae_prob, DiffEq.IDA(); save_everystep=true)

[JuhaHeiskala]

The issue does not happen with Sundials 3.8.3.
The issue does happen with Sundials 3.9

[ChrisRackauckas]

function cond1(u, t, integrator)
    if (t < 1e1)
        val = -1
    else
        val = 1
    end
    return val    
end

This doesn't make sense. The condition of a ContinuousCallback is supposed to be a rootfinding function, and this won't have a root. First of all, this should probably be a DiscreteCallback? But secondly, the condition as a rootfinding problem would be t-1e1

[JuhaHeiskala]

I wondered the same. Like I said the first post in Discourse, this problem is from Octave ode15i-function tests when Sundials is available in Octave. I am guessing that the event definitions are intended to be a "tough" case for the event finding and not a realistic problem.

[ChrisRackauckas]

In our world, it's not a tough case but an ill-defined one. @kanav99 it would be interesting to somehow support non-root cases though.

[JuhaHeiskala]

My guess is that the purpose of the discontinuous event function is to mimic some possible real world case where there would be need to differentiate two almost simultaneous events.

Anyway, the discontinuity does not seem to be the issue.
I get the same IDA error with conditions defined as:

function cond1(u, t, integrator)
    tanh(t-9.0)
end

function cond2(u, t, integrator)
    tanh(t-10.0)
end

and actually even with:

function cond1(u, t, integrator)
    t-9
end

function cond2(u, t, integrator)
    t-10
end

Both the above cases work correctly with Sundials 3.8.3

[JKRT]

Saw your latest post. t-9 should be a zero-cross right?
Might be similar to the issue I am having here
Works with Sundials from C, not via Sundials.jl might be the same root cause:)

[JuhaHeiskala]

Yes, it should detect two events at t=9 and t=10. The first t=9 event is detected, the second is not, instead the IDA error is generated.
Have you tried your example with Sundials 3.8.3? If it does work, then maybe that indicates the same root cause. (though I have no knowledge of Sundials, so just guessing.)

[ChrisRackauckas]

SciML uses its own rootfinding because it fixes a few root misdetection issues. However, it always looks for the rootfinding 0. We haven't special-cased for the chance that the sign change is at a point which is not a zero, but a shift from -Inf to Inf. We can do that though.

[ChrisRackauckas]

Oh just read the other message. That's... weird. And it's an interpolation issue that Sundials throws? That needs to be looked into. Maybe it's something to do with the reinitialization

[JuhaHeiskala]

Below is the error:
Sol 2
Event 1 at time: 8.999999999999998

[IDAS ERROR] IDAGetDky
Illegal value for k.

┌ Warning: IDAGetDky failed with error code =
│ flag = -25
└ @ Sundials ~/.julia/packages/Sundials/xAoTr/src/simple.jl:20

[copybit]

I'm experiencing the same, or a similar, issue. I want to simulate a hybrid DAE system with a zero-order-hold controller. I use ContinuousCallbacks for triggering jumps in the dynamics and a PeriodicCallback to run the controller.

The attached example is a toy example. I realize that it is an ODE I am solving, and it does work fine to formulate as an ODEproblem and solve with Tsit5, in which case I have no issues with CallbackSets that combine ContinuousCallbaks and PeriodicCallbacks. However, I eventually want to simulate a DAE that is not an ODE.

using DifferentialEquations
using Sundials

# ODE disquised as DAE
function dyns!(y,dx,x,p,t)
    y[1] = dx[1] - x[1]
end

# xmin callback
xmin= -0.5
function xmin_cond(x,t,integrator)
   x[1]-xmin
end
function xmin_affect_neg!(integrator)
    integrator.u[1] = 0.01    # Reset x1. A bit confusing: integrator.u is x, DifferentialEquations reserves u for the state
end
xmin_cb = ContinuousCallback(xmin_cond,nothing,xmin_affect_neg!)

# xmax callback
xmax= 1.0
function xmax_cond(x,t,integrator)
   -(x[1]-xmax)
end
function xmax_affect_neg!(integrator)
    integrator.u[1] = -0.01    # Reset x1. A bit confusing: integrator.u is x, DifferentialEquations reserves u for the state
end
xmax_cb = ContinuousCallback(xmax_cond,nothing,xmax_affect_neg!)

# Nonsense ZOH controller with state in struct
mutable struct Controller_state
    u
    Controller_state() = new()
end
cs = Controller_state()
cs.u = 0.01

h = 1.0
function control_affect!(integrator)
    #nothing
    integrator.u[1] = 0.01
end

control_cb = PeriodicCallback(control_affect!,h)

# Initial state
x0 = [0.5]
dx0 = 0*x0

# Simulation time
tf=10.0
tspan = (0.0,tf)

# Formulate, simulate, plot
differential_vars = [true]
prob = DAEProblem(dyns!,dx0,x0,tspan,differential_vars=differential_vars)

#cbs=CallbackSet(control_cb) # this works
#cbs=CallbackSet(xmin_cb, xmax_cb) # this works
#cbs=CallbackSet(xmin_cb,control_cb) # this works
cbs=CallbackSet(xmin_cb,xmax_cb,control_cb) # this breaks - even if I set control_affect!(integrator) = nothing

sol = solve(prob,IDA(),callback=cbs)
plot(sol,vars=(1))

[ChrisRackauckas]

I see, this issue only happens when you have a callback set and do two different reinits in the same step.