# MTK -- removing parameters from model?

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:

1. Is my suspicion correct?
2. If so, how do I get around this?
3. 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.

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,
T_aδ ~ sol_alg,
T_ah ~ sol_alg,
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:

1. 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.

2. The really, really weird results:

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…

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.