Miss executions of ContinuousCallbacks due to simultaneous timing and numerical precision issues

[sablonl]

I've encountered two related bugs with continuous events in ModelingToolkit/DifferentialEquations:
(I know this is a very bad way of doing Periodic/Discrete Callbacks but I guess that issue could occur in different cases)

Bug 1: Only first event executes when multiple events occur simultaneously

When multiple continuous events are scheduled to trigger at exactly the same time, only the first event in the array executes. The subsequent events that should fire at the same time are skipped entirely.

Expected behavior: All events that satisfy their trigger conditions at the same time should execute.

Actual behavior: Only the first event executes when multiple events have simultaneous trigger times.

Bug 2: Events skipped due to numerical precision in modulo operations

Some continuous events are being skipped when using modulo operations in event conditions, likely due to floating-point precision issues. The t % frequency ≈ 0 condition fails to trigger even when it should mathematically be true.

Expected behavior: Events should trigger reliably at their specified intervals, accounting for floating-point precision limitations.

MWE:

using Pkg
Pkg.activate(temp = true)

Pkg.add(["ModelingToolkit", "DifferentialEquations"])

using ModelingToolkit, DifferentialEquations
using ModelingToolkit: D_nounits as D, t_nounits as t

@variables begin
	# Define the variables used in the equations
	x(t)::Float64
	dum_t_a(t)::Float64
	dum_t_b(t)::Float64
end

eqs = [D(x) ~ -x
	D(dum_t_a) ~ 0
	D(dum_t_b) ~ 0] # Initial equations

function cb_a(_t)
	"Callback A executed at time: $_t" |> println
	return _t
end
@register_symbolic cb_a(_t)
function cb_b(_t)
	"Callback B executed at time: $_t" |> println
	return _t
end
@register_symbolic cb_b(_t)

tstops(frequency_a, frequency_b) = sort(unique(union(
	range(0.0, 10.0, step = frequency_a),
	range(0.0, 10.0, step = frequency_b),
)))
function test(frequency_a, frequency_b)
	evs = Any[]
	push!(evs,
		[t % frequency_a ~ 0] => [dum_t_a ~ cb_a(Pre(t))],
	)
	push!(evs,
		[t % frequency_b ~ 0] => [dum_t_b ~ cb_b(Pre(t))],
	)

	@named sys = ODESystem(eqs, t;
		continuous_events = evs)

	compiled_sys = mtkcompile(sys)

	ode_prob = ODEProblem(compiled_sys, [x => 1.0, dum_t_a => 0.0, dum_t_b => 0.0], (0.0, 10.0))

	sol = solve(ode_prob;
		tstops = tstops(frequency_a, frequency_b),
	);
	return nothing
end

function truth(frequency_a, frequency_b)
	vtstops = tstops(frequency_a, frequency_b)
	for t in vtstops
		if isapprox(t % frequency_a, 0, atol = 1e-6) # Because of reasons, t % frequency_a ≈ 0 does not work 
			println("Callback A executed at time: $t")
		end
		if isapprox(t % frequency_b, 0, atol = 1e-6)
			println("Callback B executed at time: $t")
		end
	end
end

test(1.0, 1.0)
truth(1.0, 1.0)
# output:
#       test(1.0, 1.0)                     |       truth(1.0, 1.0)

#                                          | Callback A executed at time: 0.0
#                                          | Callback B executed at time: 0.0
# Callback A executed at time: 1.0         | Callback A executed at time: 1.0
#                                          | Callback B executed at time: 1.0
# Callback A executed at time: 2.0         | Callback A executed at time: 2.0
#                                          | Callback B executed at time: 2.0
# Callback A executed at time: 3.0         | Callback A executed at time: 3.0
#                                          | Callback B executed at time: 3.0
# Callback A executed at time: 4.0         | Callback A executed at time: 4.0
#                                          | Callback B executed at time: 4.0
# Callback A executed at time: 5.0         | Callback A executed at time: 5.0
#                                          | Callback B executed at time: 5.0
# Callback A executed at time: 6.0         | Callback A executed at time: 6.0
#                                          | Callback B executed at time: 6.0
# Callback A executed at time: 7.0         | Callback A executed at time: 7.0
#                                          | Callback B executed at time: 7.0
# Callback A executed at time: 8.0         | Callback A executed at time: 8.0
#                                          | Callback B executed at time: 8.0
# Callback A executed at time: 9.0         | Callback A executed at time: 9.0
#                                          | Callback B executed at time: 9.0
# Callback A executed at time: 10.0        | Callback A executed at time: 10.0
#                                          | Callback B executed at time: 10.0

test(1.0, 1.15)
truth(1.0, 1.15)
# output:
#       test(1.0, 1.15)                    |       truth(1.0, 1.15)

#                                          | Callback A executed at time: 0.0
#                                          | Callback B executed at time: 0.0
# Callback A executed at time: 1.0         | Callback A executed at time: 1.0
# Callback B executed at time: 1.15        | Callback B executed at time: 1.15
# Callback A executed at time: 2.0         | Callback A executed at time: 2.0
# Callback B executed at time: 2.3         | Callback B executed at time: 2.3
# Callback A executed at time: 3.0         | Callback A executed at time: 3.0
#                                          | Callback B executed at time: 3.45
# Callback A executed at time: 4.0         | Callback A executed at time: 4.0
# Callback B executed at time: 4.6         | Callback B executed at time: 4.6
# Callback A executed at time: 5.0         | Callback A executed at time: 5.0
#                                          | Callback B executed at time: 5.75
# Callback A executed at time: 6.0         | Callback A executed at time: 6.0
#                                          | Callback B executed at time: 6.9
# Callback A executed at time: 7.0         | Callback A executed at time: 7.0
# Callback A executed at time: 8.0         | Callback A executed at time: 8.0
#                                          | Callback B executed at time: 8.05
# Callback A executed at time: 9.0         | Callback A executed at time: 9.0
# Callback B executed at time: 9.2         | Callback B executed at time: 9.2
# Callback A executed at time: 10.0        | Callback A executed at time: 10.0

Environment:

julia> versioninfo()
Julia Version 1.11.5
Commit 760b2e5b739 (2025-04-14 06:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 8 × 11th Gen Intel(R) Core(TM) i7-1165G7 @ 2.80GHz
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, tigerlake)
Threads: 1 default, 0 interactive, 1 GC (on 8 virtual cores)
Environment:
  JULIA_NUM_THREAD = 16
  JULIA_EDITOR = code
  JULIA_VSCODE_REPL = 1

julia> Pkg.status()
Status `/tmp/jl_wLkfu4/Project.toml`
  [0c46a032] DifferentialEquations v7.16.1
  [961ee093] ModelingToolkit v10.10.0

[ChrisRackauckas]

The t % frequency ≈ 0 condition fails to trigger even when it should mathematically be true.

This isn't going to work because it's not a zero crossing function. The modulus function never crosses zero, so the condition that sol(t) > 0 and sol(t + theta*dt) < 0 is never true, so there's no event to trigger. There's nothing floating point about this, just by the definition of the event handling functions you have a guarantee that this won't trigger.

Bug 1: Only first event executes when multiple events occur simultaneously

This is a duplicate of SciML/DiffEqBase.jl#519. It's not hard to fix but it's known to not be handled right now.

[AayushSabharwal]

If you want periodic events at a set rate, prefer to use symbolic discrete callbacks. This is how the event should be phrased:

frequency_a => [dum_t_a ~ cb_a(Pre(t))]

On a sidenote, Pre(t) is also not really a thing.