Incorrect observed variables after structural_simplify

[vvainola]

I have used the RC-circuit in the examples for simple 3-phase system that has only a voltage source and inductor in series in each phase. Depending on the order of components in "compose" after creating connections, the observed variables are different

using ModelingToolkit, DifferentialEquations

@variables t
@connector function Pin(; name)
    sts = @variables v(t) = 1.0 i(t) = 1.0
    ODESystem(Equation[], t, sts, []; name = name)
end

function ModelingToolkit.connect(::Type{Pin}, ps...)
    eqs = [
        0 ~ sum(p -> p.i, ps), # KCL
    ]
    # KVL
    for i = 1:length(ps)-1
        push!(eqs, ps[i].v ~ ps[i+1].v)
    end

    return eqs
end

function Ground(; name)
    @named g = Pin()
    eqs = [g.v ~ 0]
    compose(ODESystem(eqs, t, [], []; name = name), g)
end

function OnePort(; name)
    @named p = Pin()
    @named n = Pin()
    sts = @variables v(t) = 1.0 i(t) = 1.0
    eqs = [
        v ~ p.v - n.v
        0 ~ p.i + n.i
        i ~ p.i
    ]
    compose(ODESystem(eqs, t, sts, []; name = name), p, n)
end

function Inductor(; name, L = 1.0)
    @named oneport = OnePort()
    @unpack v, i = oneport
    ps = @parameters L = L
    D = Differential(t)
    eqs = [D(i) ~ v / L]
    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
end

function ConstantVoltage(; name, V = 1.0)
    @named oneport = OnePort()
    @unpack v = oneport
    ps = @parameters V = V
    eqs = [V ~ v]
    extend(ODESystem(eqs, t, [], ps; name = name), oneport)
end

R = 1.0
V = 1.0
L = 1.0

@named L_a = Inductor(L = L)
@named L_b = Inductor(L = L)
@named L_c = Inductor(L = L)

@named U_a = ConstantVoltage(V = V)
@named U_b = ConstantVoltage(V = V)
@named U_c = ConstantVoltage(V = V)

@named ground = Ground()

rl_eqs = [
    connect(U_a.n, U_b.n, U_c.n)
    connect(U_a.p, L_a.p)
    connect(U_b.p, L_b.p)
    connect(U_c.p, L_c.p)
    connect(L_a.n, L_b.n, L_c.n, ground.g)
]

@named _rl_model = ODESystem(rl_eqs, t)
@named rl_model = compose(
    _rl_model,
    [
        U_a,
        U_b,
        U_c,
        L_a,
        L_b,
        L_c,
        ground,
    ],
)

sys = structural_simplify(rl_model)
for eq in sys.eqs
    println(eq)
end
println()
for eq in sys.observed
    println(eq)
end

results in

Differential(t)(L_a₊i(t)) ~ (U_a₊V + L_b₊v(t) - U_b₊V) / L_a₊L
Differential(t)(L_b₊i(t)) ~ L_b₊v(t) / L_b₊L
Differential(t)(L_c₊i(t)) ~ (U_c₊V + L_b₊v(t) - U_b₊V) / L_c₊L
0 ~ (-(U_a₊V + L_b₊v(t) - U_b₊V)) / L_a₊L + (-(U_c₊V + L_b₊v(t) - U_b₊V)) / L_c₊L + (-L_b₊v(t)) / L_b₊L

U_c₊i(t) ~ -L_c₊i(t)
L_c₊n₊i(t) ~ -L_c₊i(t)
U_b₊i(t) ~ -L_b₊i(t)
L_b₊n₊i(t) ~ -L_b₊i(t)
L_b₊p₊v(t) ~ L_b₊v(t)
U_a₊i(t) ~ -L_a₊i(t)
L_a₊n₊i(t) ~ -L_a₊i(t)
U_c₊p₊i(t) ~ -L_c₊i(t)
U_c₊n₊i(t) ~ L_c₊i(t)
U_b₊p₊i(t) ~ -L_b₊i(t)
U_b₊n₊i(t) ~ L_b₊i(t)
U_a₊p₊i(t) ~ -L_a₊i(t)
U_a₊n₊i(t) ~ L_a₊i(t)
L_a₊p₊i(t) ~ L_a₊i(t)
L_c₊p₊i(t) ~ L_c₊i(t)
U_b₊p₊v(t) ~ L_b₊v(t)
L_b₊p₊i(t) ~ L_b₊i(t)
L_a₊n₊v(t) ~ 0
L_b₊n₊v(t) ~ 0
L_c₊n₊v(t) ~ 0
ground₊g₊v(t) ~ 0
ground₊g₊i(t) ~ L_a₊i(t) + L_b₊i(t) + L_c₊i(t)
U_c₊v(t) ~ U_c₊V
U_b₊v(t) ~ U_b₊V
U_a₊v(t) ~ U_a₊V
L_c₊v(t) ~ L_b₊v(t) + U_c₊v(t) - U_b₊v(t)
L_c₊p₊v(t) ~ L_c₊v(t)
U_c₊p₊v(t) ~ L_c₊v(t)
L_a₊v(t) ~ L_c₊v(t) + U_a₊v(t) - U_c₊v(t)
U_c₊n₊v(t) ~ L_c₊v(t) - U_c₊v(t)
L_a₊p₊v(t) ~ L_a₊v(t)
U_a₊p₊v(t) ~ L_a₊v(t)
U_a₊n₊v(t) ~ U_c₊n₊v(t)
U_b₊n₊v(t) ~ U_c₊n₊v(t)

which looks ok except that I think in sys.eqs there should be only 2 states since there is 3 inductors and the current of third is the sum of the 2 others. However, if

@named rl_model = compose(
    _rl_model,
    [
        U_a,
        U_b,
        U_c,
        L_a,
        L_b,
        L_c,
        ground,
    ],
)

is changed to

@named rl_model = compose(
    _rl_model,
    [
        L_a,
        L_b,
        L_c,
        U_a,
        U_b,
        U_c,
        ground,
    ],
)

the output is then

Differential(t)(L_a₊i(t)) ~ (U_a₊V + L_b₊v(t) - U_b₊V) / L_a₊L
Differential(t)(L_b₊i(t)) ~ L_b₊v(t) / L_b₊L
Differential(t)(L_c₊i(t)) ~ (U_c₊V + L_b₊v(t) - U_b₊V) / L_c₊L
0 ~ (-(U_a₊V + L_b₊v(t) - U_b₊V)) / L_a₊L + (-(U_c₊V + L_b₊v(t) - U_b₊V)) / L_c₊L + (-L_b₊v(t)) / L_b₊L

U_c₊i(t) ~ -L_c₊i(t)
U_c₊n₊i(t) ~ L_c₊i(t)
U_b₊i(t) ~ -L_b₊i(t)
U_b₊n₊i(t) ~ L_b₊i(t)
U_a₊i(t) ~ -L_a₊i(t)
U_a₊n₊i(t) ~ L_a₊i(t)
U_c₊p₊i(t) ~ -L_c₊i(t)
U_b₊p₊i(t) ~ -L_b₊i(t)
L_b₊p₊v(t) ~ L_b₊v(t)
U_a₊p₊i(t) ~ -L_a₊i(t)
L_a₊p₊i(t) ~ L_a₊i(t)
L_c₊p₊i(t) ~ L_c₊i(t)
U_b₊p₊v(t) ~ L_b₊v(t)
L_b₊p₊i(t) ~ L_b₊i(t)
L_a₊n₊v(t) ~ 0
L_b₊n₊v(t) ~ 0
L_c₊n₊v(t) ~ 0
ground₊g₊v(t) ~ 0
ground₊g₊i(t) ~ L_a₊i(t) + L_b₊i(t) + L_c₊i(t)
L_a₊n₊i(t) ~ (-3//1)*L_a₊i(t)
U_a₊v(t) ~ U_a₊V
U_b₊v(t) ~ U_b₊V
U_c₊v(t) ~ U_c₊V
L_b₊n₊i(t) ~ (-3//1)*L_b₊i(t)
L_c₊n₊i(t) ~ (-3//1)*L_c₊i(t)
L_a₊v(t) ~ L_b₊v(t) + U_a₊v(t) - U_b₊v(t)
L_a₊p₊v(t) ~ L_a₊v(t)
U_a₊p₊v(t) ~ L_a₊v(t)
L_c₊v(t) ~ L_a₊v(t) + U_c₊v(t) - U_a₊v(t)
U_c₊n₊v(t) ~ L_a₊v(t) - U_a₊v(t)
L_c₊p₊v(t) ~ L_c₊v(t)
U_c₊p₊v(t) ~ L_c₊v(t)
U_a₊n₊v(t) ~ U_c₊n₊v(t)
U_b₊n₊v(t) ~ U_c₊n₊v(t)

which has additional multiplications by 3 in some of the observed equations

L_a₊n₊i(t) ~ (-3//1)*L_a₊i(t)
L_a₊p₊i(t) ~ L_a₊i(t)

Julia 1.6.3
ModelingToolkit 7.1.1

[ChrisRackauckas]

Fixed awhile back.