Solver gets trapped when using callbacks

[BastienDassonneville]

Using Julia 1.5.3 and DifferentialEquations v6.15.0, the solver gets stuck at a time value and won't progress anymore. I'm using a VectorContinuousCallback to determine events that are root findable, defined in the condition function. There are additional boolean conditions that need to be met for the events to actually occur. I put these conditions in the affect function. To reproduce that error, you can use the following code snippet:

const torque = 1.0
const Ic = 10.0
const Iv = 0.15
const r_c = 1.0
const r_v = 0.3
const αc = 24.0 * pi/180
const e_r = 0.05
const nteeth = round(Int, 2pi/αc)
const Mv = Iv/r_v^2
const Mc = Ic/r_c^2
const Gc = Mv*(1 + e_r)/r_c/(Mv + Mc)
const Gv = Mc*(1 + e_r)/r_v/(Mv + Mc)
const αv = Gv*r_c/(1-e_r)*αc/2 

p = (torque, Ic)

function f(du,u,p,t)
    du[1] = u[3]
    du[2] = u[4]
    du[3] = p[1]/p[2] # p =(torque, crown inertia)
    du[4] = 0.0
    @show t
end

#out = vector of length 2nteeth
function condition(out,u,t,integrator)     
    for m in 1:nteeth
        out[m] = r_c*sin(u[1] - (m-1)*αc) - r_v*tan(u[2] + αv/2)#upper
        out[m + nteeth] = r_c*sin((m-1)*αc - u[1]) - r_v*tan(-u[2] + αv/2) #lower
    end
end

function affect!(integrator, m)
    c_up = r_c*integrator.u[3] - r_v*integrator.u[4] > 0.0 && ( (m-1)*αc <= mod(integrator.u[1] + αc/2, 2π) <= m*αc)#upper
    c_dn = r_c*integrator.u[3] + r_v*integrator.u[4] > 0.0 && ( (m-1)*αc <= mod(integrator.u[1] + αc/2 + π, 2π) <= m*αc ) #lower
    if  m <= nteeth  && c_up
        integrator.u[3] = -r_c*Gc*integrator.u[3] + r_v*Gc*integrator.u[4]
        integrator.u[4] = + r_c*Gv*integrator.u[3] - r_v*Gv*integrator.u[4] 
    elseif m > nteeth && c_dn
        integrator.u[3] =  -r_c*Gc*integrator.u[3] - r_v*Gc*integrator.u[4]
        integrator.u[4] =  -r_c*Gv*integrator.u[3] - r_v*Gv*integrator.u[4] 
    end
end

cb = VectorContinuousCallback(condition,affect!,2nteeth)

u0 = [0, 0.0, 0.0, 3.0]
prob = ODEProblem(f,u0,(0.0,2.0),p)
sol = solve(prob, callback=cb, abstol=1e-3, reltol=1e-3, dt = 1e-1, adaptive=false) #
plot(sol,layout=(4,1)) 

Note that I would get an "Event repeated at same time" error when using v6.12.0. If this is of any help, I'm trying to reproduce results published in Roup et al.

[ChrisRackauckas]

@kanav99 can you take a look at this?

[kanav99]

Sure, I will look into it

[kanav99]

@BastienDassonneville looks like the condition is not continuous due to tan. If you can make sure that the interval (dt) always contains a continuous branch of tan, we might be able to do something.

[BastienDassonneville]

Thanks @kanav99. The dynamics of the system should keep -π/2<u[2]+αv/2<π/2 : I expect that the first event pushes u[2]+αv/2 towards π/2 and the following event should go the other way round and bring u[2]+αv/2 towards -π/2. I should enforce that condition in the program but I'm not sure how to do it. I tried to replace u[2]+αv/2 by sign(u[2] + αv/2)*mod(u[2] + αv/2, pi/2) (and similarly for the -u[2]+αv/2 term) but that didn't fix it...
I'm not sure I understand where the problem is coming from though and I'd appreciate if you could clarify how the discontinuity is making the algorithm fail.

[ChrisRackauckas]

The callbacks are a rootfinding position, searching for the time at which condition=0. Your condition goes directly from +Inf to -Inf, so it doesn't have a zero.

[BastienDassonneville]

Thank you for your feedback, massaging the condition to avoid the discontinuity of tan (simply multiplying by cos) fixed the issue. The condition function that works reads:

function condition(out,u,t,integrator)     
    for m in 1:nteeth
        out[m] = r_c*sin(u[1] - (m-1)*αc)*cos(u[2] + αv/2) - r_v*sin(u[2] + αv/2)#upper
        out[m + nteeth] = r_c*sin((m-1)*αc - u[1])*cos(-u[2] + αv/2) - r_v*sin(-u[2] + αv/2) #lower
    end
end

[mforets]

Hi @BastienDassonneville,

I got interested in the model that you wrote, thanks for sharing. (For my research I write algorithms for hybrid systems.)

The condition function that works reads:

Sorry, isn't that function the same as the one in your original post? Perhaps you pasted the wrong function? Were you able to recover plots similar to those in Fig. 6 of the paper?

[BastienDassonneville]

oops, I copied-paste the wrong condition function. The correct one is:

function condition(out,u,t,integrator)     
    for m in 1:nteeth
        out[m] = r_c*sin(u[1] - (m-1)*αc)*cos(u[2] + αv/2) - r_v*sin(u[2] + αv/2)#upper
        out[m + nteeth] = r_c*sin((m-1)*αc - u[1])*cos(-u[2] + αv/2) - r_v*sin(-u[2] + αv/2) #lower
    end
end

I edited the previous post to fix that mistake.

Using that condition function, you can find that you reach a limit cycle that looks like a bow tie (plot(sol,vars=(3,4)) yields something similar to fig. 5 of Roup et al). When you look at the time evolution, it looks like it takes longer to reach the permanent regime (more like after 120s) using the following conditions:

u0 = [0, 0.0, 0.0, 3.0]
prob = ODEProblem(f,u0,(0.0,400.0),p)
sol = solve(prob, callback=cb, dt = 1e-2, adaptive=false)

On fig.6 it looks like it's happening on a scale of a few seconds. Not sure yet why it'd take longer, let me know if you have any insight.

[mforets]

out[m] = r_c*sin(u[1] - (m-1)αc)r_vcos(u[2] + αv/2) - r_vsin(u[2] + αv/2)#upper

i think that should be

out[m] = r_c*sin(u[1] - (m-1)*αc)cos(u[2] + αv/2) - r_vsin(u[2] + αv/2)#upper

(same in the term for out[m + nteeth])

for the crown gear velocity (plot(sol, vars=(0, 3))) i get a saw-tooth behavior similar to Fig. 7 of the paper but with a different time scale. (.. i have double checked the rest of the model and constant parameter values used in Julia and they lgtm.)

[BastienDassonneville]

you're totally right, I mistakenly multiplied the first term by r_v (that confirms that I'm not really awake this morning ;) ). I updated the previous posts with the correct condition. Thanks for double checking the model and confirming that no other typos found their way in the code.

[devmotion]

BTW you might be interested in sincos which computes sin and cos of the same argument more efficiently in one function.

[BastienDassonneville]

Thanks @devmotion for bringing sincos to my attention: I got a ~3X speedup (evaluating using @btime from BenchmarkTools) by rewriting the condition as

function condition(out,u,t,integrator) 
    SCup = sincos(u[2] + αv/2)
    SCdn = sincos(-u[2] + αv/2)
    for m in 1:nteeth
        S = sin(u[1] - (m-1)*αc)       
        out[m] = r_c*S*SCup[2] - r_v*SCup[1]#upper
        out[m + nteeth] = -r_c*S*SCdn[2] - r_v*SCdn[1] #lower
    end
end

[BastienDassonneville]

@mforets I realized I made 2 mistakes in the affect function:

1. in the indexing for the definition of c_dn, I forgot to subtract an nteeth
2. when the conditions are met, it should be a += and not an =

This gives a timescale in agreement with the article. Note that it converges toward a different verge velocity (0.647 instead of 0.813) but if you look at an animation, it looks physically realistic (in contrast with the previously posted code). Not sure where the difference is coming from though... Here's the code:

const torque = 1.0
const Ic = 10.0
const Iv = 0.15
const r_c = 1.0
const r_v = 0.3
const αc = 24.0 * pi/180
const e_r = 0.05
const nteeth = round(Int, 2pi/αc)
const Mv = Iv/r_v^2
const Mc = Ic/r_c^2
const Gc = Mv*(1 + e_r)/r_c/(Mv + Mc)
const Gv = Mc*(1 + e_r)/r_v/(Mv + Mc)
const αv =  Gv*r_c/(1-e_r)*αc/2 

p = (torque, Ic)

function f(du,u,p,t)
    du[1] = u[3]
    du[2] = u[4]
    du[3] = p[1]/p[2] # p =(torque, crown inertia)
    du[4] = 0.0
    
end

#out = vector of length 2nteeth
function condition(out,u,t,integrator) 
    SCup = sincos(u[2] + αv/2)
    SCdn = sincos(-u[2] + αv/2)
    for m in 1:nteeth
        S = sin(u[1] - (m-1)*αc)       
        out[m] = r_c*S*SCup[2] - r_v*SCup[1]#upper
        out[m + nteeth] = -r_c*S*SCdn[2] - r_v*SCdn[1] #lower 
    end
end

function affect!(integrator, i)
    # i is the index of the out vector and runs from 1 to 2nteeth
    # We chose to offset the index by nteeth for the lower paddle in the condition function.
    # We should keep track of that offset in the remaining conditions for the lower paddle.
    c_up = r_c*integrator.u[3] - r_v*integrator.u[4] > 0.0 && ( (i-1)*αc <= mod(integrator.u[1] + αc/2, 2π) <= i*αc)#upper   
    c_dn = r_c*integrator.u[3] + r_v*integrator.u[4] > 0.0 && ( (i-1-nteeth)*αc <= mod(integrator.u[1] + αc/2 + π, 2π) <= (i-nteeth)*αc ) #lower

    if  i <= nteeth  && c_up
        #  `+=` because there is a transfert of momentum from the paddle to the tooth
        integrator.u[3] += -r_c*Gc*integrator.u[3] + r_v*Gc*integrator.u[4] 
        integrator.u[4] += r_c*Gv*integrator.u[3] - r_v*Gv*integrator.u[4]
        
    elseif i > nteeth && c_dn
        integrator.u[3] +=  -r_c*Gc*integrator.u[3] - r_v*Gc*integrator.u[4]
        integrator.u[4] +=  -r_c*Gv*integrator.u[3] - r_v*Gv*integrator.u[4] 
        
    end
end

cb = VectorContinuousCallback(condition,affect!,2nteeth)

u0 = [0.0, 0.0, 0.0, 3.0]
prob = ODEProblem(f,u0,(0.0,60.0),p)
sol = solve(prob, callback=cb, dt = 1e-1, adaptive=false) 
plot(sol,layout=(4,1)) 

for an animation:

colors = [m == 0 ? "black" : "gray" for m in 0:nteeth-1]
@gif for i in 1:length(sol.t)
    θc, θv, ωc, ωv = sol.u[i]
    teeth = [θc + m*αc for m in 0:nteeth-1] 
    scatter(sin.(teeth), cos.(teeth), aspect_ratio=:equal, color= colors, xlims = (-1.25, 1.25), ylims = (-1.25, 1.25))
    
    plot!([r_v*tan(θv + αv/2), r_v*tan(θv + αv/2)], [0.9, 1.1], aspect_ratio=:equal, color= "red")
    plot!([r_v*tan(θv - αv/2), r_v*tan(θv - αv/2)], [-0.9, -1.1], aspect_ratio=:equal, color= "red")
    annotate!(1.2, 1.0, Plots.text("$(sol.t[i])"))
    plot!(legend=nothing)
    end  

[ChrisRackauckas]

Try a lower tolerance?