MTK -- removing parameters from model?

BLI · March 31, 2023, 11:58am

My MTK implementation of a DAE for doing Monte Carlo study of parameter uncertainty doesn’t work. I think it has to do with how MTK removes parameters that do not show up explicitly in the model equations.

Here is some code:

using ModelingToolkit
using DifferentialEquations

@register_symbolic u_If(t)
@register_symbolic u_It(t)
@register_symbolic u_Twc(t)

# Give the inputs some constant values, e.g.,
u_If(t) = 800
u_It(t) = 5e3
u_TWc(t) = 4

# 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, 
    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)

# 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))

# Input equations/input system
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/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]
@named sys_p = ODESystem(eqs_p)

# Total system
sys = extend(sys_i, sys_p)
sys_simp = structural_simplify(sys)

# Converting symbolic model to numeric model
tspan = (0, 600)
#
prob = ODEProblem(sys_simp, [], tspan)

The last statement crashes with error message:

ArgumentError: Any[UA_x / (ĉ_pw*ṁ_w), 0.017674089784376106UA_x, UA_x / (ĉ_pw*ṁ_w) - 0.017674089784376106UA_x] are either missing from the variable map or missing from the system's states/parameters list.

My suspicion is that MTK has removed parameters in the @parameters list that do not explicitly appear in the model. This includes UA_x, ĉ_pw, ṁ_w.

Questions:

Is my suspicion correct?
If so, how do I get around this?
I could of course remove the parameters UA_x, ĉ_pw, ṁ_w and instead just specify the two Stanton numbers N_St_a and N_St_w. Then I think the model would work. However, I want to study uncertainty in the heat transfer coefficient UA_x, which is a single parameter.

BLI · March 31, 2023, 2:54pm

Hm… I partially got around the problem by defining the Stanton numbers as variables – this means that I compute N_St_a, etc. as time functions by adding the relationship between “parameters” (turned variables) N_St_a, N_St_w, and N_St_Δ and the real parameters:

# 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))
#...
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)

Still some weird behavior, though. Observing that the 3 algebraic relationships give linear equations for the temperatures T_ac, T_aδ, and T_ah, I use Symbolics to solve the 3 algebraic equations for these temperatures:

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]
# ...
sol_alg = simplify.(Symbolics.solve_for(eqs_alg, [T_ac, T_aδ, T_ah]))
# ...
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(sys_ode_simp)

Again, I can solve both sys_simp and sys_ode_simp without getting error messages.

One strange observation and one really weird result:

Solving sys_ode_simp (which does not have a mass matrix) gives identical solutions as to solving sys_simp. HOWEVER, sys_ode_simp is some 20% slower to solve than the version with mass matrix, and allocates some 40% more memory. That is somewhat surprising.
The really, really weird results:

When I solve sys_simp, N_St_a, N_St_w, and N_St_Δ appear when I do observed(sys_simp), and I can plot the results for all three – they are obviously constants.
When I solve sys_ode_simp, N_St_a, N_St_w, and N_St_Δ appear when I do observed(sys_simp). I can plot N_St_a and N_St_w,

image396×639 31.2 KB

but I get an error message when I try to plot N_St_Δ:

In fact, I cannot plot any variable below N_St_a in the observed(sys_ode_simp) list…

BLI · March 31, 2023, 2:56pm

Anyway, my postings have gotten somewhat convoluted… Perhaps I should move the discussion to another forum (Issues) and provide the experts with my complete code.

My reason for raising the question here was in case some of you users would have tricks to share.