ModelingToolkit experiments with functional composition of models

Here's an attempt to functionally compose MTK models where each model is an ODESystem:

using ModelingToolkit, Symbolics


struct NodeSet end

isnode(x::Num) = isnode(Symbolics.value(x))
isnode(x::Symbolics.Symbolic) = hasmetadata(x, NodeSet)
isnode(x) = false

"""
    tonode(s::Sym)
Maps the variable to a node variable.
"""
tonode(s::Symbolics.Symbolic) = setmetadata(s, NodeSet, Set())
tonode(s::Num) = Num(tonode(Symbolics.value(s)))

"""
Define one or more known node variables.
"""
macro nodes(xs...)
    Symbolics._parse_vars(:nodes,
                          Real,
                          xs,
                          x -> x isa Array ? tonode.(x) : tonode(x)
                         ) |> esc
end

RefBranch(x, i) = nothing
function RefBranch(n::Num, i)
    push!(getmetadata(n.val, NodeSet), i)
    nothing
end

function Branch(n1, n2, v, i)
    RefBranch(n1, i)
    RefBranch(n2, -i)
    v ~ n1 - n2
end

function addflows!(sys::ODESystem)
    for s in states(sys)
        if isnode(s)
            push!(sys.eqs, 0 ~ sum(getmetadata(s, NodeSet)))
        end
    end
end

@parameters t
const D = Differential(t)

function VoltageSource(n1, n2, V; name) 
    @variables i(t)
    @variables v(t)
    eqs = [
        Branch(n1, n2, v, i)
        v ~ V
    ]
    ODESystem(eqs, t, [i, v], [], name=name)
end

function Resistor(n1, n2, R; name) 
    @variables i(t)
    @variables v(t)
    eqs = [
        Branch(n1, n2, v, i)
        v ~ R * i
    ]
    ODESystem(eqs, t, [i, v], [], name=name)
end

function Capacitor(n1, n2, C; name) 
    @variables i(t)
    @variables v(t)
    eqs = [
        Branch(n1, n2, v, i)
        D(v) ~ i / C
    ]
    ODESystem(eqs, t, [i, v], [], name=name)
end

function Subsystem(n1, n2; name)
    @nodes vs(t)
    g = 0.0  # A ground has zero volts; it's not a variable.
    systems = [
        Resistor(n1, vs, 10.0, name = :r1)
        Capacitor(vs, g, 5.0e-3, name = :c1)
        Resistor(vs, n2, 10.0, name = :r2)
    ]
    ODESystem(Equation[], t, [vs], [], systems=systems, name=name)
end

function Circuit(; name)
    @nodes v1(t) v2(t)
    g = 0.0  # A ground has zero volts; it's not a variable.
    systems = [
        VoltageSource(v1, g, sin(2pi * 60 * t), name = :vsrc)
        Subsystem(v1, v2, name = :ss)
        Resistor(v2, g, 5.0, name = :r)
    ]
    ODESystem(Equation[], t, [v1, v2], [], systems=systems, name=name)
end

@named ckt = Circuit()

addflows!(ckt)

The main issues with this are:

Some of the variable names in the equations do not have the right levels (SciML/ModelingToolkit.jl#971).
The equations for the sum of the flow variables are not right. Multiple variables with the name i(t) (but different instances) are canceling each other.

tshort/summary.md