@mtkmodel Simple begin
    @variables begin
        x(t) = 0
    end
    @equations begin
        D(x) ~ 1
    end
end

@mtkbuild sys = Simple()
initprob = ModelingToolkit.InitializationProblem(sys, 0.0)
initsol = solve(initprob)
u0 = initsol[unknowns(sys)]
prob = ODEProblem(sys, u0, (0, 0))  #OK, no warning

function System()
    pars = @parameters begin
        rho_0 = 1000
        beta = 2e9
        A = 0.1
        m = 100
        L = 1
        p_s = 100e5
        p_r = 10e5
        C = 2.7*0.5
        c = 1000
        A_p = 0.00094
    end   
    vars = @variables begin
        p_1(t) = p_s
        p_2(t) = p_r
        x(t)=0
        dx(t)=0
        ddx(t), [guess=0]
        rho_1(t), [guess=rho_0]
        rho_2(t), [guess=rho_0]
        drho_1(t), [guess=0]
        drho_2(t), [guess=0]
        dm_1(t), [guess=0]
        dm_2(t), [guess=0]
    end
    
    # let -----
    u_1 = dm_1/(rho_0*A_p)
    u_2 = dm_2/(rho_0*A_p)

    eqs = [
        D(x) ~ dx
        D(dx) ~ ddx
        D(rho_1) ~ drho_1
        D(rho_2) ~ drho_2
        +dm_1 ~ drho_1*(L+x)*A + rho_1*dx*A
        -dm_2 ~ drho_2*(L-x)*A - rho_2*dx*A
        rho_1 ~ rho_0*(1 + p_1/beta)
        rho_2 ~ rho_0*(1 + p_2/beta)
        m*ddx ~ (p_1 - p_2)*A - c*dx
        (p_s - p_1) ~ C*rho_0*(u_1)*abs(u_1)
        (p_2 - p_r) ~ C*rho_0*(u_2)*abs(u_2)
    ]

    return ODESystem(eqs, t, vars, pars; name)
end

@mtkbuild sys = System()
initprob = ModelingToolkit.InitializationProblem(sys, 0.0)
initsol = solve(initprob)
u0 = initsol[unknowns(sys)]
prob = ODEProblem(sys, u0, (0, 0)) # Warning: Initialization system is overdetermined. 4 equations for 0 unknowns.

julia> prob = ODEProblem(sys, [], (0, 0))
ODEProblem with uType Vector{Float64} and tType Int64. In-place: true
timespan: (0, 0)
u0: 6-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.0
 0.0