Get observed variable from solution error

[hzgzh]

fails to get the observed variable for "sol[sink.p]", the error code :"Cannot find the parent of var"pipe2₊p(t)[1]ˍt"
MWE

using ModelingToolkit, LinearAlgebra, DifferentialEquations
using Plots

@variables t
der = Differential(t)
function cheb(n)
    @assert n>0 "N should great than 1"
    
    x=cos.(π*(0:n)./n)
    c=vcat(2,ones(n-1),2).*(-1).^(0:n)
    X=repeat(x,1,n+1)
    ΔX=X-X'
    D=(c*(1 ./c)')./(ΔX+I)
    D-=diagm(sum(D,dims=2)[:])
    D,x
end

@connector function Port(;name, p = 101325, m_dot = 0.0)
    sts = @variables p(t) = p m_dot(t) = m_dot [connect = Flow]
    ODESystem(Equation[], t, sts, []; name=name)
end

function Source(;name, p)
    @named oup = Port()
    sts = @variables m_dot(t) = 0.0
    eqs = [
        oup.p ~ p
        m_dot ~ -oup.m_dot
        ]
    compose(ODESystem(eqs, t, sts, []; name=name), [oup])
end

function Sink(;name, m_dot)
    @named inp = Port()
    sts = @variables p(t) = 0.0
    eqs = [
        inp.m_dot ~ m_dot
        p ~ inp.p
        ]
    compose(ODESystem(eqs, t, sts, []; name=name), [inp])
end

function Pipe(; name, L, N, d, λ)
    @named inp = Port()
    @named oup = Port()
    # generate chebyshev matr
    D, x = cheb(N)
    x=(L/2)*x
    D=D/(L/2)
    
    Dp=D[1:end-1,:]
    Dg=D[2:end,:]
    R = 287.0 #gas constant
    T = 300.0 #air 300℃
    A = π*d^2.0/4.0
    #calculate pipe properties
    #ρ = 3e5/R/T

    @variables m_dot(t)[1:N + 1] p(t)[1:N + 1] ρ(t)
    sts = [m_dot..., p..., ρ]
    eqs = [
        ρ ~ (inp.p + oup.p) / (2R*T)
        p[end] ~ inp.p
        m_dot[1] ~ -oup.m_dot
        [der(p[i]) ~ -R*T/A * dot(Dg[i,:], [m_dot...]) for i in 1 : N]
        [der(m_dot[i+1]) ~ -A * dot(Dp[i,:], [p...]) - λ/(2A*d*ρ) * m_dot[i+1] * abs(m_dot[i]) for i in 1:N]
        inp.m_dot ~ m_dot[end]
        oup.p ~ p[1]
    ]
    compose(ODESystem(eqs, t, sts, [];name = name, defaults = [p[6] => 101325]), [inp, oup])
end

@named source = Source(p = 3e5)
@named sink = Sink(m_dot = 1.00)

@named pipe1 = Pipe(; N = 5, L = 30e3, λ = 0.1, d = 0.5)
@named pipe2 = Pipe(; N = 5, L = 30e3, λ = 0.1, d = 0.5)
@named pipe3 = Pipe(; N = 5, L = 30e3, λ = 0.1, d = 0.5)
@named pipe4 = Pipe(; N = 5, L = 30e3, λ = 0.1, d = 0.5)

subs = [source, pipe1, pipe2, pipe3, pipe4, sink]

eqs = [
        connect(source.oup, pipe1.inp)
        connect(pipe1.oup,  pipe2.inp, pipe3.inp)
        connect(pipe2.oup,  pipe3.oup, pipe4.inp)
        connect(pipe4.oup,  sink.inp)
    ]
    
@named system = ODESystem(eqs, t, systems = subs)
sys = structural_simplify(system)

tspan = (0.0, 2*3600.0)

u0 = [
    collect(pipe1.m_dot) .=> 1.13
    collect(pipe1.p) .=> 3e5
    collect(pipe2.m_dot) .=> 0.56
    collect(pipe2.p) .=> 3e5
    collect(pipe3.m_dot) .=> 0.56
    collect(pipe3.p) .=> 3e5
    collect(pipe4.m_dot) .=> 1.13
    collect(pipe4.p) .=> 3e5
]

prob = ODAEProblem(sys, u0, tspan)

prob1 = SteadyStateProblem(prob)
sol1 = solve(prob1, SSRootfind())

prob2 = remake(prob, u0 = sol1)
@time sol = solve(prob2, Rodas42())

after run "observed(sys)"
I get

sink₊p(t) ~ (pipe4₊p(t))[1]
 pipe2₊inp₊m_dot(t) ~ (pipe2₊m_dot(t))[6]
 pipe1₊oup₊p(t) ~ (pipe1₊p(t))[1]
 source₊oup₊m_dot(t) ~ -(pipe1₊m_dot(t))[6]
 pipe3₊inp₊m_dot(t) ~ (pipe3₊m_dot(t))[6]
 pipe2₊oup₊m_dot(t) ~ -(pipe2₊m_dot(t))[1]
 source₊m_dot(t) ~ (pipe1₊m_dot(t))[6]
 pipe3₊oup₊p(t) ~ (pipe3₊p(t))[1]
 pipe4₊inp₊m_dot(t) ~ (pipe4₊m_dot(t))[6]
 pipe4₊oup₊p(t) ~ (pipe4₊p(t))[1]
 pipe1₊inp₊m_dot(t) ~ (pipe1₊m_dot(t))[6]
 sink₊inp₊p(t) ~ (pipe4₊p(t))[1]
 (pipe4₊p(t))[6] ~ (pipe3₊p(t))[1]
 (pipe2₊p(t))[1] ~ (pipe3₊p(t))[1]
 (pipe1₊m_dot(t))[1] ~ (pipe2₊m_dot(t))[6] + (pipe3₊m_dot(t))[6]      
 (pipe2₊p(t))[6] ~ (pipe1₊p(t))[1]
 (pipe4₊m_dot(t))[1] ~ 1.0
 var"pipe2₊p(t)[1]ˍt" ~ 34.227271569789735(pipe2₊m_dot(t))[2] + 58.46715989423868(pipe2₊m_dot(t))[3] + 18.067346020157103(pipe2₊m_dot(t))[5] - 76.5345059143958(pipe2₊m_dot(t))[1] - 26.14730879497342(pipe2₊m_dot(t))[4] - 8.079962774816314(pipe2₊m_dot(t))[6]
 (pipe3₊m_dot(t))[1] ~ 0.4472135954999579(pipe3₊m_dot(t))[2] + 0.7639320225002102(pipe3₊m_dot(t))[3] + 0.23606797749978964(pipe3₊m_dot(t))[5] - 0.01742133587201403var"pipe3₊p(t)[1]ˍt" - 0.3416407864998738(pipe3₊m_dot(t))[4] - 0.1055728090000841(pipe3₊m_dot(t))[6]
 (pipe1₊p(t))[6] ~ 300000.0
 pipe4₊inp₊p(t) ~ (pipe4₊p(t))[6]
 pipe4₊ρ(t) ~ 5.807200929152149e-6(pipe4₊p(t))[1] + 5.807200929152149e-6(pipe4₊p(t))[6]
 pipe2₊oup₊p(t) ~ (pipe2₊p(t))[1]
 pipe1₊oup₊m_dot(t) ~ -(pipe1₊m_dot(t))[1]
 pipe2₊inp₊p(t) ~ (pipe2₊p(t))[6]
 (pipe3₊p(t))[6] ~ (pipe2₊p(t))[6]
 pipe2₊ρ(t) ~ 5.807200929152149e-6(pipe2₊p(t))[1] + 5.807200929152149e-6(pipe2₊p(t))[6]
 sink₊inp₊m_dot(t) ~ (pipe4₊m_dot(t))[1]
 pipe4₊oup₊m_dot(t) ~ -(pipe4₊m_dot(t))[1]
 var"pipe3₊p(t)[1]ˍt" ~ var"pipe2₊p(t)[1]ˍt"
 pipe3₊oup₊m_dot(t) ~ -(pipe3₊m_dot(t))[1]
 source₊oup₊p(t) ~ (pipe1₊p(t))[6]
 pipe1₊inp₊p(t) ~ (pipe1₊p(t))[6]
 pipe1₊ρ(t) ~ 5.807200929152149e-6(pipe1₊p(t))[1] + 5.807200929152149e-6(pipe1₊p(t))[6]
 pipe3₊inp₊p(t) ~ (pipe3₊p(t))[6]
 pipe3₊ρ(t) ~ 5.807200929152149e-6(pipe3₊p(t))[1] + 5.807200929152149e-6(pipe3₊p(t))[6]

but there are var"pipe2+p(t)[1]_t" in observed equations

var"pipe2₊p(t)[1]ˍt" ~ 34.227271569789735(pipe2₊m_dot(t))[2] + 58.46715989423868(pipe2₊m_dot(t))[3] + 18.067346020157103(pipe2₊m_dot(t))[5] - 76.5345059143958(pipe2₊m_dot(t))[1] - 26.14730879497342(pipe2₊m_dot(t))[4] - 8.079962774816314(pipe2₊m_dot(t))[6]

[YingboMa]

Why is this issue closed?

[YingboMa]

I think it's the same problem with #1726