Propagate input domain back to the variable and add multirate tests

Co-authored-by: Fredrik Bagge Carlson <[email protected]>

Author: YingboMa

src/clock.jl

@@ -89,6 +89,14 @@ See also [`has_continuous_domain`](@ref)
 """
 is_continuous_domain(x) = !has_discrete_domain(x) && has_continuous_domain(x)
 
+struct ClockInferenceException <: Exception
+    msg::Any
+end
+
+function Base.showerror(io::IO, cie::ClockInferenceException)
+    print(io, "ClockInferenceException: ", cie.msg)
+end
+
 abstract type AbstractClock <: AbstractDiscrete end
 
 """

src/systems/clock_inference.jl

@@ -46,8 +46,12 @@ function infer_clocks!(ci::ClockInference)
         c = BitSet(c′)
         idxs = intersect(c, inferred)
         isempty(idxs) && continue
+        for i in idxs
+            @show var_domain[i]
+        end
         if !allequal(var_domain[i] for i in idxs)
-            throw(ClockInferenceException("Clocks are not consistent in connected component $(fullvars[c])"))
+            display(fullvars[c′])
+            throw(ClockInferenceException("Clocks are not consistent in connected component $(fullvars[c′])"))
         end
         vd = var_domain[first(idxs)]
         for v in c′

src/systems/systemstructure.jl

@@ -236,13 +236,15 @@ function TearingState(sys; quick_cancel = false, check = true)
         isalgeq = true
         statevars = []
         for var in vars
+            set_incidence = true
+            @label ANOTHER_VAR
             _var, _ = var_from_nested_derivative(var)
             any(isequal(_var), ivs) && continue
             if isparameter(_var) || (istree(_var) && isparameter(operation(_var)))
                 continue
             end
             varidx = addvar!(var)
-            push!(statevars, var)
+            set_incidence && push!(statevars, var)
 
             dvar = var
             idx = varidx
@@ -254,6 +256,14 @@ function TearingState(sys; quick_cancel = false, check = true)
                 dvar = arguments(dvar)[1]
                 idx = addvar!(dvar)
             end
+
+            if istree(var) && operation(var) isa Symbolics.Operator &&
+               !isdifferential(var) && (it = input_timedomain(var)) !== nothing
+                set_incidence = false
+                var = only(arguments(var))
+                var = setmetadata(var, ModelingToolkit.TimeDomain, it)
+                @goto ANOTHER_VAR
+            end
         end
         push!(symbolic_incidence, copy(statevars))
         empty!(statevars)

test/clock.jl

@@ -1,4 +1,12 @@
 using ModelingToolkit, Test
+
+function infer_clocks(sys)
+    ts = TearingState(sys)
+    ci = ModelingToolkit.ClockInference(ts)
+    ModelingToolkit.infer_clocks!(ci), Dict(ci.ts.fullvars .=> ci.var_domain)
+end
+
+@info "Testing hybrid system"
 dt = 0.1
 @variables t x(t) y(t) u(t) yd(t) ud(t) r(t)
 @parameters kp
@@ -14,10 +22,7 @@ eqs = [yd ~ Sample(t, dt)(y)
 @named sys = ODESystem(eqs)
 # compute equation and variables' time domains
 
-ts = TearingState(sys)
-ci = ModelingToolkit.ClockInference(ts)
-ModelingToolkit.infer_clocks!(ci)
-varmap = Dict(ci.ts.fullvars .=> ci.var_domain)
+ci, varmap = infer_clocks(sys)
 eqmap = ci.eq_domain
 
 d = Clock(t, dt)
@@ -34,3 +39,38 @@ d = Clock(t, dt)
 @test varmap[x] == Continuous()
 @test varmap[y] == Continuous()
 @test varmap[u] == Continuous()
+
+@info "Testing multi-rate hybrid system"
+dt = 0.1
+dt2 = 0.2
+@variables t x(t) y(t) u(t) r(t) yd1(t) ud1(t) yd2(t) ud2(t)
+@parameters kp
+D = Differential(t)
+
+eqs = [
+       # controller (time discrete part `dt=0.1`)
+       yd1 ~ Sample(t, dt)(y)
+       ud1 ~ kp * (Sample(t, dt)(r) - yd1)
+       yd2 ~ Sample(t, dt2)(y)
+       ud2 ~ kp * (Sample(t, dt2)(r) - yd2)
+
+       # plant (time continuous part)
+       u ~ Hold(ud1) + Hold(ud2)
+       D(x) ~ -x + u
+       y ~ x]
+@named sys = ODESystem(eqs)
+ci, varmap = infer_clocks(sys)
+
+d = Clock(t, dt)
+d2 = Clock(t, dt2)
+#@test get_eq_domain(eqs[1]) == d
+#@test get_eq_domain(eqs[3]) == d2
+
+@test varmap[yd1] == d
+@test varmap[ud1] == d
+@test varmap[yd2] == d2
+@test varmap[ud2] == d2
+@test varmap[r] == Continuous()
+@test varmap[x] == Continuous()
+@test varmap[y] == Continuous()
+@test varmap[u] == Continuous()