MTK mtkcompile "fully_determined" kwarg unexpected behavior

Specific Domains Modelling & Simulations
1

Hi all,

I’m not sure if bug or not, but setting the fully_determined keyword argument to true on mtkcompile results in a compilation error, whereas setting it to false results in what seems to be a valid system, with an equal number of equations and unknowns, which can then be solved.

Here is code to reproduce.

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEq
using DataInterpolationsND

itp_dims = (LinearInterpolationDimension([0.1, 10.0]), LinearInterpolationDimension([300.0, 310.0]))

pars = @parameters begin
    V_internal = 10.0, []
    (liq_PT_2_U_itp::NDInterpolation)(..) = NDInterpolation([112.55 154.35; 111.74 153.25], itp_dims), []
    (liq_PT_2_D_itp::NDInterpolation)(..) = NDInterpolation([996.56 993.38; 1001.0 997.7], itp_dims), []
    (gas_PT_2_D_itp::NDInterpolation)(..) = NDInterpolation([939.84 971.00; 940.13 971.30], itp_dims), []
    (gas_PT_2_U_itp::NDInterpolation)(..) = NDInterpolation([0.16039 0.15522; 1.5971 1.5458], itp_dims), []
end

p_init = 0.3
T_liq_init = 300
T_gas_init = 300
m_liq_init = 1000
d_liq_init = liq_PT_2_D_itp(p_init, T_liq_init)
V_liq_init = m_liq_init / d_liq_init
V_ull_init = V_internal - V_liq_init
U_liq_init = m_liq_init * liq_PT_2_U_itp(p_init, T_liq_init)
d_gas_init = gas_PT_2_D_itp(p_init, T_gas_init)
m_gas_init = d_gas_init * V_ull_init
U_gas_init = m_gas_init * gas_PT_2_U_itp(p_init, T_gas_init)

vars = @variables begin
    V_ull(t), [description = "ullage volume"]
    V_liq(t), [description = "liquid volume"]
    p(t) = p_init, [description = "pressure"]
    m_gas(t) = m_gas_init, [description = "fluid mass in volume"]
    U_gas(t) = U_gas_init, [description = "extensive internal energy of gas"]
    T_gas(t), [description = "temperature of gas"]
    d_gas(t), [description = "density of gas"]
    m_liq(t) = m_liq_init, [description = "fluid mass in volume"]
    U_liq(t) = U_liq_init, [description = "extensive internal energy of liquid"]
    T_liq(t), [description = "temperature of liquid"]
    d_liq(t), [description = "density of liquid"]
    V_liq_dot(t), [description = "liquid volume rate of change"]
end

eqs = [
    V_internal ~ V_ull + V_liq,
    d_gas ~ m_gas / V_ull,
    d_gas ~ gas_PT_2_D_itp(p, T_gas),
    U_gas ~ m_gas * gas_PT_2_U_itp(p, T_gas),
    D(m_gas) ~ 0.0,
    D(U_gas) ~ 1.0e3 * p * V_liq_dot,
    d_liq ~ liq_PT_2_D_itp(p, T_liq),
    U_liq ~ m_liq * liq_PT_2_U_itp(p, T_liq),
    D(m_liq) ~ 10.0,
    D(U_liq) ~ 10.0 * 150.0 - 1.0e3 * p * V_liq_dot,

    # chain rule D(m_liq/d_liq) to get volume rate of change
    V_liq_dot ~ (D(m_liq) * d_liq - m_liq * D(d_liq)) / d_liq^2,
    V_liq ~ m_liq / d_liq,
]

@named _sys = System(eqs, t, vars, pars)
# this works, resulting in system with 10 equations and 10 unknowns.
sys = mtkcompile(_sys; fully_determined=false)

# this errors out, resulting in taking derivative wrt T_gas and p_gas instead of time:
# Encountered operator `Differential(p(t), 1)` which has different independent variable than the one used in the system `t`.

# sys = mtkcompile(_sys; fully_determined=true)

guesses = [T_gas => 300, T_liq => 300, d_liq => 1000, V_ull => 9.0]
prob = ODEProblem(sys, [], (0, 10); guesses=guesses, use_scc=false)
# problem solves
sol = solve(prob; abstol=1e-6)
sol.retcode

the stack trace of error

ERROR: Encountered operator `Differential(p(t), 1)` which has different independent variable than the one used in the system `t`.

Context:
-U_gasˍt(t) + (T_gasˍt(t)*Differential(T_gas(t), 1)(gas_PT_2_U_itp(p(t), T_gas(t))) + pˍt(t)*Differential(p(t), 1)(gas_PT_2_U_itp(p(t), T_gas(t))))*m_gas(t)
Stacktrace:
  [1] validate_operator(op::Differential, args::SymbolicUtils.SmallVec{…}, iv::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{…}; context::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{…})
    @ ModelingToolkitBase C:\Users\mark.garnett\.julia\packages\ModelingToolkitBase\Rcink\src\utils.jl:722
  [2] validate_operator
    @ C:\Users\mark.garnett\.julia\packages\ModelingToolkitBase\Rcink\src\utils.jl:721 [inlined]
  [3] collect_vars!(unknowns::OrderedCollections.OrderedSet{…}, parameters::OrderedCollections.OrderedSet{…}, expr::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{…}, iv::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{…}, ::Type{…}; depth::Int64)
    @ ModelingToolkitBase C:\Users\mark.garnett\.julia\packages\ModelingToolkitBase\Rcink\src\utils.jl:833
  [4] collect_vars!
    @ C:\Users\mark.garnett\.julia\packages\ModelingToolkitBase\Rcink\src\utils.jl:798 [inlined]
  [5] #collect_vars!#37
    @ C:\Users\mark.garnett\.julia\packages\ModelingToolkitBase\Rcink\src\utils.jl:892 [inlined]
  [6] collect_vars!
    @ C:\Users\mark.garnett\.julia\packages\ModelingToolkitBase\Rcink\src\utils.jl:887 [inlined]
  [7] collect_scoped_vars!(unknowns::OrderedCollections.OrderedSet{…}, parameters::OrderedCollections.OrderedSet{…}, sys::System, iv::SymbolicUtils.BasicSymbolicImpl.var"typeof(BasicSymbolicImpl)"{…}, ::Type{…}; depth::Int64)
    @ ModelingToolkitBase C:\Users\mark.garnett\.julia\packages\ModelingToolkitBase\Rcink\src\utils.jl:688
  [8] discover_globalscoped(sys::System)
    @ ModelingToolkitBase C:\Users\mark.garnett\.julia\packages\ModelingToolkitBase\Rcink\src\systems\abstractsystem.jl:643
  [9] complete(sys::System; split::Bool, flatten::Bool, add_initial_parameters::Bool)
    @ ModelingToolkitBase C:\Users\mark.garnett\.julia\packages\ModelingToolkitBase\Rcink\src\systems\abstractsystem.jl:680
 [10] complete
    @ C:\Users\mark.garnett\.julia\packages\ModelingToolkitBase\Rcink\src\systems\abstractsystem.jl:677 [inlined]
 [11] mtkcompile(sys::System; additional_passes::Tuple{}, inputs::Vector{…}, outputs::Vector{…}, disturbance_inputs::Vector{…}, split::Bool, kwargs::@Kwargs{…})
    @ ModelingToolkitBase C:\Users\mark.garnett\.julia\packages\ModelingToolkitBase\Rcink\src\systems\systems.jl:104
 [12] top-level scope
    @ c:\Users\mark.garnett\EngineModeling\ignore\tank_debug_standalone.jl:65

It seems to me like it ought to be possible for the terms of which the differentials need to be taken

-U_gasˍt(t) + (T_gasˍt(t)*Differential(T_gas(t), 1)(gas_PT_2_U_itp(p(t), T_gas(t))) + pˍt(t)*Differential(p(t), 1)(gas_PT_2_U_itp(p(t), T_gas(t))))*m_gas(t)

the chain rule is being applied, and

Differential(T_gas(t), 1)(gas_PT_2_U_itp(p(t), T_gas(t)))

should be available via

gas_PT_2_U_itp(p(t), T_gas(t); derivative_orders=(0,1)

and the derivative w.r.t. p(t) similarly.

It doesn’t seem like in principle this is wrong.

installed packages

  [a93c6f00] DataFrames v1.8.2
  [4f1ef021] DataInterpolationsND v0.1.2
  [8bb1440f] DelimitedFiles v1.9.1
⌅ [0c46a032] DifferentialEquations v7.17.0
⌅ [5b8099bc] DomainSets v0.7.18
⌃ [e9467ef8] GLMakie v0.13.10
⌃ [c27321d9] Glob v1.4.0
  [a98d9a8b] Interpolations v0.16.2
⌃ [7ed4a6bd] LinearSolve v3.76.0
⌃ [961ee093] ModelingToolkit v11.26.0
  [16a59e39] ModelingToolkitStandardLibrary v2.28.0
  [15e1cf62] NPZ v0.4.3
⌃ [8913a72c] NonlinearSolve v4.19.0
⌅ [1dea7af3] OrdinaryDiffEq v6.111.0
⌅ [127b3ac7] OrdinaryDiffEqNonlinearSolve v1.28.0
  [91a5bcdd] Plots v1.41.6
  [efcf1570] Setfield v1.1.2
  [1986cc42] Unitful v1.28.0
  [de0858da] Printf v1.11.0
2

Investigating further, it seems to me that this is a bug, or perhaps un-implemented functionality in DataInterpolationsND, Symbolics, or MTK. Not sure exactly where.

If we replace the linear DataInterpolationND functions with a hand-rolled 4 point bilinear interpolation, MTK can successfully apply the chain rule and compile the system. The exact same data is passed to it as was passed to NDInterpolation.

So in principle this ought to be able to make this work for an actual interpolation function, with more data points.

Can anyone help pinpoint better where the problem is occurring?

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using OrdinaryDiffEq
using DataInterpolationsND

p_data = [0.1, 10.0]
T_data = [300.0, 310.0]
liq_U_data = [112.55 154.35; 111.74 153.25]
liq_D_data = [996.56 993.38; 1001.0 997.7]
gas_U_data = [939.84 971.00; 940.13 971.30]
gas_D_data = [0.16039 0.15522; 1.5971 1.5458]

function itp(x, y, X, Y, Q)
    ([X[2] - x, x - X[1]]' * Q * [Y[2] - y, y - Y[1]]) / ((X[2] - X[1]) * (Y[2] - Y[1]))
end

liq_PT_2_U_itp(p, T) = itp(p, T, p_data, T_data, liq_U_data)
liq_PT_2_D_itp(p, T) = itp(p, T, p_data, T_data, liq_D_data)
gas_PT_2_D_itp(p, T) = itp(p, T, p_data, T_data, gas_U_data)
gas_PT_2_U_itp(p, T) = itp(p, T, p_data, T_data, gas_D_data)


# extrap=(ExtendExtrap(), ExtendExtrap())
pars = @parameters begin
    V_internal = 10.0, []
end

p_init = 0.3
T_liq_init = 300
T_gas_init = 300
m_liq_init = 1000
d_liq_init = liq_PT_2_D_itp(p_init, T_liq_init)
V_liq_init = m_liq_init / d_liq_init
V_ull_init = V_internal - V_liq_init
U_liq_init = m_liq_init * liq_PT_2_U_itp(p_init, T_liq_init)
d_gas_init = gas_PT_2_D_itp(p_init, T_gas_init)
m_gas_init = d_gas_init * V_ull_init
U_gas_init = m_gas_init * gas_PT_2_U_itp(p_init, T_gas_init)

vars = @variables begin
    V_ull(t), [description = "ullage volume"]
    V_liq(t), [description = "liquid volume"]
    p(t) = p_init, [description = "pressure"]
    m_gas(t) = m_gas_init, [description = "fluid mass in volume"]
    U_gas(t) = U_gas_init, [description = "extensive internal energy of gas"]
    T_gas(t), [description = "temperature of gas"]
    d_gas(t), [description = "density of gas"]
    m_liq(t) = m_liq_init, [description = "fluid mass in volume"]
    U_liq(t) = U_liq_init, [description = "extensive internal energy of liquid"]
    T_liq(t), [description = "temperature of liquid"]
    d_liq(t), [description = "density of liquid"]
    # V_liq_dot(t), [description = "liquid volume rate of change"]

end

eqs = [
    V_internal ~ V_ull + V_liq,
    d_gas ~ m_gas / V_ull,
    d_gas ~ gas_PT_2_D_itp(p, T_gas),
    U_gas ~ m_gas * gas_PT_2_U_itp(p, T_gas),
    D(m_gas) ~ 0.1,
    D(U_gas) ~ 1.0e3 * p * D(V_liq),
    d_liq ~ liq_PT_2_D_itp(p, T_liq),
    U_liq ~ m_liq * liq_PT_2_U_itp(p, T_liq),
    D(m_liq) ~ 10.0,
    D(U_liq) ~ 10.0 * 150.0 - 1.0e3 * p * D(V_liq),
    # MTK applies chain rule, don't need this
    # D(V_liq) ~ (D(m_liq) * d_liq - m_liq * D(d_liq)) / d_liq^2,
    V_liq ~ m_liq / d_liq,
]

@named _sys = System(eqs, t, vars, pars)
# now it works, resulting in system with 9 equations and 9 unknowns.

sys = mtkcompile(_sys; fully_determined=true)
3

Following up in Differentiating through DataInterpolationsND functions throws error · Issue #69 · SciML/DataInterpolationsND.jl · GitHub

4

Thank you Chris, I read through the issues generated. Looks like we are on the way to a fix so that’s great!

Looks like to make it work with callable parameters, I’ll need to wait for this issue in Symbolics to be resolved, but in the meanwhile I can hard-code the interps in like claude-bot suggests or maybe figure out some macro-magic to make it so that the right interp functions for a given fluid are substituted in.