diff --git a/src/clock.jl b/src/clock.jl index 292577ee49..67f6197d33 100644 --- a/src/clock.jl +++ b/src/clock.jl @@ -98,11 +98,14 @@ end abstract type AbstractClock <: AbstractDiscrete end """ - Clock <: AbstractClock - Clock([t]; dt) + Clock() + Clock(dt) + Clock(t, dt) The default periodic clock with independent variables `t` and tick interval `dt`. If `dt` is left unspecified, it will be inferred (if possible). + +See also [`RealtimeClock`](@ref). """ struct Clock <: AbstractClock "Independent variable" @@ -122,3 +125,25 @@ function Base.:(==)(c1::Clock, c2::Clock) end is_concrete_time_domain(x) = x isa Union{AbstractClock, Continuous} + +""" + RealtimeClock() + RealtimeClock(dt) + RealtimeClock(t, dt) + +Similar to [`Clock`](@ref), but with with the additional property that the simulation is run no faster than real-time, i.e., a simulation with `tspan = (0, 10)` will take at least 10 seconds of wallclock time to run. + +To achieve real-time execution, the wall-clock duration of each integration step is measured, and a call to `Libc.systemsleep` is made to ensure that the next integration step is not started before `dt` seconds of wall-clock time has elapsed. This leads to a simulation that takes at least as long as the length of the time span, but may take longer if the computational load is high enough that one integration step takes more than `dt` seconds. +""" +struct RealtimeClock <: AbstractClock + clock::Clock + RealtimeClock(args...) = new(Clock(args...)) +end + +sampletime(c) = sampletime(c.clock) +Base.hash(c::RealtimeClock, seed::UInt) = hash(c.dt, seed ⊻ 0x9d3d7a9a18874b90) +Base.:(==)(c1::RealtimeClock, c2::RealtimeClock) = c1.clock == c2.clock +function Base.getproperty(c::RealtimeClock, s::Symbol) + s ∈ fieldnames(typeof(c)) && return getfield(c, s) + getproperty(getfield(c, :clock), s) +end diff --git a/src/systems/clock_inference.jl b/src/systems/clock_inference.jl index 9ffeb56319..a8190b8c88 100644 --- a/src/systems/clock_inference.jl +++ b/src/systems/clock_inference.jl @@ -236,3 +236,28 @@ function generate_discrete_affect(syss, inputs, continuous_id, id_to_clock; defaults = Dict{Any, Any}(v => 0.0 for v in Iterators.flatten(inputs)) return affects, clocks, svs, appended_parameters, defaults end + +function generate_discrete_callback(affect, clock, sv) + error("$clock is not a supported clock type.") +end # Fallback + +function generate_discrete_callback(affect, clock::Clock, sv) + PeriodicCallback(DiscreteSaveAffect(affect, sv), clock.dt) +end + +function generate_discrete_callback(affect, clock::RealtimeClock, sv) + start_time = time() # Set start time before simulation starts + function realtime_affect!(integrator, saved_values) + # This is a closure that closes over start_time + affect(integrator, saved_values) + execution_time = time() - start_time + sleep_time = clock.dt - execution_time + if sleep_time < 0 + @debug "RealtimeClock failed to meet deadline dt = $(clock.dt), execution time = $execution_time at simulation time t = $(t)" + end + Libc.systemsleep(max(0, sleep_time)) + start_time = time() # Update start time to current time + nothing + end + PeriodicCallback(DiscreteSaveAffect(realtime_affect!, sv), clock.dt) +end diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index 88a1ed9bb3..85da4829af 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -947,13 +947,7 @@ function DiffEqBase.ODEProblem{iip, specialize}(sys::AbstractODESystem, u0map = cbs = process_events(sys; callback, has_difference, kwargs...) if has_discrete_subsystems(sys) && (dss = get_discrete_subsystems(sys)) !== nothing affects, clocks, svs = ModelingToolkit.generate_discrete_affect(dss...) - discrete_cbs = map(affects, clocks, svs) do affect, clock, sv - if clock isa Clock - PeriodicCallback(DiscreteSaveAffect(affect, sv), clock.dt) - else - error("$clock is not a supported clock type.") - end - end + discrete_cbs = map(ModelingToolkit.generate_discrete_callback, affects, clocks, svs) if cbs === nothing if length(discrete_cbs) == 1 cbs = only(discrete_cbs) diff --git a/test/clock.jl b/test/clock.jl index 74946d21d7..f80c325a9c 100644 --- a/test/clock.jl +++ b/test/clock.jl @@ -310,3 +310,26 @@ if VERSION >= v"1.7" @test sol.u ≈ sol2.u end + +# RealtimeClock +dt = 0.01 +@variables t x(t)=0 y(t)=0 u(t)=0 yd(t)=0 ud(t)=0 r(t)=1 +@parameters kp = 1 +D = Differential(t) +d = ModelingToolkit.RealtimeClock(dt) +# u(n + 1) := f(u(n)) + +eqs = [yd ~ Sample(d)(y) + ud ~ kp * (r - yd) + r ~ 1.0 + +# plant (time continuous part) + u ~ Hold(ud) + D(x) ~ -x + u + y ~ x] +@named sys = ODESystem(eqs) +prob = ODEProblem(structural_simplify(sys), [], (0.0, 3.0)) +solve(prob, Tsit5(), kwargshandle = KeywordArgSilent) +execution_time = @elapsed solve(prob, Tsit5(), kwargshandle = KeywordArgSilent) +@test_skip execution_time≈3 atol=10 * dt # Very crude test that can be evaluated locally but not on CI +@test execution_time >= 3 - dt # This can be evaluated on CI as well since it only tests a lower bound