Add split_system that splits the system by their time domain

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

Author: YingboMa

src/systems/abstractsystem.jl

@@ -1031,15 +1031,20 @@ This will convert all `inputs` to parameters and allow them to be unconnected, i
 simplification will allow models where `n_states = n_equations - n_inputs`.
 """
 function structural_simplify(sys::AbstractSystem, io = nothing; simplify = false,
-                             simplify_constants = true, kwargs...)
+                             kwargs...)
     sys = expand_connections(sys)
     sys isa DiscreteSystem && return sys
     state = TearingState(sys)
+    structural_simplify!(state, io; simplify, kwargs...)
+end
+
+function structural_simplify!(state::TearingState, io = nothing; simplify = false,
+                              kwargs...)
     has_io = io !== nothing
     has_io && markio!(state, io...)
     state, input_idxs = inputs_to_parameters!(state, io)
     sys, ag = alias_elimination!(state; kwargs...)
-    check_consistency(state, ag)
+    #check_consistency(state, ag)
     sys = dummy_derivative(sys, state, ag; simplify)
     fullstates = [map(eq -> eq.lhs, observed(sys)); states(sys)]
     @set! sys.observed = topsort_equations(observed(sys), fullstates)
@@ -1148,8 +1153,8 @@ end
 
 function markio!(state, inputs, outputs; check = true)
     fullvars = state.fullvars
-    inputset = Dict(inputs .=> false)
-    outputset = Dict(outputs .=> false)
+    inputset = Dict{Any, Bool}(i => false for i in inputs)
+    outputset = Dict{Any, Bool}(o => false for o in outputs)
     for (i, v) in enumerate(fullvars)
         if v in keys(inputset)
             v = setio(v, true, false)

src/systems/clock_inference.jl

@@ -28,7 +28,7 @@ function infer_clocks!(ci::ClockInference)
     @unpack ts, eq_domain, var_domain, inferred = ci
     @unpack fullvars = ts
     @unpack graph = ts.structure
-    # TODO: add a graph time to do this lazily
+    # TODO: add a graph type to do this lazily
     var_graph = SimpleGraph(ndsts(graph))
     for eq in 𝑠vertices(graph)
         vvs = 𝑠neighbors(graph, eq)
@@ -64,3 +64,78 @@ function infer_clocks!(ci::ClockInference)
 
     return ci
 end
+
+function resize_or_push!(v, val, idx)
+    n = length(v)
+    if idx > n
+        for i in (n + 1):idx
+            push!(v, Int[])
+        end
+        resize!(v, idx)
+    end
+    push!(v[idx], val)
+end
+
+function split_system(ci::ClockInference)
+    @unpack ts, eq_domain, var_domain, inferred = ci
+    @unpack fullvars = ts
+    @unpack graph = ts.structure
+    continuous_id = 0
+    clock_to_id = Dict{TimeDomain, Int}()
+    id_to_clock = TimeDomain[]
+    eq_to_cid = Vector{Int}(undef, nsrcs(graph))
+    cid_to_eq = Vector{Int}[]
+    var_to_cid = Vector{Int}(undef, ndsts(graph))
+    cid_to_var = Vector{Int}[]
+    cid = 0
+    for (i, d) in enumerate(eq_domain)
+        cid = get!(clock_to_id, d) do
+            cid += 1
+            push!(id_to_clock, d)
+            if d isa Continuous
+                continuous_id = cid
+            end
+            cid
+        end
+        eq_to_cid[i] = cid
+        resize_or_push!(cid_to_eq, i, cid)
+    end
+    input_discrete = Int[]
+    inputs = []
+    for (i, d) in enumerate(var_domain)
+        cid = get(clock_to_id, d, 0)
+        @assert cid!==0 "Internal error!"
+        var_to_cid[i] = cid
+        v = fullvars[i]
+        #TODO: remove Inferred*
+        if cid == continuous_id && istree(v) && (o = operation(v)) isa Operator &&
+           !(input_timedomain(o) isa Continuous)
+            push!(input_discrete, i)
+            push!(inputs, fullvars[i])
+        end
+        resize_or_push!(cid_to_var, i, cid)
+    end
+
+    eqs = equations(ts)
+    tss = similar(cid_to_eq, TearingState)
+    for (id, ieqs) in enumerate(cid_to_eq)
+        vars = cid_to_var[id]
+        ts_i = ts
+        fadj = Vector{Int}[]
+        eqs_i = Equation[]
+        var_set_i = BitSet(vars)
+        ne = 0
+        for eq_i in ieqs
+            vars = copy(graph.fadjlist[eq_i])
+            ne += length(vars)
+            push!(fadj, vars)
+            push!(eqs_i, eqs[eq_i])
+        end
+        @set! ts_i.structure.graph = complete(BipartiteGraph(ne, fadj, ndsts(graph)))
+        @set! ts_i.sys.eqs = eqs_i
+        tss[id] = ts_i
+    end
+    return tss, (; inputs, outputs = ())
+
+    #id_to_clock, cid_to_eq, cid_to_var
+end

test/clock.jl

@@ -21,9 +21,51 @@ eqs = [yd ~ Sample(t, dt)(y)
        y ~ x]
 @named sys = ODESystem(eqs)
 # compute equation and variables' time domains
+#TODO: test linearize
+
+#=
+ Differential(t)(x(t)) ~ u(t) - x(t)
+ 0 ~ Sample(Clock(t, 0.1))(y(t)) - yd(t)
+ 0 ~ kp*(r(t) - yd(t)) - ud(t)
+ 0 ~ Hold()(ud(t)) - u(t)
+ 0 ~ x(t) - y(t)
+
+====
+By inference:
+
+ Differential(t)(x(t)) ~ u(t) - x(t)
+ 0 ~ Hold()(ud(t)) - u(t) # Hold()(ud(t)) is constant except in an event
+ 0 ~ x(t) - y(t)
+
+ 0 ~ Sample(Clock(t, 0.1))(y(t)) - yd(t)
+ 0 ~ kp*(r(t) - yd(t)) - ud(t)
+
+====
+
+ Differential(t)(x(t)) ~ u(t) - x(t)
+ 0 ~ Hold()(ud(t)) - u(t)
+ 0 ~ x(t) - y(t)
+
+ yd(t) := Sample(Clock(t, 0.1))(y(t))
+ ud(t) := kp*(r(t) - yd(t))
+=#
+
+#=
+     D(x) ~ Shift(x, 0, dt) + 1 # this should never meet with continous variables
+=>   (Shift(x, 0, dt) - Shift(x, -1, dt))/dt ~ Shift(x, 0, dt) + 1
+=>   Shift(x, 0, dt) - Shift(x, -1, dt) ~ Shift(x, 0, dt) * dt + dt
+=>   Shift(x, 0, dt) - Shift(x, 0, dt) * dt ~ Shift(x, -1, dt) + dt
+=>   (1 - dt) * Shift(x, 0, dt) ~ Shift(x, -1, dt) + dt
+=>   Shift(x, 0, dt) := (Shift(x, -1, dt) + dt) / (1 - dt) # Discrete system
+=#
 
 ci, varmap = infer_clocks(sys)
 eqmap = ci.eq_domain
+tss, io = ModelingToolkit.split_system(deepcopy(ci))
+ts_c = deepcopy(tss[1])
+@set! ts_c.structure.solvable_graph = nothing
+sss, = ModelingToolkit.structural_simplify!(ts_c, io)
+@test equations(sss) == [D(x) ~ u - x]
 
 d = Clock(t, dt)
 # Note that TearingState reorders the equations