DifferentialEquations.jl

DifferentialEquations.jl copied to clipboard

MWE:

using DifferentialEquations

function usode!(du,u,p,t)
    C1 = 1.63e-9
    RS = 790.0
    C2 = 3.3e-9
    C0 = 0.3e-9
    L1 = 1e-3
    L2 = 82e-3
    CM = 4e-12
    LM = 8.0
    RM = 7500.0
    VC1, VC2, VC0, VCM, IL1, IL2, ILM, Vp = u
    fv, VAmag = p
    VA = VAmag*sinpi(2*Vp)
    du[1] = 1/C1*IL1
    du[2] = -1/C2*IL1
    du[3] = 1/C0*(IL1-IL2-ILM)
    du[4] = 1/CM*ILM
    du[5] = 1/L1*(VA-VC1-VC0-VC2-IL1*RS)
    du[6] = 1/L2*VC0
    du[7] = 1/LM*(VC0-VCM-ILM*RM)
    du[8] = fv
end

mutable struct Controller
    f::Float64
end
function(c::Controller)(integrator)
    integrator.p[1] = c.f
    if c.f < 30000.0
        c.f += 15.0
    end
    println(join([integrator.t, c.f, integrator.p[3], integrator.p[4]], ", "))
end

condVzero(u,t,integrator) = integrator.p[2]*sinpi(2*u[8])-u[1] # Vmod
function affectVzero!(integrator)
    integrator.p[3]=integrator.t
end
cbVz = ContinuousCallback(condVzero, affectVzero!, nothing)
condIzero(u,t,integrator) = u[5]
function affectIzero!(integrator)
    integrator.p[4]=integrator.t
end
cbIz = ContinuousCallback(condIzero, affectIzero!, nothing)

function sim()
    p = [0.0, 100.0, 0.0, 0.0]
    cb1 = PeriodicCallback(Controller(27000.0),0.005)
    cbs = CallbackSet(cb1,cbVz,cbIz)
    #cbs = CallbackSet(cb1,cbVz)
    #cbs = CallbackSet(cb1,cbIz)
    u0 = [0.0,0.0,0.0,0.0, 0.0,0.0,0.0, 0.0]
    tspan = [0.0, 0.5]
    prob = ODEProblem(usode!,u0,tspan,p)
    @time sol = solve(prob,Tsit5(), callback=cbs, reltol=1e-8, abstol=1e-8, dense=false, maxiters=10_000_000)
    sol
end

sol = sim()

Looking at function sim, there are three lines that begin with cbs = CallbackSet. If I run as is, it stalls after 0.2 s (simulation time in print statements) and eventually dies wanting larger maxiters… If I instead use either of the commented cbs = lines, each only having one of the two ContinuousCallbacks, then it succeeds.

Substituting Trapezoid() for Tsit5() gives an interesting error message rather than failing silently:

0.20500000000000002, 27630.0, 0.20498183392116934, 0.20498234571777418
┌ Warning: dt <= dtmin. Aborting. There is either an error in your model specification or the true solution is unstable.
└ @ DiffEqBase C:\Users\klaff\.julia\packages\DiffEqBase\XoVg5\src\integrator_interface.jl:343

The println output indicates that there were several very closely spaced callbacks just prior to the failure. I suspect the closely spaced callbacks are significant because I see a similar pattern in a non-NWE version.

Using VectorContinuousCallback instead of two separate ContinuousCallbacks as follows appears to work:

using DifferentialEquations

function usode!(du,u,p,t)
    C1 = 1.63e-9
    RS = 790.0
    C2 = 3.3e-9
    C0 = 0.3e-9
    L1 = 1e-3
    L2 = 82e-3
    CM = 4e-12
    LM = 8.0
    RM = 7500.0
    VC1, VC2, VC0, VCM, IL1, IL2, ILM, Vp = u
    fv, VAmag = p
    VA = VAmag*sinpi(2*Vp)
    du[1] = 1/C1*IL1
    du[2] = -1/C2*IL1
    du[3] = 1/C0*(IL1-IL2-ILM)
    du[4] = 1/CM*ILM
    du[5] = 1/L1*(VA-VC1-VC0-VC2-IL1*RS)
    du[6] = 1/L2*VC0
    du[7] = 1/LM*(VC0-VCM-ILM*RM)
    du[8] = fv
end

mutable struct Controller
    f::Float64
end
function(c::Controller)(integrator)
    integrator.p[1] = c.f
    if c.f < 30000.0
        c.f += 15.0
    end
    println(join([integrator.t, c.f, integrator.p[3], integrator.p[4]], ", "))
end

condVzero(u,t,integrator) = integrator.p[2]*sinpi(2*u[8])-u[1] # Vmod
function affectVzero!(integrator)
    integrator.p[3]=integrator.t
end
cbVz = ContinuousCallback(condVzero, affectVzero!, nothing)
condIzero(u,t,integrator) = u[5]
function affectIzero!(integrator)
    integrator.p[4]=integrator.t
end
cbIz = ContinuousCallback(condIzero, affectIzero!, nothing)

function condition(out,u,t,integrator)
    out[1] =  integrator.p[2]*sinpi(2*u[8])-u[1]
    out[2] =  u[5]
end
function affect!(integrator, idx)
    if idx == 1
        integrator.p[3]=integrator.t
    elseif idx ==2
        integrator.p[4]=integrator.t
    end
end
vccb = VectorContinuousCallback(condition,affect!,2)

function sim()
    p = [0.0, 100.0, 0.0, 0.0]
    cb1 = PeriodicCallback(Controller(27000.0),0.005)
    #cbs = CallbackSet(cb1,cbVz,cbIz)
    cbs = CallbackSet(cb1,vccb)
    #cbs = CallbackSet(cb1,cbVz)
    #cbs = CallbackSet(cb1,cbIz)
    u0 = [0.0,0.0,0.0,0.0, 0.0,0.0,0.0, 0.0]
    tspan = [0.0, 0.5]
    prob = ODEProblem(usode!,u0,tspan,p)
    @time sol = solve(prob,Tsit5(), callback=cbs, reltol=1e-8, abstol=1e-8, dense=false, maxiters=10_000_000)
    sol
end

sol = sim()

Nice that it worked for you, but still we need to fix that, I am busy till Tuesday, I will get back to the issue once I am free.

Nice that it worked for you, but still we need to fix that, I am busy till Tuesday, I will get back to the issue once I am free.

I try to leave breadcrumbs in case someone else has the same issue and also let you know I'm not dead in the water so you can prioritize appropriately.

As a note of appreciation - this is an amazing library - I'm beating on it pretty hard and in the process showing coworkers what is possible. Thanks!