# MTK -- removing parameters from model?

**URL:** <https://discourse.julialang.org/t/mtk-removing-parameters-from-model/96887>\
**Category:** Modelling & Simulations\
**Created:** [March 31, 2023, 11:58am UTC](https://discourse.julialang.org/t/mtk-removing-parameters-from-model/96887 "2023-03-31T11:58:19Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [March 31, 2023, 11:58am UTC](https://discourse.julialang.org/t/mtk-removing-parameters-from-model/96887/1 "2023-03-31T11:58:19Z")

</div>

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:

```julia
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, ṁ_a=49.2, ṁ_w=54.2,
    ĉ_pCu=0.385, ĉ_pFe=0.465, ĉ_pa=1.15, ĉ_pw = 4.18, 
    N_St_a=UA_x/(ĉ_pa*ṁ_a), N_St_w=UA_x/(ĉ_pw*ṁ_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*ĉ_pCu),
            Dt(T_s) ~ (3*R_s*I_t^2 - UA_s2Fe*(T_s-T_Fe))/(m_s*ĉ_pCu),
            Dt(T_Fe) ~ (UA_s2Fe*(T_s-T_Fe) - UA_Fe2a*(T_Fe-T_ah) + Q̇_Fe_σ)/(m_Fe*ĉ_pFe),
            0 ~ ṁ_a*ĉ_pa*(T_ac-T_aδ) + UA_r2δ*(T_r-T_aδ) + Q̇_f_σ,
            0 ~ ṁ_a*ĉ_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:

```julia
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`, `ĉ_pw`, `ṁ_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`, `ĉ_pw`, `ṁ_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.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [March 31, 2023, 2:54pm UTC](https://discourse.julialang.org/t/mtk-removing-parameters-from-model/96887/2 "2023-03-31T14:54:03Z")

</div>

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:

```julia
# 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/(ĉ_pa*ṁ_a), N_St_w(t) = UA_x/(ĉ_pw*ṁ_w), N_St_Δ(t) = UA_x/(ĉ_pw*ṁ_w)-UA_x/(ĉ_pa*ṁ_a))
#...
eqs_p = [Dt(T_r) ~ (1.1*R_r*I_f^2 - UA_r2δ*(T_r-T_aδ))/(m_r*ĉ_pCu),
            Dt(T_s) ~ (3*R_s*I_t^2 - UA_s2Fe*(T_s-T_Fe))/(m_s*ĉ_pCu),
            Dt(T_Fe) ~ (UA_s2Fe*(T_s-T_Fe) - UA_Fe2a*(T_Fe-T_ah) + Q̇_Fe_σ)/(m_Fe*ĉ_pFe),
            0 ~ ṁ_a*ĉ_pa*(T_ac-T_aδ) + UA_r2δ*(T_r-T_aδ) + Q̇_f_σ,
            0 ~ ṁ_a*ĉ_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/(ĉ_pa*ṁ_a),
            0 ~ N_St_w - UA_x/(ĉ_pw*ṁ_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:

```julia
eqs_alg = [0 ~ ṁ_a*ĉ_pa*(T_ac-T_aδ) + UA_r2δ*(T_r-T_aδ) + Q̇_f_σ,
    0 ~ ṁ_a*ĉ_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*ĉ_pCu),
            Dt(T_s) ~ (3*R_s*I_t^2 - UA_s2Fe*(T_s-T_Fe))/(m_s*ĉ_pCu),
            Dt(T_Fe) ~ (UA_s2Fe*(T_s-T_Fe) - UA_Fe2a*(T_Fe-T_ah) + Q̇_Fe_σ)/(m_Fe*ĉ_pFe),
            T_ac ~ sol_alg[1],
            T_aδ ~ sol_alg[2],
            T_ah ~ sol_alg[3],
            N_St_a ~ UA_x/(ĉ_pa*ṁ_a),
            N_St_w ~ UA_x/(ĉ_pw*ṁ_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:

- 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`,  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/d/a/da50cc19e9dcf3c40c1bb5d4b8a3d7c63bcc6d7e.png)

but I get an error message when I try to plot `N_St_Δ`:  
 ![image](https://global.discourse-cdn.com/julialang/original/3X/f/f/ffd8f7b3a2153d9913dd53ff86172b62e65a9558.png)

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

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [March 31, 2023, 2:56pm UTC](https://discourse.julialang.org/t/mtk-removing-parameters-from-model/96887/3 "2023-03-31T14:56:44Z")

</div>

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.