Merge pull request #1534 from SciML/myb/param

Add but_ordered_incidence

Author: YingboMa

examples/electrical_components.jl

@@ -66,3 +66,41 @@ function Inductor(; name, L = 1.0)
           ]
     extend(ODESystem(eqs, t, [], ps; name=name), oneport)
 end
+
+@connector function HeatPort(;name)
+    @variables T(t)=293.15 Q_flow(t)=0.0 [connect = Flow]
+    ODESystem(Equation[], t, [T, Q_flow], [], name=name)
+end
+
+function HeatingResistor(;name, R=1.0, TAmbient=293.15, alpha=1.0)
+    @named p = Pin()
+    @named n = Pin()
+    @named h = HeatPort()
+    @variables v(t) RTherm(t)
+    @parameters R=R TAmbient=TAmbient alpha=alpha
+    eqs = [
+           RTherm ~ R*(1 + alpha*(h.T - TAmbient))
+           v ~ p.i * RTherm
+           h.Q_flow ~ -v * p.i # -LossPower
+           v ~ p.v - n.v
+           0 ~ p.i + n.i
+          ]
+    compose(ODESystem(
+        eqs, t, [v, RTherm], [R, TAmbient, alpha],
+        name=name,
+    ), p, n, h)
+end
+
+function HeatCapacitor(;name, rho=8050, V=1, cp=460, TAmbient=293.15)
+    @parameters rho=rho V=V cp=cp
+    C = rho*V*cp
+    @named h = HeatPort()
+    D = Differential(t)
+    eqs = [
+           D(h.T) ~ h.Q_flow / C
+          ]
+    compose(ODESystem(
+        eqs, t, [], [rho, V, cp],
+        name=name,
+    ), h)
+end

src/structural_transformation/StructuralTransformations.jl

@@ -45,6 +45,7 @@ export sorted_incidence_matrix, pantelides!, tearing_reassemble, find_solvables!
 export tearing_assignments, tearing_substitution
 export torn_system_jacobian_sparsity
 export full_equations
+export but_ordered_incidence, lowest_order_variable_mask, highest_order_variable_mask
 
 include("utils.jl")
 include("pantelides.jl")

src/structural_transformation/utils.jl

@@ -208,6 +208,38 @@ function find_solvables!(state::TearingState; kwargs...)
     return nothing
 end
 
+highest_order_variable_mask(ts) = let v2d = ts.structure.var_to_diff
+    v->isempty(outneighbors(v2d, v))
+end
+
+lowest_order_variable_mask(ts) = let v2d = ts.structure.var_to_diff
+    v->isempty(outneighbors(v2d, v))
+end
+
+function but_ordered_incidence(ts::TearingState, varmask=highest_order_variable_mask(ts))
+    graph = complete(ts.structure.graph)
+    var_eq_matching = complete(maximal_matching(graph, _->true, varmask))
+    scc = find_var_sccs(graph, var_eq_matching)
+    vordering = Vector{Int}(undef, 0)
+    bb = Int[1]
+    sizehint!(vordering, ndsts(graph))
+    sizehint!(bb, ndsts(graph))
+    l = 1
+    for c in scc
+        isemptyc = true
+        for v in c
+            if varmask(v)
+                push!(vordering, v)
+                l += 1
+                isemptyc = false
+            end
+        end
+        isemptyc || push!(bb, l)
+    end
+    mm = incidence_matrix(graph)
+    mm[[var_eq_matching[v] for v in vordering if var_eq_matching[v] isa Int], vordering], bb
+end
+
 # debugging use
 function reordered_matrix(sys, torn_matching)
     s = TearingState(sys)

test/components.jl

@@ -212,3 +212,38 @@ end
 
 @named foo = Circuit()
 @test structural_simplify(foo) isa ModelingToolkit.AbstractSystem
+
+# BLT tests
+using LinearAlgebra
+function parallel_rc_model(i; name, source, ground, R, C)
+    resistor = HeatingResistor(name=Symbol(:resistor, i), R=R)
+    capacitor = Capacitor(name=Symbol(:capacitor, i), C=C)
+    heat_capacitor = HeatCapacitor(name=Symbol(:heat_capacitor, i))
+
+    rc_eqs = [
+              connect(source.p, resistor.p)
+              connect(resistor.n, capacitor.p)
+              connect(capacitor.n, source.n, ground.g)
+              connect(resistor.h, heat_capacitor.h)
+             ]
+
+    compose(ODESystem(rc_eqs, t, name=Symbol(name, i)),
+            [resistor, capacitor, source, ground, heat_capacitor])
+end
+V = 2.0
+@named source = ConstantVoltage(V=V)
+@named ground = Ground()
+N = 50
+Rs = 10 .^range(0, stop=-4, length=N)
+Cs = 10 .^range(-3, stop=0, length=N)
+rc_systems = map(1:N) do i
+    parallel_rc_model(i; name=:rc, source=source, ground=ground, R=Rs[i], C=Cs[i])
+end;
+@variables E(t)=0.0
+eqs = [
+       D(E) ~ sum(((i, sys),)->getproperty(sys, Symbol(:resistor, i)).h.Q_flow, enumerate(rc_systems))
+      ]
+@named _big_rc = ODESystem(eqs, t, [E], [])
+@named big_rc = compose(_big_rc, rc_systems)
+ts = TearingState(expand_connections(big_rc))
+@test istriu(but_ordered_incidence(ts)[1])