DiffEq skips over event in the first `solve` but works after `reinit!`

cuihantao · November 28, 2022, 5:08am

The code below simulates a series resistor and inductor supplied by an AC source through a diode. The diode turns off when the current is <=0 and turns back on when the voltage across the diode is >=1.

Skipping over events

My first issue with the code is that it skips over a turn-on event the first time it’s run. The second run yields the expected result.

using DifferentialEquations
using Plots
using BenchmarkTools

function I_load!(du, u, p, t)
    du[1] = (sqrt(2) * 230 * sin(2*pi*50.0*t) - u[1] * p[1]) / 0.5
end

u0 = [0.0]
p = [100, 1] # resistance, diode flag
tspan = (0.0, 60.0E-3)


# test if diode needs to be turned off
function condition_off(u, t, integrator)
    # println("test affect off at $(integrator.t), u=$(integrator.u)")
    if integrator.p[2] == 1
        u[1]
    else
        return 1
    end
    
end


function affect_off!(integrator)
    # println("applied affect off at t=$(integrator.t)")
    integrator.p[1] = 1e10
    integrator.p[2] = 0
end

# test if diode needs to be turned on
function condition_on(u, t, integrator)
    if integrator.p[2] == 0
        # @show integrator.t
        # @show (sqrt(2)*230*sin(2*pi*50.0*t) - (0.5 * get_du(integrator)[1] + u[1] * 100)) - 1
        return (sqrt(2)*230*sin(2*pi*50.0*t) - (0.5 * get_du(integrator)[1] + u[1] * 100)) - 1
    else
        return -1
    end
end


function affect_on!(integrator)
    # println("applied affect on at t=$(integrator.t)")
    integrator.p[1] = 100
    integrator.p[2] = 1
end

cb_off = ContinuousCallback(condition_off, nothing, affect_off!)  # only when condition goes negative
cb_on = ContinuousCallback(condition_on, affect_on!, nothing)  # only when condition goes positive
cb_set = CallbackSet(cb_off, cb_on);

prob = ODEProblem(I_load!, u0, tspan, p)

solver = AutoTsit5(Rosenbrock23())

integrator = init(prob, solver,
    callback=cb_set,
);
sol = solve!(integrator)
plot(sol)

Doing a reinit! yields the expected solution

reinit!(integrator, u0)
prob.p[:] = [100, 1.0];  # restore parameters
sol = solve!(integrator)

plot(sol.t, sol')

If you uncomment the two lines in function condition_on, the following output can be seen:

integrator.t = 0.03984942934492531
(sqrt(2) * 230 * sin(2 * pi * 50.0 * t) - (0.5 * (get_du(integrator))[1] + u[1] * 100)) - 1 = -14.524901788326165
integrator.t = 0.04421192005119001
(sqrt(2) * 230 * sin(2 * pi * 50.0 * t) - (0.5 * (get_du(integrator))[1] + u[1] * 100)) - 1 = 51.91978477325262

This means that the function reports a rising zero crossing, but the solve did not enter the time interpolation to detect the zero crossing. I spent several hours in troubleshooting my event functions but didn’t see obvious issues here.

Trapezoid() gives wrong results - Stiffness?

Another issue about the solver. If the solver is changed to Trapezoid() or AutoTsit5(Trapezoid()), the results are wrong.

AutoTsit5(Trapezoid()) plot using plot(sol.t, sol')

I speculate that the Trapezoid method also skips over events. With a really tiny step size, the solution is expected.

integrator = init(prob, solver,
    callback=cb_set,
    dt=5e-7, adaptive=false
);

This step size is obviously too small to be practical. The same problem is solved by a hand-baked trapezoid method in a textbook using a step size of 5e-5 (although that example is used to demonstrate numerical chattering). But if I use dtmax=5e-5, adaptive=true, the last turn-on event is again skipped, producing the same result as in the first plot.

Different Plots

Also, plot(sol.t, sol') and plot(sol) gives seemingly unrelated plots.

AutoTsit5(Trapezoid()) plot using plot(sol):

julia> maximum(sol.u)

1-element Vector{Float64}:
 2.0098042506364657

Please let me know if I missed something obvious. Thank you for your time.

cuihantao · November 28, 2022, 7:48pm

I tested the code with solver = Rosenbrock23(), and it works fine upon the first solve!.

I’m using DifferentialEquations v7.6.0.

ChrisRackauckas · November 29, 2022, 7:22pm

First thing, did you try lowering the tolerances?

That’s going to be ill-conditioned, so you’ll need a lower tolerance solution in order to resolve it well. My guess is that a lower tolerance is all you need, and that at higher tolerances you may need to adjust

https://diffeq.sciml.ai/stable/features/callback_functions/#SciMLBase.ContinuousCallback

repeat_nudge=1//100 to something like repeat_nudge=1//1000. But anyways,

integrator = init(prob, solver,
    callback=cb_set,
    abstol = 1e-8, reltol = 1e-8
);

is probably all that’s needed. And thn I wouldn’t use Rosenbrock23 for those kinds of tolerances.

cuihantao · November 29, 2022, 7:24pm

Thanks a lot, Chris. I was aware of the ill condition but did not understand that a low tolerance solution is achievable in such a condition. I will test it and report back.

cuihantao · December 2, 2022, 6:17pm

Lowering the tolerance and using a different solver works. It’s just that the total time is greater than a fixed-stepping method.

Oscar_Smith · December 2, 2022, 6:40pm

How does the accuracy compare? It’s really easy to be faster if you are a lot less accurate.

ChrisRackauckas · December 2, 2022, 7:21pm

And what about repeat_nudge=1//10?

cuihantao · December 2, 2022, 8:13pm

It does cut the total time by a lot. From destats, I can tell the number of root-finding calls has been reduced. Currently, I’m using

solver = Rodas4P()

integrator = init(prob, solver,
    callback=cb_set,
    abstol=1e-3,
    reltol = 1e-5,
);

with repeat_nudge=1//10.

The total time is now 116.616 μs (643 allocations: 30.81 KiB). The custom integrator takes about 2 us (without allocation). That’s not an apples-to-apples comparison but still two orders apart. I couldn’t identify allocations in my code. How could I further reduce allocations?

cuihantao · December 2, 2022, 8:16pm

The DiffEq results are more accurate right after the event, but I’m willing to trade reasonable accuracy for speed.

Below shows the results using a custom integrator (based on the trapezoidal rule) and using Rodas4P. What settings can I change to further improve the speed?

Oscar_Smith · December 2, 2022, 8:26pm

If you want to trade accuracy for speed, raise your tolerances. Rodas4P is also usually not the best solver. Rosenbrock23 will be probably be better with loose tolerances, and Rodas5P will usually be better at stricter tolerances.

cuihantao · December 2, 2022, 8:30pm

I understand that a lower tolerance leads to a higher number of iterations and thus slower speed.

Using Rosenbrock23, the results are not correct. It skips over the events. That’s why I posted the original problem.

repeat_nudge=1//10 and the following is what I used to get the figure above.

solver = Rosenbrock23()

integrator = init(prob, solver,
    callback=cb_set,
    abstol=1e-3,
    reltol = 1e-5,
);

ChrisRackauckas · December 3, 2022, 1:09am

If it’s missing an event, you may need to decrease it. Can you have it print out integrator.dt in the affect!? I think this may be a case where you have a bit of an adverse effect because the ODE might be simple enough that the error is essentially zero for any step, so the dt always grows, and so it grows so large that numerical errors become easy. In which case, fixing maxdt needs to be done. Show integrator.dt to see if that’s the case.

cuihantao · December 3, 2022, 3:28am

Thanks again! You are right, dt keeps increasing. Setting dtmax solves the problem. Below is the solver and tolerance I’m using now:

solver = Rodas5P()

integrator = init(prob, solver,
    callback=cb_set,
    abstol=1e-4,
    reltol = 1e-5,
    dtmax=1/60,
);

I need to build a larger problem to compare DiffEq’s solver time over a manual fixed-step variant of trapezoidal rule. Currently, DiffEq’s time is 122 us and the manual solver just takes 4 us. It’s probably because the manually written solver is non-allocating, and variables happen to fit into the cache.

The current destats are

DiffEqBase.DEStats
Number of function 1 evaluations:                  384
Number of function 2 evaluations:                  0
Number of W matrix evaluations:                    39
Number of linear solves:                           312
Number of Jacobians created:                       35
Number of nonlinear solver iterations:             0
Number of nonlinear solver convergence failures:   0
Number of rootfind condition calls:                843
Number of accepted steps:                          35
Number of rejected steps:                          4

There is a total of 645 allocations: 32.42 KiB. The optimization example here (https://diffeq.sciml.ai/stable/tutorials/faster_ode_example/#Example-Accelerating-a-Non-Stiff-Equation:-The-Lorenz-Equation) only has 41 allocations. Is my allocation number normal?

Also, since my actual problem is going to be large, I won’t be able to use the optimization techniques for small matrices.

Thank you again!

Oscar_Smith · December 3, 2022, 3:33am

You can use saveat to get preallocated output storage. That said, if your actual problem will be large, benchmarking small problems isn’t that useful. The larger problem will probably be stiffer, and require more of the machinery that the fixed timestep version is missing.

cuihantao · December 3, 2022, 3:37am

Noted. Thank you!!

ChrisRackauckas · December 3, 2022, 3:41am

A lot of those allocations are just saving the time series. I think there’s some more in the callback code that I need to eliminate though.

Oscar_Smith · December 3, 2022, 3:42am

Specifically, for large systems all the time is in jacobian solves and function evaluation. Allocating a length 3 array is roughly 1/3d as expensive as allocating a length 100 array, but the jacobians will be thousands of times slower.