Solution for algebraic variables in DAE with callbacks

[baggepinnen]

The following example defines a DAE on mass-matrix form with callbacks. The last state p is the only algebraic state and the equation p-w(t) forces this state to be equal to the forcing function w which is a sinusoid. When solving this system, the solution looks strange unless dtmax is specified. Note, the solution for the differential states are correct in both cases, it's only the algebraic state solution that looks strange. If the callbacks are removed, the problem goes away. It looks like there are very few points in the trajectory for the algebraic states, and many of them do not satisfy the algebraic equation. Is this a problem with the error control for algebraic variables, do we perhaps need to specify additional tolerances for DAEs?

using OrdinaryDiffEq, LinearAlgebra, Plots, Unitful.DefaultSymbols, AutomaticDocstrings

function dynamics!(dx,x,θ,t)
    zmin, zmax, M, g = θ
    z, dz, p = x

    dx[1] = dz  
    dx[2] = -M*g -z + p
    dx[3] = p-w(t) # the last state (algebraic) p has to be equal to w

    if z <= zmin
        dx[2] = max(dx[2], 0.0)
    elseif z >= zmax
        dx[2] = min(dx[2], 0.0)
    end
end

w(t) = .1*(7.0+3.0sin(t))
zmin_cond(x,t,integrator) =  x[1]-zmin
zmax_cond(x,t,integrator) =  zmax-x[1]

function zmin_affect_neg!(integrator)
    integrator.u[1] = zmin
    integrator.u[2] = 0.0
end

function zmax_affect_neg!(integrator)
    integrator.u[1] = zmax
    integrator.u[2] = 0.0
end

sbs = CallbackSet(
    ContinuousCallback(zmin_cond,zmin_affect_neg!),
    ContinuousCallback(zmax_cond,zmax_affect_neg!)
)

tf = 20.0
tspan = (0.0,tf)

zmin = 1.0e-3
zmax = 2.0e-2
M = 0.05
g = 9.81
θ = zmin, zmax, M, g

x0 = [(zmin+zmax)/2,0.0,w(0.0)]

E = diagm([1.0,M,0.0])

f = ODEFunction(dynamics!,mass_matrix=E)
prob = ODEProblem(f,x0,tspan,θ)
sol1 = solve(prob,Rosenbrock23(),callback=cbs)
sol2 = solve(prob,Rosenbrock23(),callback=cbs, dtmax=0.1)
plot(sol1, layout=(3, 1), link=:x, lab="free dt")
plot!(sol2, lab="dtmax")

[YingboMa]

I think this is mostly an interpolation problem. The problem is that algebraic variables would be interpreted as having zero derivatives by the dense output.

[ChrisRackauckas]

See what happens with Rodas4 and Rodas5P.

[baggepinnen]

With Rodas4 the algebraic state looks much better, only some slight artifacts. With Rodas5P the entire solution including differential states is a disaster, shown below.

[ChrisRackauckas]

Rodas4, Rodas5, and Rodas5P are all fine now in https://github.com/SciML/OrdinaryDiffEq.jl/pull/1660 . Rosenbrock23 still is a bit weird, I'm trying to figure out what that could be because it doesn't look "wrong". The interpolation of algebraic variables seems fine without events, but only interpolation with algebraic variables is off after the event. This suggest something off with how add_steps! on Rosenbrock23 is doing the mass matrix, but I don't quite see it.

This model is silly though. Note that the event will never fire if it's computed at a high enough accuracy because it saturates asymptotically to the value, so the only way it can ever trigger is due to numerical error 😓. So it's a bad model, you want to trigger say -100eps() below to make it robust or something, otherwise the result as tolerances go to zero is different from what you think is correct, and what you think is correct is not correct.

[baggepinnen]

I'm not sure I follow you, the events do trigger at

zmin = 1.0e-3
zmax = 2.0e-2

and the model also has the term -M*g that drives negative without the event.

[ChrisRackauckas]

Keeping open to track for Rosenbrock23.