Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

callbacks_buckExample #998

Open
mongibellili opened this issue Nov 15, 2023 · 15 comments
Open

callbacks_buckExample #998

mongibellili opened this issue Nov 15, 2023 · 15 comments

Comments

@mongibellili
Copy link

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:
Tsit5()34_buck_ft0005

However, after simulating the system via 3 other tools:

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
buck_Csolver
ltspice
buck1ms
QS_Solver.jl
plot_buck_nmLiqss2_0 0001_ _ft_0 001_2023_11_14_16_13_44
Am I misusing callbacks?
thank you.

@oscardssmith
Copy link
Contributor

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

@mongibellili
Copy link
Author

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

@pepijndevos
Copy link

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
Copy link
Author

mongibellili commented Nov 20, 2023

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
Copy link

What do you mean "runs infinitely"?

@mongibellili
Copy link
Author

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
Copy link

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
Copy link

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.

Rodas5P(; linsolve = nothing, precs = DEFAULT_PRECS,)34_buck_ft0005

@pepijndevos
Copy link

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
Copy link

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
Copy link
Member

Doing this with a VectorContinuousCallback is fine?

@mongibellili
Copy link
Author

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
Copy link
Member

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

@mongibellili
Copy link
Author

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
Copy link
Member

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

5 participants