MTK Bug? Some `observed` are not computed?

B-LIE commented 1 year ago

I came across a weird observation with MTK -- a bug? related to observed values that are not computed.

A birds eye view of what happened is the following...

I have an index 1 DAE consisting of 3 differential equations and 3 algebraic equations.
I posed the model as this DAE using MTK. I then added "artificial" variables for some parameters (as discussed in a separate Julia Discourse post) for doing uncertainty analysis.
MTK produces a structural_simplify model where the states are the differential variables + one of the algebraic variables.
Everything seems to work

Then I observed that the 3 algebraic equations, in fact, are linear in the algebraic variables. So I wanted to to "manually" reduce the DAE to index 0 so that I could avoid using a solver with mass matrix, and hopefully speed up the simulation.

I took out the algebraic equations, and used Symbolics.solve_for() to find the algebraic variables expressed by the differential variables.
I used the same model as for the index 1 DAE, except that I replaced the implicit algebraic equations by explicit expressions for the algebraic variables. I included "equations" for giving values to some "parameters" (as discussed in separate Julia Discourse post).c
Everything appears to work.

I then start to plot the results from the two models to confirm that the results are the same. Then the "shock":

With the manually reduced index model, some of the variables that appears as observed can not be evaluated or plotted.

QUESTION: is this a bug? What is going wrong?

[I provide the codes in a follow-up message below.]

B-LIE commented 1 year ago

Here is the index 1 DAE -- this seems to work as it should (sorry for complex example...)

using ModelingToolkit
using DifferentialEquations

# Registering input variables
@register_symbolic u_If(t)
@register_symbolic u_It(t)
@register_symbolic u_Twc(t)

# Example of possible input functions that should work
u_If(t) = 900
u_It(t) = 5000
u_TWc(t) = 4

# Parameters, variables, etc.
# -- Parameters with defaults
@parameters (UA_r2δ=2.7, UA_s2Fe=20, UA_Fe2a=14.3, UA_x=44.4, Q̇_Fe_σ=212, Q̇_f_σ=422.4, 
    m_r=9260, m_s=6827,  m_Fe=71200, R_r=0.127e-3,  R_s=1.95e-6, ṁ_a=49.2, ṁ_w=54.2,
    ĉ_pCu=0.385,  ĉ_pFe=0.465,  ĉ_pa=1.15, ĉ_pw = 4.18)

# -- Independent variable and differentiation operator
@variables t 
Dt = Differential(t)

# -- Dependent variables
@variables (T_r(t)=28, T_s(t)=28, T_Fe(t)=28, T_ac(t)=14, T_aδ(t)=14, T_ah(t)=22, I_f(t), I_t(t), T_wc(t),
            N_St_a(t) = UA_x/(ĉ_pa*ṁ_a), N_St_w(t) = UA_x/(ĉ_pw*ṁ_w), N_St_Δ(t) = UA_x/(ĉ_pw*ṁ_w)-UA_x/(ĉ_pa*ṁ_a))

# Input equations
eqs_i = [I_f ~ u_If(t), I_t ~ u_It(t), T_wc ~ u_Twc(t)]
@named sys_i = ODESystem(eqs_i)

# Model equations and System
eqs_p = [Dt(T_r) ~ (1.1*R_r*I_f^2 - UA_r2δ*(T_r-T_aδ))/(m_r*ĉ_pCu),
            Dt(T_s) ~ (3*R_s*I_t^2 - UA_s2Fe*(T_s-T_Fe))/(m_s*ĉ_pCu),
            Dt(T_Fe) ~ (UA_s2Fe*(T_s-T_Fe) - UA_Fe2a*(T_Fe-T_ah) + Q̇_Fe_σ)/(m_Fe*ĉ_pFe),
            0 ~ ṁ_a*ĉ_pa*(T_ac-T_aδ) + UA_r2δ*(T_r-T_aδ) + Q̇_f_σ,
            0 ~ ṁ_a*ĉ_pa*(T_aδ-T_ah) + UA_Fe2a*(T_Fe-T_ah),
            0 ~ T_ac*(N_St_w-N_St_a*exp(-N_St_Δ)) - N_St_Δ*T_ah - N_St_a*(1-exp(-N_St_Δ))*T_wc,
            0 ~ N_St_a - UA_x/(ĉ_pa*ṁ_a),
            0 ~ N_St_w - UA_x/(ĉ_pw*ṁ_w),
            0 ~ N_St_Δ - (N_St_w - N_St_a)]
@named sys_p = ODESystem(eqs_p)

sys = extend(sys_i, sys_p)
sys_simp = structural_simplify(sys)

# States: OBSERVE that even though algebraic equations are *linear*, one of the algebraic variables
#  is kept as a "state"
states(sys_simp)

# Solving the index 1 DAE
tspan = (0.0, 500)
#
prob = ODEProblem(sys_simp, [], tspan)
#
sol = solve(prob)

NOTE: for this implementation, I can plot all states and all observed variables.

B-LIE commented 1 year ago

And here is the problematic case where I manually reduce the above index 1 DAE to an index 0 DAE. The formulation can still be solved. And I can plot (and evaluate) the states. But some variables listed as observed can not be plotted (not evaluated)... WHY, OH WHY? :-) [First part of code is quite similar to the above, but I include the complete code...]

using ModelingToolkit
using DifferentialEquations
using Symbolics

# Registering input variables
@register_symbolic u_If(t)
@register_symbolic u_It(t)
@register_symbolic u_Twc(t)

# Example of possible input functions that should work
u_If(t) = 900
u_It(t) = 5000
u_TWc(t) = 4

# Parameters, variables, etc.
# -- Parameters with defaults
@parameters (UA_r2δ=2.7, UA_s2Fe=20, UA_Fe2a=14.3, UA_x=44.4, Q̇_Fe_σ=212, Q̇_f_σ=422.4, 
    m_r=9260, m_s=6827,  m_Fe=71200, R_r=0.127e-3,  R_s=1.95e-6, ṁ_a=49.2, ṁ_w=54.2,
    ĉ_pCu=0.385,  ĉ_pFe=0.465,  ĉ_pa=1.15, ĉ_pw = 4.18)

# -- Independent variable and differentiation operator
@variables t 
Dt = Differential(t)

# -- Dependent variables
@variables (T_r(t)=28, T_s(t)=28, T_Fe(t)=28, T_ac(t)=14, T_aδ(t)=14, T_ah(t)=22, I_f(t), I_t(t), T_wc(t),
            N_St_a(t) = UA_x/(ĉ_pa*ṁ_a), N_St_w(t) = UA_x/(ĉ_pw*ṁ_w), N_St_Δ(t) = UA_x/(ĉ_pw*ṁ_w)-UA_x/(ĉ_pa*ṁ_a))

# Input equations
eqs_i = [I_f ~ u_If(t), I_t ~ u_It(t), T_wc ~ u_Twc(t)]
@named sys_i = ODESystem(eqs_i)

#== Code that is different from index 1 DAE ==

# Algebraic equations
eqs_alg = [0 ~ ṁ_a*ĉ_pa*(T_ac-T_aδ) + UA_r2δ*(T_r-T_aδ) + Q̇_f_σ,
    0 ~ ṁ_a*ĉ_pa*(T_aδ-T_ah) + UA_Fe2a*(T_Fe-T_ah),
    0 ~ T_ac*(N_St_w-N_St_a*exp(-N_St_Δ)) - N_St_Δ*T_ah - N_St_a*(1-exp(-N_St_Δ))*T_wc]

# Solving algebraic equations wrt. algebraic variables
sol_alg = simplify.(Symbolics.solve_for(eqs_alg, [T_ac, T_aδ, T_ah]))

# Equations, etc.  for index-reduced case
eqs_ode_p = [Dt(T_r) ~ (1.1*R_r*I_f^2 - UA_r2δ*(T_r-T_aδ))/(m_r*ĉ_pCu),
            Dt(T_s) ~ (3*R_s*I_t^2 - UA_s2Fe*(T_s-T_Fe))/(m_s*ĉ_pCu),
            Dt(T_Fe) ~ (UA_s2Fe*(T_s-T_Fe) - UA_Fe2a*(T_Fe-T_ah) + Q̇_Fe_σ)/(m_Fe*ĉ_pFe),
            T_ac ~ sol_alg[1],
            T_aδ ~ sol_alg[2],
            T_ah ~ sol_alg[3],
            N_St_a ~ UA_x/(ĉ_pa*ṁ_a),
            N_St_w ~ UA_x/(ĉ_pw*ṁ_w),
            N_St_Δ ~ N_St_w - N_St_a]
@named sys_ode_p = ODESystem(eqs_ode_p)

sys_ode = extend(sys_i, sys_ode_p)
sys_ode_simp = structural_simplify(sys_ode)

# States: OBSERVE that this time, the states only contain the differential variables
states(sys_ode_simp)

# Solving the index 0 DAE
tspan = (0.0, 500)
#
prob_ode = ODEProblem(sys_ode_simp, [], tspan)
#
sol_ode= solve(prob_ode)

Also this time, I can plot/evaluate all states(sys_ode_simp). But I can only plot/evaluate the 5 first elements of observed(sys_ode_simp) -- if I try to plot the 4 last elements of observed(sys_ode_simp), I get an error message...

Etc.

ChrisRackauckas commented 9 months ago

@YingboMa was this handled?

YingboMa commented 9 months ago

Yes,

julia> sol = solve(prob, Rodas5P());

julia> sol[T_ah]
14-element Vector{Float64}:
 22.05082832708408
 22.050828329810468
 22.05082835707445
 22.050828629714264
 22.050831356117264
 22.0508586206319
 22.05113131422503
 22.053863077758404
 22.081646785805003
 22.21711801414909
 22.45612655243214
 22.76886779735181
 23.194042603627985
 23.717656203469453

Note that the MWE is wrong. It should be u_Twc(t) = 4 not u_TWc(t) = 4.

SciML / ModelingToolkit.jl

MTK Bug? Some `observed` are not computed? #2136