Merge pull request #2690 from SciML/inferred_sampletime

add `sampletime` operator

Author: YingboMa

src/ModelingToolkit.jl

@@ -269,7 +269,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
+export Sample, Hold, Shift, ShiftIndex, sampletime, SampleTime
 export Clock #, InferredDiscrete,
 
 end # module

src/discretedomain.jl

@@ -1,5 +1,11 @@
 using Symbolics: Operator, Num, Term, value, recursive_hasoperator
 
+struct SampleTime <: Operator end
+SymbolicUtils.promote_symtype(::Type{SampleTime}, t...) = Real
+function SampleTime()
+    SymbolicUtils.term(SampleTime, type = Real)
+end
+
 # Shift
 
 """
@@ -15,8 +21,6 @@ $(FIELDS)
 ```jldoctest
 julia> using Symbolics
 
-julia> @variables t;
-
 julia> Δ = Shift(t)
 (::Shift) (generic function with 2 methods)
 ```
@@ -176,16 +180,18 @@ end
 function (xn::Num)(k::ShiftIndex)
     @unpack clock, steps = k
     x = value(xn)
-    t = clock.t
     # Verify that the independent variables of k and x match and that the expression doesn't have multiple variables
     vars = Symbolics.get_variables(x)
     length(vars) == 1 ||
         error("Cannot shift a multivariate expression $x. Either create a new unknown and shift this, or shift the individual variables in the expression.")
     args = Symbolics.arguments(vars[]) # args should be one element vector with the t in x(t)
     length(args) == 1 ||
         error("Cannot shift an expression with multiple independent variables $x.")
-    isequal(args[], t) ||
-        error("Independent variable of $xn is not the same as that of the ShiftIndex $(k.t)")
+    t = args[]
+    if hasfield(typeof(clock), :t)
+        isequal(t, clock.t) ||
+            error("Independent variable of $xn is not the same as that of the ShiftIndex $(k.t)")
+    end
 
     # d, _ = propagate_time_domain(xn)
     # if d != clock # this is only required if the variable has another clock

src/systems/clock_inference.jl

@@ -21,6 +21,52 @@ function ClockInference(ts::TransformationState)
     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)
+    @unpack ts, eq_domain = ci
+    eqs = copy(equations(ts))
+    @assert length(eqs) == length(eq_domain)
+    for i in eachindex(eqs)
+        eq = eqs[i]
+        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
+    end
+    @set! ts.sys.eqs = eqs
+    @set! ci.ts = ts
+end
+
+function substitute_sample_time(eq::Equation, dt)
+    substitute_sample_time(eq.lhs, dt) ~ substitute_sample_time(eq.rhs, dt)
+end
+
+function substitute_sample_time(ex, dt)
+    istree(ex) || return ex
+    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
+        end
+        similarterm(ex, op, new_args; metadata = metadata(ex))
+    end
+end
+
 function infer_clocks!(ci::ClockInference)
     @unpack ts, eq_domain, var_domain, inferred = ci
     @unpack var_to_diff, graph = ts.structure
@@ -66,6 +112,7 @@ function infer_clocks!(ci::ClockInference)
         end
     end
 
+    ci = substitute_sample_time(ci)
     return ci
 end
 

src/systems/systemstructure.jl

@@ -627,7 +627,7 @@ function structural_simplify!(state::TearingState, io = nothing; simplify = fals
         kwargs...)
     if state.sys isa ODESystem
         ci = ModelingToolkit.ClockInference(state)
-        ModelingToolkit.infer_clocks!(ci)
+        ci = ModelingToolkit.infer_clocks!(ci)
         time_domains = merge(Dict(state.fullvars .=> ci.var_domain),
             Dict(default_toterm.(state.fullvars) .=> ci.var_domain))
         tss, inputs, continuous_id, id_to_clock = ModelingToolkit.split_system(ci)

test/clock.jl

@@ -330,7 +330,7 @@ using ModelingToolkitStandardLibrary.Blocks
 
 dt = 0.05
 d = Clock(t, dt)
-k = ShiftIndex(d)
+k = ShiftIndex()
 
 @mtkmodel DiscretePI begin
     @components begin
@@ -347,7 +347,7 @@ k = ShiftIndex(d)
         y(t)
     end
     @equations begin
-        x(k) ~ x(k - 1) + ki * u(k)
+        x(k) ~ x(k - 1) + ki * u(k) * SampleTime() / dt
         output.u(k) ~ y(k)
         input.u(k) ~ u(k)
         y(k) ~ x(k - 1) + kp * u(k)
@@ -364,21 +364,18 @@ end
     end
 end
 
-@mtkmodel Holder begin
-    @components begin
-        input = RealInput()
-        output = RealOutput()
-    end
+@mtkmodel ZeroOrderHold begin
+    @extend u, y = siso = Blocks.SISO()
     @equations begin
-        output.u ~ Hold(input.u)
+        y ~ Hold(u)
     end
 end
 
 @mtkmodel ClosedLoop begin
     @components begin
         plant = FirstOrder(k = 0.3, T = 1)
         sampler = Sampler()
-        holder = Holder()
+        holder = ZeroOrderHold()
         controller = DiscretePI(kp = 2, ki = 2)
         feedback = Feedback()
         ref = Constant(k = 0.5)
@@ -444,7 +441,7 @@ prob = ODEProblem(ssys,
     [model.plant.x => 0.0; model.controller.kp => 2.0; model.controller.ki => 2.0],
     (0.0, Tf))
 int = init(prob, Tsit5(); kwargshandle = KeywordArgSilent)
-@test int.ps[Hold(ssys.holder.input.u)] == 2 # constant output * kp issue https://github.com/SciML/ModelingToolkit.jl/issues/2356
+@test_broken int.ps[Hold(ssys.holder.input.u)] == 2 # constant output * kp issue https://github.com/SciML/ModelingToolkit.jl/issues/2356
 @test int.ps[ssys.controller.x] == 1 # c2d
 @test int.ps[Sample(d)(ssys.sampler.input.u)] == 0 # disc state
 sol = solve(prob,