MTK with MonteCarloMeasurements Error

Specific Domains Modelling & Simulations
, , ,
1

Hi out there!

I am trying to use MonteCarloMeasurements.jl together with MTK to solve a differential equation with uncertain parameters.

using ModelingToolkit, DifferentialEquations
using MonteCarloMeasurements
using CairoMakie

@independent_variables t
D = Differential(t)
@parameters my_param = 1.5
@variables xp(t) Tp(t)  

function my_problematic_function(xp)
    xp^1.8 
end

eqs = [
    D(xp) ~ my_param * my_problematic_function(xp),
    D(Tp) ~ (Tp + my_param * my_problematic_function(xp))
]

@named sys_model = System(eqs, t)
sys = mtkcompile(sys_model)

tspan = (0.0, 30.0)
u0_nom = [
    sys.xp => 0.75,
    sys.Tp => 293.15]
prob = ODEProblem(sys, u0_nom, tspan; jac=true)
# Redefine prob with uncertain parameters
prob_uncertain = remake(prob; u0 = u0_nom, p=[my_param=> 1.6 ± 0.3])
sol_uncertain = solve(prob_uncertain, Tsit5())

I think I ran into a bug. Whenever I take the power of a variable, I get a stackoverflow-error.

Warning: detected a stack overflow; program state may be corrupted, so further execution might be unreliable.
ERROR: StackOverflowError:

This might need a fix. I would also be interested in workarounds!
Thanks for your efforts in advance.

2

This probably needs an issue and an overload.

3

A temporary workaround could be to register the function:

function my_problematic_function(xp)
    # xp^1.8 -> domain error
    exp(1.8*log(xp)) # works
end
@register_symbolic my_problematic_function(xp)
4

Until the issue is sorted out, does it work if you avoid changing the type of the parameters when you remake? You could try to add my_param=> 1.6 ± 0.0 to the u0_nom map to make sure the problem has the intended type from the start.

5

Thanks for the Reply!

If I put my_param into the u0_nom map like you recommended, I get an

ArgumentError: invalid index: ModelingToolkitBase.ParameterIndex{SciMLStructures.Tunable, Int64}(SciMLStructures.Tunable(), 1, false) of type ModelingToolkitBase.ParameterIndex{SciMLStructures.Tunable, Int64}

If I start off with

@parameters my_param = 1.5 ± 0.1

I get a

ERROR: MethodError: no method matching Float64(::Particles{Float64, 2000})
The type `Float64` exists, but no method is defined for this combination of argument types when trying to construct it.
6

Thanks for your reply!

This workaround indeed works!

7

Here’s an example where it works

the difference is that the example defines the parameter with the type annotated, in your case it would be

T = typeof(1.0 ± 0.0)
@parameters my_param::T = 1.5 ± 0.0
8

Thanks for the link! This certainly is more elegant as it avoids using remake, but I still get the stack overflow error