callbacks_buckExample

[mongibellili]

Greetings,
I tried to simulate the electronic buck converter using the following code:

using DifferentialEquations
using Plots
function odeDiffEquPackage()
    function f(du, u, p, t)
        C = 1e-4; L = 1e-4; R = 10;U = 24.0; T = 1e-4; DC = 0.5; ROn = 1e-5;ROff = 1e5;
        rd=p[1];rs=p[2];
        il=u[1] ;uc=u[2]
        id=(il*rs-U)/(rd+rs) # diode's current
        du[1] =(-id*rd-uc)/L # inductor's current
        du[2]=(il-uc/R)/C    # capacitor's voltage
    end
    function condition1( u, t, integrator) 
        (t-p[3])
    end
    function condition2( u, t, integrator) 
        (t-p[4]-0.5*1e-4)
    end
    function condition3( u, t, integrator) 
         (p[5]*((u[1]*p[2]-24.0)/(p[1]+p[2]))+(1.0-p[5])*((u[1]*p[2]-24.0)*p[1]/(p[1]+p[2])))
    end
    function affect1!(integrator)
                p[4]=p[3]
                p[3]=p[3]+1e-4
                p[2]=1e-5
    end
    function affect2!(integrator)
                p[2]=1e5
    end
    function affect3!(integrator)
            p[1]=1e-5
            p[5]=1.0
    end
    function affect33!(integrator)
        p[1]=1e5
        p[5]=0.0
    end
    cb1 = ContinuousCallback(condition1, affect1!,nothing;  )
    cb2 = ContinuousCallback(condition2, affect2!,nothing; )
    cb3 = ContinuousCallback(condition3, affect3!,affect33!;  )
    cbs = CallbackSet(cb1, cb2,cb3)
    u0 = [0.0, 0.0]
    tspan = (0.0, 0.001)
    p = [1e5,1e-5,1e-4,0.0,0.0,0.0]
    prob = ODEProblem(f, u0, tspan, p)
    sol = solve(prob, Tsit5(), callback = cbs, reltol=1e-3,abstol=1e-4#= dt = 1e-3, adaptive = false =#)
    p1=plot!(sol);
    savefig(p1, "Tsit5()_34_buck_ft0005_")
 end
 odeDiffEquPackage() 

These are the results:

However, after simulating the system via 3 other tools:

• ltspice: https://www.analog.com/en/design-center/design-tools-and-calculators/ltspice-simulator.html
• c solver https://github.com/CIFASIS/qss-solver
• julia qssolver https://github.com/mongibellili/QS_Solver using the same system

using qss
function test()
    odeprob = @NLodeProblem begin
          name=(buck,)
          C = 1e-4; L = 1e-4; R = 10.0;U = 24.0; T = 1e-4; DC = 0.5; ROn = 1e-5;ROff = 1e5;
          discrete = [1e5,1e-5,1e-4,0.0,0.0]
          u = [0.0,0.0]
          rd=discrete[1];rs=discrete[2];nextT=discrete[3];lastT=discrete[4];diodeon=discrete[5]
          il=u[1] ;uc=u[2]
          id=(il*rs-U)/(rd+rs) # diode's current
          du[1] =(-id*rd-uc)/L
          du[2]=(il-uc/R)/C
          if t-nextT>0.0 
            lastT=nextT
            nextT=nextT+T
            rs=ROn
          end
          if t-lastT-DC*T>0.0 
            rs=ROff
          end                          
          if diodeon*(id)+(1.0-diodeon)*(id*rd)>0
          rd=ROn
          diodeon=1.0
          else
          rd=ROff
          diodeon=0.0
          end     
    end
    tspan = (0.0, 0.0002)
    sol= solve(odeprob,nmliqss2(),tspan,abstol=1e-4,reltol=1e-3)    
    save_Sol(sol)
end
test()

I received the following results:
c solver

ltspice

QS_Solver.jl

Am I misusing callbacks?
thank you.

[oscardssmith]

what happens if you simulate this without adaptive false? turning off adaptivity could easily lead to loss of precision

[mongibellili]

I tried both (with and without adaptive), and tried to change dt. still no correct output.

[pepijndevos]

From the graph it seems like Rd isn't working correctly. Could you add some printing to verify Rd has the value you expect and condition3 is producing the expected result and affect3 and affect33 are running when that happens?

[mongibellili]

From the graph it seems like Rd isn't working correctly. Could you add some printing to verify Rd has the value you expect and condition3 is producing the expected result and affect3 and affect33 are running when that happens?

I am printing condition3 and rd from inside affect3 and affect33.

 function affect3!(integrator)
            zcfPOS= (p[5]*((integrator.u[1]*p[2]-24.0)/(p[1]+p[2]))+(1.0-p[5])*((integrator.u[1]*p[2]-24.0)*p[1]/(p[1]+p[2])))
            @show zcfPOS,p[1]
            p[1]=1e-5
            p[5]=1.0
    end
    function affect33!(integrator)
        zcfNEG= (p[5]*((integrator.u[1]*p[2]-24.0)/(p[1]+p[2]))+(1.0-p[5])*((integrator.u[1]*p[2]-24.0)*p[1]/(p[1]+p[2])))
        @show zcfNEG,p[1]
        p[1]=1e5
        p[5]=0.0
    end

The output:
(zcfNEG, p[1]) = (3.254285729781259e-12, 100000.0)
(zcfPOS, p[1]) = (-2.225775119768514e-12, 100000.0)
(zcfNEG, p[1]) = (6.901359483200565e-15, 1.0e-5)
(zcfNEG, p[1]) = (9.991474314574589e-11, 100000.0)
(zcfNEG, p[1]) = (1.872813015779684e-10, 100000.0)
(zcfNEG, p[1]) = (5.0528470296740124e-11, 100000.0)
(zcfNEG, p[1]) = (3.581703822419513e-10, 100000.0)
(zcfNEG, p[1]) = (4.676419251836705e-10, 100000.0)
(zcfNEG, p[1]) = (5.983977757750836e-10, 100000.0)
(zcfNEG, p[1]) = (6.473275249163635e-10, 100000.0)

So you are right. it looks like the issue is in condition3: the zero crossing function is positive inside affect33. when I flip the sign of condition3 or switch positions of affect3 and affect33, the code runs infinitely.

[pepijndevos]

What do you mean "runs infinitely"?

[mongibellili]

Well, I waited 5 minutes for it and it did not finish. with a verbose mode, it was printing continuously. ie condition3 was being triggered a lot.

[pepijndevos]

That's probably because you change the condition equation in the affect function, so it ends up oscillating unable to find the zero point. Why isn't it just using id?

[contradict]

I suspect Tsit5 is the wrong solver for this problem. Using the automatic algorithm selection by changing your solver call to

sol = solve(prob; alg_hints=[:stiff],  callback = cbs)

Results in a plot that looks more like the other solvers, but it appears that the inductor current is in mA rather than A.

[pepijndevos]

The following code seems to miss the positive callback entirely. Surely that's a bug right? What goes down must come up.

using DifferentialEquations
using OrdinaryDiffEq
using Plots
function f(du, u, p, t)
    C = 1e-4; L = 1e-4; R = 10;U = 24.0; T = 1e-4; DC = 0.5; ROn = 1e-5;ROff = 1e5;
    rd=p[1];rs=p[2];
    il=u[1] ;uc=u[2]
    id=(il*rs-U)/(rd+rs) # diode's current
    du[1] =(-id*rd-uc)/L # inductor's current
    du[2]=(il-uc/R)/C    # capacitor's voltage
end
function condition1( u, t, integrator) 
    (t-p[3])
end
function condition2( u, t, integrator) 
    (t-p[4]-0.5*1e-4)
end
function condition3( u, t, integrator) 
    C = 1e-4; L = 1e-4; R = 10;U = 24.0; T = 1e-4; DC = 0.5; ROn = 1e-5;ROff = 1e5;
    rd=p[1];rs=p[2];
    il=u[1] ;uc=u[2]
    id=(il*rs-U)/(rd+rs) # diode's current
    diodeon=p[5]
    res=diodeon*(id)+(1.0-diodeon)*(id*rd)
end
function affect1!(integrator)
            p[4]=p[3]
            p[3]=p[3]+1e-4
            p[2]=1e-5
end
function affect2!(integrator)
            p[2]=1e5
end
function affect3!(integrator)
        @show "pos", p[1], condition3(integrator.u, integrator.t, integrator)
        p[1]=1e-5
        p[5]=1.0
        @show "pos", p[1], condition3(integrator.u, integrator.t, integrator)
end
function affect33!(integrator)
    @show "neg", p[1], condition3(integrator.u, integrator.t, integrator)
    p[1]=1e5
    p[5]=0.0
    @show "neg", p[1], condition3(integrator.u, integrator.t, integrator)
end
cb1 = ContinuousCallback(condition1, affect1!,nothing;  )
cb2 = ContinuousCallback(condition2, affect2!,nothing; )
cb3 = ContinuousCallback(condition3, affect3!,affect33!;  )
cbs = CallbackSet(cb1, cb2,cb3)
u0 = [0.0, 0.0]
tspan = (0.0, 0.001)
p = [1e5,1e-5,1e-4,0.0,0.0,0.0]
prob = ODEProblem(f, u0, tspan, p)
sol = solve(prob, Rodas5P(), callback = cbs, reltol=1e-9,abstol=1e-9#= dt = 1e-3, adaptive = false =#)
plot(sol)

[contradict]

After playing with this example a little bit, it looks like the problem has something to do with one callback creating the condition for the other. As soon as the switch opens, the voltage rises above threshold (in one integrator timestep). It looks like the diode condition3 is not re-checked after affect2!. This persists even if the clocking scheme (cb1 and cb2) here is replaced with PresetTimeCallbacks.

[ChrisRackauckas]

Doing this with a VectorContinuousCallback is fine?

[mongibellili]

Doing this with a VectorContinuousCallback is fine?

 function f(du, u, p, t)
        C = 1e-4; L = 1e-4; R = 10;U = 24.0; T = 1e-4; DC = 0.5; ROn = 1e-5;ROff = 1e5;
        rd=p[1];rs=p[2];
        il=u[1] ;uc=u[2]
        id=(il*rs-U)/(rd+rs) # diode's current
        du[1] =(-id*rd-uc)/L # inductor's current
        du[2]=(il-uc/R)/C    # capacitor's voltage
    end
    function condition(out, u, t, integrator) # Event when condition(out,u,t,integrator) == 0
        out[1] = (t-p[3])
        out[2] =  (t-p[4]-0.5*1e-4)
        out[3] = (p[5]*((u[1]*p[2]-24.0)/(p[1]+p[2]))+(1.0-p[5])*((u[1]*p[2]-24.0)*p[1]/(p[1]+p[2])))
        out[4] = -(p[5]*((u[1]*p[2]-24.0)/(p[1]+p[2]))+(1.0-p[5])*((u[1]*p[2]-24.0)*p[1]/(p[1]+p[2])))
    end
  
    function affect!(integrator, idx)
        if idx == 1
            p[4]=p[3]
            p[3]=p[3]+1e-4
            p[2]=1e-5
        elseif idx == 2
            p[2]=1e5
        elseif idx == 3
            p[1]=1e-5
            p[5]=1.0
        elseif idx == 4
            p[1]=1e5
            p[5]=0.0
        end
    end

  
    cbs = VectorContinuousCallback(condition, affect!, 4)
    u0 = [0.0, 0.0]
    tspan = (0.0, 0.0002)
    p = [1e5,1e-5,1e-4,0.0,0.0,0.0]
    prob = ODEProblem(f, u0, tspan, p)
    sol = solve(prob, Rodas5P(), callback = cbs, reltol=1e-4,abstol=1e-6,dtmax=1e-9,maxiters=Int(1e12)#= dt = 1e-3, adaptive = false =#)

ERROR: LoadError: Double callback crossing floating pointer reducer errored. Report this issue.
Stacktrace:
[1] error(s::String)
@ Base .\error.jl:35
[2] find_callback_time(integrator::OrdinaryDiffEq.ODEIntegrator{Rodas5P{1, true, Nothing, typeof(OrdinaryDiffEq.DEFAULT_PRECS), Val{:forward}, true, nothing}, true, Vector{Float64}, Nothing, Float64, Vector{Float64}, Float64, Float64, Float64, Float64, Vector{Vector{Float64}}, ODESolution{Float64, 2, Vector{Vector{Float64}}, Nothing, Nothing, Vector{Float64}, Vector{Vector{Vector{Float64}}}, ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Vector{Float64}, ......
more error details are in attached file.
errordetails.txt

[ChrisRackauckas]

Do you have multiple events at the same time with this?

[mongibellili]

Yes. There are 2 events at the same zero crossing function. 1 event for the rising transition and 1 event for the falling transition.

[ChrisRackauckas]

Okay yes, that case is not handled right now. We should get to it soon, but currently it's not soon enough. SciML/DiffEqBase.jl#519 is tracking this.