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?