From 4f1fe9bd72cc2fc2c3aabcc49c9132d25f14e387 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Mon, 27 May 2024 12:20:39 +0200 Subject: [PATCH 1/5] add ClockChange operator --- src/discretedomain.jl | 29 ++++++++++++++++++++++++++++- test/clock.jl | 21 +++++++++++++++++++++ 2 files changed, 49 insertions(+), 1 deletion(-) diff --git a/src/discretedomain.jl b/src/discretedomain.jl index 723612b083..30908f593e 100644 --- a/src/discretedomain.jl +++ b/src/discretedomain.jl @@ -147,6 +147,29 @@ Returns true if the expression or equation `O` contains [`Hold`](@ref) terms. """ hashold(O) = recursive_hasoperator(Hold, unwrap(O)) +# ClockChange + +""" +$(TYPEDEF) + +Change the clock of a discrete-time variable by sub or super sampling +``` +cont_x = ClockChange(from, to)(disc_x) +``` +""" +@kwdef struct ClockChange <: Operator + from::Any + to::Any +end +(D::ClockChange)(x) = Term{symtype(x)}(D, Any[x]) +(D::ClockChange)(x::Num) = Num(D(value(x))) +SymbolicUtils.promote_symtype(::ClockChange, x) = x + +function Base.:(==)(D1::ClockChange, D2::ClockChange) + isequal(D1.to, D2.to) && isequal(D1.from, D2.from) +end +Base.hash(D::ClockChange, u::UInt) = hash(D.from, hash(D.to, xor(u, 0xa5b640d6d952f101))) + # ShiftIndex """ @@ -242,10 +265,14 @@ function input_timedomain(h::Hold, arg = nothing) end output_timedomain(::Hold, arg = nothing) = Continuous() +input_timedomain(cc::ClockChange, arg = nothing) = cc.from +output_timedomain(cc::ClockChange, arg = nothing) = cc.to + sampletime(op::Sample, arg = nothing) = sampletime(op.clock) sampletime(op::ShiftIndex, arg = nothing) = sampletime(op.clock) +sampletime(op::ClockChange, arg = nothing) = sampletime(op.to) -changes_domain(op) = isoperator(op, Union{Sample, Hold}) +changes_domain(op) = isoperator(op, Union{Sample, Hold, ClockChange}) function output_timedomain(x) if isoperator(x, Operator) diff --git a/test/clock.jl b/test/clock.jl index f847f94188..963771d8a1 100644 --- a/test/clock.jl +++ b/test/clock.jl @@ -528,3 +528,24 @@ prob = ODEProblem(sys, [], (0.0, 10.0), [x(k - 1) => 2.0]) int = init(prob, Tsit5(); kwargshandle = KeywordArgSilent) @test int.ps[x] == 3.0 @test int.ps[x(k - 1)] == 2.0 + +## ClockChange +using ModelingToolkit: D_nounits as D +k1 = ShiftIndex(Clock(t, 1)) +k2 = ShiftIndex(Clock(t, 2)) +@mtkmodel CC begin + @variables begin + x(t) = 0 + y(t) = 0 + dummy(t) = 0 + end + @equations begin + x(k1) ~ x(k1 - 1) + + ModelingToolkit.ClockChange(from = k2.clock, to = k1.clock)(y(k2)) + y(k2) ~ y(k2 - 1) + 1 + D(dummy) ~ 0 + end +end +@mtkbuild cc = CC() +prob = ODEProblem(cc, [], (0, 10)) +sol = solve(prob, Tsit5(); kwargshandle = KeywordArgSilent) From 0badc165ad235cd3c2cf832ffd837544892af7b8 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Mon, 27 May 2024 14:56:51 +0200 Subject: [PATCH 2/5] handle new operator in shift_discrete_system --- src/systems/systemstructure.jl | 9 +++++---- test/clock.jl | 22 ++++++++++++++++++++++ 2 files changed, 27 insertions(+), 4 deletions(-) diff --git a/src/systems/systemstructure.jl b/src/systems/systemstructure.jl index 1c5f3ca28b..4991ad2096 100644 --- a/src/systems/systemstructure.jl +++ b/src/systems/systemstructure.jl @@ -457,23 +457,24 @@ function lower_order_var(dervar, t) end function shift_discrete_system(ts::TearingState) + dops = Union{Sample, Hold, ClockChange} @unpack fullvars, sys = ts discvars = OrderedSet() eqs = equations(sys) for eq in eqs - vars!(discvars, eq; op = Union{Sample, Hold}) + vars!(discvars, eq; op = dops) end iv = get_iv(sys) discmap = Dict(k => StructuralTransformations.simplify_shifts(Shift(iv, 1)(k)) for k in discvars - if any(isequal(k), fullvars) && !isa(operation(k), Union{Sample, Hold})) + if any(isequal(k), fullvars) && !isa(operation(k), dops)) for i in eachindex(fullvars) fullvars[i] = StructuralTransformations.simplify_shifts(fast_substitute( - fullvars[i], discmap; operator = Union{Sample, Hold})) + fullvars[i], discmap; operator = dops)) end for i in eachindex(eqs) eqs[i] = StructuralTransformations.simplify_shifts(fast_substitute( - eqs[i], discmap; operator = Union{Sample, Hold})) + eqs[i], discmap; operator = dops)) end @set! ts.sys.eqs = eqs @set! ts.fullvars = fullvars diff --git a/test/clock.jl b/test/clock.jl index 963771d8a1..4d4ab56816 100644 --- a/test/clock.jl +++ b/test/clock.jl @@ -549,3 +549,25 @@ end @mtkbuild cc = CC() prob = ODEProblem(cc, [], (0, 10)) sol = solve(prob, Tsit5(); kwargshandle = KeywordArgSilent) +d1 = reduce(vcat, sol.prob.kwargs[:disc_saved_values][1].saveval) +d2 = reduce(vcat, sol.prob.kwargs[:disc_saved_values][2].saveval) + +if length(d1) < length(d2) + d1, d2 = d2, d1 +end + +@test d2 == 1:6 # y + +X = [0, 0] # xy +ti = 0 +for ti in 0:10 + i = ti + 1 + x, y = X + if ti % 2 == 0 + y = y + 1 + end + x = x + y + + @test d1[i] == x + X = [x, y] +end From edd2f85b2ba4755cdb810529fcd36067937f50b2 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Tue, 28 May 2024 14:47:21 +0200 Subject: [PATCH 3/5] export operator --- src/ModelingToolkit.jl | 2 +- src/clock.jl | 2 +- src/discretedomain.jl | 4 ++-- src/systems/clock_inference.jl | 3 ++- src/systems/discrete_system/discrete_system.jl | 3 ++- 5 files changed, 8 insertions(+), 6 deletions(-) diff --git a/src/ModelingToolkit.jl b/src/ModelingToolkit.jl index 7d85f05336..a5da59e3d9 100644 --- a/src/ModelingToolkit.jl +++ b/src/ModelingToolkit.jl @@ -268,7 +268,7 @@ export debug_system #export Continuous, Discrete, sampletime, input_timedomain, output_timedomain #export has_discrete_domain, has_continuous_domain #export is_discrete_domain, is_continuous_domain, is_hybrid_domain -export Sample, Hold, Shift, ShiftIndex, sampletime, SampleTime +export Sample, Hold, Shift, ShiftIndex, sampletime, SampleTime, ClockChange export Clock #, InferredDiscrete, end # module diff --git a/src/clock.jl b/src/clock.jl index 7ca1707724..87539d81ea 100644 --- a/src/clock.jl +++ b/src/clock.jl @@ -55,7 +55,7 @@ See also [`is_discrete_domain`](@ref) """ function has_discrete_domain(x) issym(x) && return is_discrete_domain(x) - hasshift(x) || hassample(x) || hashold(x) + hasshift(x) || hassample(x) || hashold(x) || hasclockchange(x) end """ diff --git a/src/discretedomain.jl b/src/discretedomain.jl index 30908f593e..f6083b320d 100644 --- a/src/discretedomain.jl +++ b/src/discretedomain.jl @@ -166,10 +166,10 @@ end SymbolicUtils.promote_symtype(::ClockChange, x) = x function Base.:(==)(D1::ClockChange, D2::ClockChange) - isequal(D1.to, D2.to) && isequal(D1.from, D2.from) + ==(D1.to, D2.to) && ==(D1.from, D2.from) end Base.hash(D::ClockChange, u::UInt) = hash(D.from, hash(D.to, xor(u, 0xa5b640d6d952f101))) - +hasclockchange(O) = recursive_hasoperator(ClockChange, unwrap(O)) # ShiftIndex """ diff --git a/src/systems/clock_inference.jl b/src/systems/clock_inference.jl index 66af67d026..112231f21e 100644 --- a/src/systems/clock_inference.jl +++ b/src/systems/clock_inference.jl @@ -250,9 +250,10 @@ function generate_discrete_affect( cont_to_disc_obs = build_explicit_observed_function( use_index_cache ? osys : syss[continuous_id], needed_cont_to_disc_obs, - throw = false, + throw = true, expression = true, output_type = SVector) + disc_to_cont_obs = build_explicit_observed_function(sys, needed_disc_to_cont_obs, throw = false, expression = true, diff --git a/src/systems/discrete_system/discrete_system.jl b/src/systems/discrete_system/discrete_system.jl index bd72c533d0..f8b5336c46 100644 --- a/src/systems/discrete_system/discrete_system.jl +++ b/src/systems/discrete_system/discrete_system.jl @@ -139,7 +139,8 @@ function DiscreteSystem(eqs::AbstractVector{<:Equation}, iv, dvs, ps; iv′ = value(iv) dvs′ = value.(dvs) ps′ = value.(ps) - if any(hasderiv, eqs) || any(hashold, eqs) || any(hassample, eqs) || any(hasdiff, eqs) + if any(hasderiv, eqs) || any(hashold, eqs) || any(hassample, eqs) || + any(hasdiff, eqs) || any(hasclockchange, eqs) error("Equations in a `DiscreteSystem` can only have `Shift` operators.") end if !(isempty(default_u0) && isempty(default_p)) From 811a6bb065c7660e08d202e9d6fa3b051d4c456a Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Wed, 29 May 2024 10:45:31 +0200 Subject: [PATCH 4/5] add manual test function --- test/clock.jl | 27 ++++++++++++++++----------- 1 file changed, 16 insertions(+), 11 deletions(-) diff --git a/test/clock.jl b/test/clock.jl index 4d4ab56816..c49f6d3578 100644 --- a/test/clock.jl +++ b/test/clock.jl @@ -558,16 +558,21 @@ end @test d2 == 1:6 # y -X = [0, 0] # xy -ti = 0 -for ti in 0:10 - i = ti + 1 - x, y = X - if ti % 2 == 0 - y = y + 1 +# Manual implementation of the dynamics +function manualtest() + x, y = 0, 0 + X = [[x, y]] # xy + ti = 0 + for ti in 0:10 + i = ti + 1 + if ti % 2 == 0 + y = y + 1 + end + x = x + y + + @test d1[i] == x + push!(X, [x, y]) end - x = x + y - - @test d1[i] == x - X = [x, y] + X end +X = manualtest() From 6c135d8c9021c8a59dab6cd37970f021f94e3bf1 Mon Sep 17 00:00:00 2001 From: Fredrik Bagge Carlson Date: Wed, 29 May 2024 10:45:51 +0200 Subject: [PATCH 5/5] add phase option to clock --- src/clock.jl | 17 ++++++++++------- src/systems/clock_inference.jl | 1 - src/systems/diffeqs/abstractodesystem.jl | 3 ++- test/clock.jl | 7 +++---- 4 files changed, 15 insertions(+), 13 deletions(-) diff --git a/src/clock.jl b/src/clock.jl index 87539d81ea..d12921b688 100644 --- a/src/clock.jl +++ b/src/clock.jl @@ -101,26 +101,29 @@ abstract type AbstractClock <: AbstractDiscrete end """ Clock <: AbstractClock - Clock([t]; dt) + Clock([t]; dt, phase = 0) The default periodic clock with independent variables `t` and tick interval `dt`. If `dt` is left unspecified, it will be inferred (if possible). +`phase` is an optional phase delay, such that the clock ticks at `(t_0 + phase : dt : t_final). """ struct Clock <: AbstractClock "Independent variable" t::Union{Nothing, Symbolic} "Period" dt::Union{Nothing, Float64} - Clock(t::Union{Num, Symbolic}, dt = nothing) = new(value(t), dt) - Clock(t::Nothing, dt = nothing) = new(t, dt) + phase::Float64 + Clock(t::Union{Num, Symbolic}, dt = nothing; phase = 0) = new(value(t), dt, phase) + Clock(t::Nothing, dt = nothing; phase = 0) = new(t, dt, phase) end -Clock(dt::Real) = Clock(nothing, dt) -Clock() = Clock(nothing, nothing) +Clock(dt::Real; phase = 0) = Clock(nothing, dt, phase) +Clock() = Clock(nothing, nothing, 0) sampletime(c) = isdefined(c, :dt) ? c.dt : nothing -Base.hash(c::Clock, seed::UInt) = hash(c.dt, seed ⊻ 0x953d7a9a18874b90) +Base.hash(c::Clock, seed::UInt) = hash(c.phase, hash(c.dt, seed ⊻ 0x953d7a9a18874b90)) function Base.:(==)(c1::Clock, c2::Clock) - ((c1.t === nothing || c2.t === nothing) || isequal(c1.t, c2.t)) && c1.dt == c2.dt + ((c1.t === nothing || c2.t === nothing) || isequal(c1.t, c2.t)) && c1.dt == c2.dt && + c1.phase == c2.phase end is_concrete_time_domain(x) = x isa Union{AbstractClock, Continuous} diff --git a/src/systems/clock_inference.jl b/src/systems/clock_inference.jl index 112231f21e..e5da514720 100644 --- a/src/systems/clock_inference.jl +++ b/src/systems/clock_inference.jl @@ -253,7 +253,6 @@ function generate_discrete_affect( throw = true, expression = true, output_type = SVector) - disc_to_cont_obs = build_explicit_observed_function(sys, needed_disc_to_cont_obs, throw = false, expression = true, diff --git a/src/systems/diffeqs/abstractodesystem.jl b/src/systems/diffeqs/abstractodesystem.jl index b976945f79..6b59f69004 100644 --- a/src/systems/diffeqs/abstractodesystem.jl +++ b/src/systems/diffeqs/abstractodesystem.jl @@ -1092,8 +1092,9 @@ function DiffEqBase.ODEProblem{iip, specialize}(sys::AbstractODESystem, u0map = affects, inits, clocks, svs = ModelingToolkit.generate_discrete_affect(sys, dss...) discrete_cbs = map(affects, clocks, svs) do affect, clock, sv if clock isa Clock + tsteps = (tspan[1] + clock.phase):(clock.dt):tspan[2] PeriodicCallback(DiscreteSaveAffect(affect, sv), clock.dt; - final_affect = true, initial_affect = true) + final_affect = tspan[2] ∈ tsteps, initial_affect = iszero(clock.phase), phase = clock.phase) elseif clock isa SolverStepClock affect = DiscreteSaveAffect(affect, sv) DiscreteCallback(Returns(true), affect, diff --git a/test/clock.jl b/test/clock.jl index c49f6d3578..b6da0d693c 100644 --- a/test/clock.jl +++ b/test/clock.jl @@ -531,7 +531,7 @@ int = init(prob, Tsit5(); kwargshandle = KeywordArgSilent) ## ClockChange using ModelingToolkit: D_nounits as D -k1 = ShiftIndex(Clock(t, 1)) +k1 = ShiftIndex(Clock(t, 1; phase = 100eps())) # TODO: this should not be required k2 = ShiftIndex(Clock(t, 2)) @mtkmodel CC begin @variables begin @@ -540,14 +540,13 @@ k2 = ShiftIndex(Clock(t, 2)) dummy(t) = 0 end @equations begin - x(k1) ~ x(k1 - 1) + - ModelingToolkit.ClockChange(from = k2.clock, to = k1.clock)(y(k2)) + x(k1) ~ x(k1 - 1) + ClockChange(from = k2.clock, to = k1.clock)(y(k2)) y(k2) ~ y(k2 - 1) + 1 D(dummy) ~ 0 end end @mtkbuild cc = CC() -prob = ODEProblem(cc, [], (0, 10)) +prob = ODEProblem(cc, [], (0, 10.1)) # TODO: delaying end point by 0.1 should not be required without phase shift sol = solve(prob, Tsit5(); kwargshandle = KeywordArgSilent) d1 = reduce(vcat, sol.prob.kwargs[:disc_saved_values][1].saveval) d2 = reduce(vcat, sol.prob.kwargs[:disc_saved_values][2].saveval)