Add phase option to Clock

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

src/clock.jl

@@ -55,7 +55,7 @@
 """
 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
 
 """
@@ -101,26 +101,29 @@
 
 """
     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}
@@ -130,7 +133,7 @@
     SolverStepClock()
     SolverStepClock(t)
 
 A clock that ticks at each solver step (sometimes referred to as "continuous sample time"). This clock **does generally not have equidistant tick intervals**, instead, the tick interval depends on the adaptive step-size slection of the continuous solver, as well as any continuous event handling. If adaptivity of the solver is turned off and there are no continuous events, the tick interval will be given by the fixed solver time step `dt`. 
 
 Due to possibly non-equidistant tick intervals, this clock should typically not be used with discrete-time systems that assume a fixed sample time, such as PID controllers and digital filters.
 """

src/discretedomain.jl

@@ -147,6 +147,29 @@
 """
 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)
+    ==(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
 
 """
@@ -242,10 +265,14 @@
 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)

src/systems/clock_inference.jl

@@ -21,9 +21,9 @@
     ClockInference(ts, eq_domain, var_domain, inferred)
 end
 
 struct NotInferedTimeDomain end
 function error_sample_time(eq)
     error("$eq\ncontains `SampleTime` but it is not an infered discrete equation.")
 end
 function substitute_sample_time(ci::ClockInference, ts::TearingState)
     @unpack eq_domain = ci
@@ -34,7 +34,7 @@
         domain = eq_domain[i]
         dt = sampletime(domain)
         neweq = substitute_sample_time(eq, dt)
         if neweq isa NotInferedTimeDomain
             error_sample_time(eq)
         end
         eqs[i] = neweq
@@ -52,13 +52,13 @@
     op = operation(ex)
     args = arguments(ex)
     if op == SampleTime
         dt === nothing && return NotInferedTimeDomain()
         return dt
     else
         new_args = similar(args)
         for (i, arg) in enumerate(args)
             ex_arg = substitute_sample_time(arg, dt)
             if ex_arg isa NotInferedTimeDomain
                 return ex_arg
             end
             new_args[i] = ex_arg
@@ -250,7 +250,7 @@
         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,

src/systems/diffeqs/abstractodesystem.jl

@@ -1092,8 +1092,9 @@
         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,

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))

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

test/clock.jl

@@ -499,10 +499,10 @@
 @mtkmodel FirstOrderWithStepCounter begin
     @components begin
         counter = CounterSys()
         fo = FirstOrderSys()
     end
     @equations begin
         counter.u ~ fo.x
     end
 end
 
@@ -528,3 +528,50 @@
 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; phase = 100eps())) # TODO: this should not be required
+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) + 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.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)
+
+if length(d1) < length(d2)
+    d1, d2 = d2, d1
+end
+
+@test d2 == 1:6 # y
+
+# 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
+end
+X = manualtest()