Merge pull request #935 from SciML/myb/ss

Refactor & fix some minor bugs in alias elimination & add consistency check

Author: YingboMa

test/rc_model.jl renamed to examples/electrical_components.jl

@@ -8,13 +8,13 @@ function Pin(;name)
     ODESystem(Equation[], t, [v, i], [], name=name, defaults=[v=>1.0, i=>1.0])
 end
 
-function Ground(name)
+function Ground(;name)
     @named g = Pin()
     eqs = [g.v ~ 0]
     ODESystem(eqs, t, [], [], systems=[g], name=name)
 end
 
-function ConstantVoltage(name; V = 1.0)
+function ConstantVoltage(;name, V = 1.0)
     val = V
     @named p = Pin()
     @named n = Pin()
@@ -26,7 +26,7 @@ function ConstantVoltage(name; V = 1.0)
     ODESystem(eqs, t, [], [V], systems=[p, n], defaults=Dict(V => val), name=name)
 end
 
-function Resistor(name; R = 1.0)
+function Resistor(;name, R = 1.0)
     val = R
     @named p = Pin()
     @named n = Pin()
@@ -40,7 +40,7 @@ function Resistor(name; R = 1.0)
     ODESystem(eqs, t, [v], [R], systems=[p, n], defaults=Dict(R => val), name=name)
 end
 
-function Capacitor(name; C = 1.0)
+function Capacitor(;name, C = 1.0)
     val = C
     @named p = Pin()
     @named n = Pin()
@@ -55,15 +55,23 @@ function Capacitor(name; C = 1.0)
     ODESystem(eqs, t, [v], [C], systems=[p, n], defaults=Dict(C => val), name=name)
 end
 
-R = 1.0
-C = 1.0
-V = 1.0
-resistor = Resistor(:resistor, R=R)
-capacitor = Capacitor(:capacitor, C=C)
-source = ConstantVoltage(:source, V=V)
-ground = Ground(:ground)
+function Inductor(; name, L = 1.0)
+    val = L
+    @named p = Pin()
+    @named n = Pin()
+    @variables v(t) i(t)
+    @parameters L
+    D = Differential(t)
+    eqs = [
+           v ~ p.v - n.v
+           0 ~ p.i + n.i
+           i ~ p.i
+           D(i) ~ v / L
+          ]
+    ODESystem(eqs, t, [v, i], [L], systems=[p, n], defaults=Dict(L => val), name=name)
+end
 
-function connect(ps...)
+function connect_pins(ps...)
     eqs = [
            0 ~ sum(p->p.i, ps) # KCL
           ]
@@ -74,10 +82,3 @@ function connect(ps...)
 
     return eqs
 end
-rc_eqs = [
-          connect(source.p, resistor.p)
-          connect(resistor.n, capacitor.p)
-          connect(capacitor.n, source.n, ground.g)
-         ]
-
-rc_model = ODESystem(rc_eqs, t, systems=[resistor, capacitor, source, ground], name=:rc)

examples/rc_model.jl

@@ -0,0 +1,17 @@
+include("electrical_components.jl")
+
+R = 1.0
+C = 1.0
+V = 1.0
+@named resistor = Resistor(R=R)
+@named capacitor = Capacitor(C=C)
+@named source = ConstantVoltage(V=V)
+@named ground = Ground()
+
+rc_eqs = [
+          connect_pins(source.p, resistor.p)
+          connect_pins(resistor.n, capacitor.p)
+          connect_pins(capacitor.n, source.n, ground.g)
+         ]
+
+@named rc_model = ODESystem(rc_eqs, t, systems=[resistor, capacitor, source, ground])

examples/serial_inductor.jl

@@ -0,0 +1,16 @@
+include("electrical_components.jl")
+
+@named source = ConstantVoltage(V=10.0)
+@named resistor = Resistor(R=1.0)
+@named inductor1 = Inductor(L=1.0e-2)
+@named inductor2 = Inductor(L=2.0e-2)
+@named ground = Ground()
+
+eqs = [
+       connect_pins(source.p, resistor.p)
+       connect_pins(resistor.n, inductor1.p)
+       connect_pins(inductor1.n, inductor2.p)
+       connect_pins(source.n, inductor2.n, ground.g)
+      ]
+
+@named ll_model = ODESystem(eqs, t, systems=[source, resistor, inductor1, inductor2, ground])

src/bipartite_graph.jl

@@ -66,13 +66,13 @@ badjlist = [[1,2,5,6],[3,4,6]]
 bg = BipartiteGraph(7, fadjlist, badjlist)
 ```
 """
-mutable struct BipartiteGraph{I<:Integer,M} <: LightGraphs.AbstractGraph{I}
+mutable struct BipartiteGraph{I<:Integer,F<:Vector{Vector{I}},B<:Union{Vector{Vector{I}},I},M} <: LightGraphs.AbstractGraph{I}
     ne::Int
-    fadjlist::Vector{Vector{I}} # `fadjlist[src] => dsts`
-    badjlist::Vector{Vector{I}} # `badjlist[dst] => srcs`
+    fadjlist::F # `fadjlist[src] => dsts`
+    badjlist::B # `badjlist[dst] => srcs` or `ndsts`
     metadata::M
 end
-BipartiteGraph(ne::Integer, fadj::AbstractVector, badj::AbstractVector) = BipartiteGraph(ne, fadj, badj, nothing)
+BipartiteGraph(ne::Integer, fadj::AbstractVector, badj::Union{AbstractVector,Integer}=maximum(maximum, fadj); metadata=nothing) = BipartiteGraph(ne, fadj, badj, metadata)
 
 """
 ```julia
@@ -93,16 +93,16 @@ $(SIGNATURES)
 
 Build an empty `BipartiteGraph` with `nsrcs` sources and `ndsts` destinations.
 """
-function BipartiteGraph(nsrcs::T, ndsts::T; metadata=nothing) where T
+function BipartiteGraph(nsrcs::T, ndsts::T, backedge::Val{B}=Val(true); metadata=nothing) where {T,B}
     fadjlist = map(_->T[], 1:nsrcs)
-    badjlist = map(_->T[], 1:ndsts)
+    badjlist = B ? map(_->T[], 1:ndsts) : ndsts
     BipartiteGraph(0, fadjlist, badjlist, metadata)
 end
 
 Base.eltype(::Type{<:BipartiteGraph{I}}) where I = I
 function Base.empty!(g::BipartiteGraph)
     foreach(empty!, g.fadjlist)
-    foreach(empty!, g.badjlist)
+    g.badjlist isa AbstractVector && foreach(empty!, g.badjlist)
     g.ne = 0
     if g.metadata !== nothing
         foreach(empty!, g.metadata)
@@ -111,17 +111,22 @@ function Base.empty!(g::BipartiteGraph)
 end
 Base.length(::BipartiteGraph) = error("length is not well defined! Use `ne` or `nv`.")
 
+@noinline throw_no_back_edges() = throw(ArgumentError("The graph has no back edges."))
+
 if isdefined(LightGraphs, :has_contiguous_vertices)
     LightGraphs.has_contiguous_vertices(::Type{<:BipartiteGraph}) = false
 end
 LightGraphs.is_directed(::Type{<:BipartiteGraph}) = false
 LightGraphs.vertices(g::BipartiteGraph) = (𝑠vertices(g), 𝑑vertices(g))
 𝑠vertices(g::BipartiteGraph) = axes(g.fadjlist, 1)
-𝑑vertices(g::BipartiteGraph) = axes(g.badjlist, 1)
+𝑑vertices(g::BipartiteGraph) = g.badjlist isa AbstractVector ? axes(g.badjlist, 1) : Base.OneTo(g.badjlist)
 has_𝑠vertex(g::BipartiteGraph, v::Integer) = v in 𝑠vertices(g)
 has_𝑑vertex(g::BipartiteGraph, v::Integer) = v in 𝑑vertices(g)
 𝑠neighbors(g::BipartiteGraph, i::Integer, with_metadata::Val{M}=Val(false)) where M = M ? zip(g.fadjlist[i], g.metadata[i]) : g.fadjlist[i]
-𝑑neighbors(g::BipartiteGraph, j::Integer, with_metadata::Val{M}=Val(false)) where M = M ? zip(g.badjlist[j], (g.metadata[i][j] for i in g.badjlist[j])) : g.badjlist[j]
+function 𝑑neighbors(g::BipartiteGraph, j::Integer, with_metadata::Val{M}=Val(false)) where M
+    g.badjlist isa AbstractVector || throw_no_back_edges()
+    M ? zip(g.badjlist[j], (g.metadata[i][j] for i in g.badjlist[j])) : g.badjlist[j]
+end
 LightGraphs.ne(g::BipartiteGraph) = g.ne
 LightGraphs.nv(g::BipartiteGraph) = sum(length, vertices(g))
 LightGraphs.edgetype(g::BipartiteGraph{I}) where I = BipartiteEdge{I}
@@ -145,7 +150,6 @@ const NO_METADATA = NoMetadata()
 LightGraphs.add_edge!(g::BipartiteGraph, i::Integer, j::Integer, md=NO_METADATA) = add_edge!(g, BipartiteEdge(i, j), md)
 function LightGraphs.add_edge!(g::BipartiteGraph, edge::BipartiteEdge, md=NO_METADATA)
     @unpack fadjlist, badjlist = g
-    verts = vertices(g)
     s, d = src(edge), dst(edge)
     (has_𝑠vertex(g, s) && has_𝑑vertex(g, d)) || error("edge ($edge) out of range.")
     @inbounds list = fadjlist[s]
@@ -157,15 +161,21 @@ function LightGraphs.add_edge!(g::BipartiteGraph, edge::BipartiteEdge, md=NO_MET
     end
 
     g.ne += 1
-    @inbounds list = badjlist[d]
-    index = searchsortedfirst(list, s)
-    insert!(list, index, s)
+    if badjlist isa AbstractVector
+        @inbounds list = badjlist[d]
+        index = searchsortedfirst(list, s)
+        insert!(list, index, s)
+    end
     return true  # edge successfully added
 end
 
 function LightGraphs.add_vertex!(g::BipartiteGraph{T}, type::VertType) where T
     if type === DST
-        push!(g.badjlist, T[])
+        if g.badjlist isa AbstractVector
+            push!(g.badjlist, T[])
+        else
+            g.badjlist += 1
+        end
     elseif type === SRC
         push!(g.fadjlist, T[])
     else

src/structural_transformation/StructuralTransformations.jl

@@ -16,7 +16,7 @@ using ModelingToolkit: ODESystem, var_from_nested_derivative, Differential,
                        states, equations, vars, Symbolic, diff2term, value,
                        operation, arguments, Sym, Term, simplify, solve_for,
                        isdiffeq, isdifferential,
-                       get_structure, get_reduced_states, defaults
+                       get_structure, defaults, InvalidSystemException
 
 using ModelingToolkit.BipartiteGraphs
 using ModelingToolkit.SystemStructures
@@ -31,7 +31,7 @@ using SparseArrays
 
 using NonlinearSolve
 
-export tearing, dae_index_lowering
+export tearing, dae_index_lowering, check_consistency
 export build_torn_function, build_observed_function, ODAEProblem
 export sorted_incidence_matrix
 

src/structural_transformation/codegen.jl

@@ -40,26 +40,27 @@ function torn_system_jacobian_sparsity(sys)
     # dependencies.
     avars2dvars = Dict{Int,Set{Int}}()
     c = 0
-    for (_, _, teqs, tvars) in partitions
+    for partition in partitions
+        @unpack e_residual, v_residual = partition
         # initialization
-        for tvar in tvars
+        for tvar in v_residual
             avars2dvars[tvar] = Set{Int}()
         end
-        for teq in teqs
+        for teq in e_residual
             c += 1
             for var in 𝑠neighbors(graph, teq)
                 # Skip the tearing variables in the current partition, because
                 # we are computing them from all the other states.
-                LightGraphs.insorted(var, tvars) && continue
+                LightGraphs.insorted(var, v_residual) && continue
                 deps = get(avars2dvars, var, nothing)
                 if deps === nothing # differential variable
                     @assert !isalgvar(s, var)
-                    for tvar in tvars
+                    for tvar in v_residual
                         push!(avars2dvars[tvar], var)
                     end
                 else # tearing variable from previous partitions
                     @assert isalgvar(s, var)
-                    for tvar in tvars
+                    for tvar in v_residual
                         union!(avars2dvars[tvar], avars2dvars[var])
                     end
                 end
@@ -97,22 +98,22 @@ function partitions_dag(s::SystemStructure)
     @unpack partitions, graph = s
 
     # `partvars[i]` contains all the states that appear in `partitions[i]`
-    partvars = map(partitions) do (_, _, reqs, tvars)
+    partvars = map(partitions) do partition
         ipartvars = Set{Int}()
-        for req in reqs
+        for req in partition.e_residual
             union!(ipartvars, 𝑠neighbors(graph, req))
         end
         ipartvars
     end
 
     I, J = Int[], Int[]
     n = length(partitions)
-    for i in 1:n
+    for (i, partition) in enumerate(partitions)
         for j in i+1:n
             # The only way for a later partition `j` to depend on an earlier
             # partition `i` is when `partvars[j]` contains one of tearing
             # variables of partition `i`.
-            if !isdisjoint(partvars[j], partitions[i][4])
+            if !isdisjoint(partvars[j], partition.v_residual)
                 # j depends on i
                 push!(I, i)
                 push!(J, j)
@@ -170,8 +171,8 @@ function get_torn_eqs_vars(sys)
     vars = s.fullvars
     eqs = equations(sys)
 
-    torn_eqs  = map(idxs-> eqs[idxs], map(x->x[3], partitions))
-    torn_vars = map(idxs->vars[idxs], map(x->x[4], partitions))
+    torn_eqs  = map(idxs-> eqs[idxs], map(x->x.e_residual, partitions))
+    torn_vars = map(idxs->vars[idxs], map(x->x.v_residual, partitions))
 
     gen_nlsolve.((sys,), torn_eqs, torn_vars)
 end
@@ -197,7 +198,7 @@ function build_torn_function(
     )
 
     s = structure(sys)
-    states = s.fullvars[diffvars_range(s)]
+    states = map(i->s.fullvars[i], diffvars_range(s))
     syms = map(Symbol, states)
 
     expr = SymbolicUtils.Code.toexpr(
@@ -243,9 +244,9 @@ given a set of `vars`, find the groups of equations we need to solve for
 to obtain the solution to `vars`
 """
 function find_solve_sequence(partitions, vars)
-    subset = filter(x -> !isdisjoint(x[4], vars), partitions)
+    subset = filter(x -> !isdisjoint(x.v_residual, vars), partitions)
     isempty(subset) && return []
-    vars′ = mapreduce(x->x[4], union, subset)
+    vars′ = mapreduce(x->x.v_residual, union, subset)
     if vars′ == vars
         return subset
     else
@@ -267,8 +268,8 @@ function build_observed_function(
     syms_set = Set(syms)
     s = structure(sys)
     @unpack partitions, fullvars, graph = s
-    diffvars = fullvars[diffvars_range(s)]
-    algvars = fullvars[algvars_range(s)]
+    diffvars = map(i->fullvars[i], diffvars_range(s))
+    algvars = map(i->fullvars[i], algvars_range(s))
 
     required_algvars = Set(intersect(algvars, syms_set))
     obs = observed(sys)
@@ -290,8 +291,8 @@ function build_observed_function(
     if !isempty(subset)
         eqs = equations(sys)
 
-        torn_eqs  = map(idxs-> eqs[idxs[3]], subset)
-        torn_vars = map(idxs->fullvars[idxs[4]], subset)
+        torn_eqs  = map(idxs-> eqs[idxs.e_residual], subset)
+        torn_vars = map(idxs->fullvars[idxs.v_residual], subset)
 
         solves = gen_nlsolve.((sys,), torn_eqs, torn_vars; checkbounds=checkbounds)
     else
@@ -338,7 +339,7 @@ function ODAEProblem{iip}(
                          ) where {iip}
     s = structure(sys)
     @unpack fullvars = s
-    dvs = fullvars[diffvars_range(s)]
+    dvs = map(i->fullvars[i], diffvars_range(s))
     ps = parameters(sys)
     defs = defaults(sys)