MonteCarloMeasurements & MTK, again

Specific Domains Modelling & Simulations
1

I’m trying to use MonteCarloMeasurements with MTK. I have some problems wrt. (i) the number of particles, and (ii) interpretation/presentation of the results.

Here is a simple RC circuit example:

using Plots, LaTeXStrings
using ModelingToolkit, DifferentialEquations
using MonteCarloMeasurements

# 1. Independent variables
@variables t 
# 2. Differential operator
Dt = Differential(t)
# 3. Parameters
pars_RC = @parameters begin 
    C=8,            [description="Capacitance, C"]
    R=0.1,          [description="Resistance, Ω"]
end
# 4. Dependent variables 
vars_RC = @variables begin
    i(t),           [description="Total current, A"]
    i_c(t),         [description="Capacitor current, A"]
    i_r(t),         [description="Resistor current, A"]
    (u(t)=0),       [description="Voltage, V"]
end


# RC circuit
eqs_a = [Dt(u) ~ i_c/C,
        i_r ~ u/R,
        i ~ i_r + i_c]

@named mod_RC_a = ODESystem(eqs_a,t,vars_RC,pars_RC)

# Specifying inputs
@register_symbolic i_u(t)

eqs_i = [i ~ i_u(t)]
@named i_input = ODESystem(eqs_i,t)

# Causal system
mod_RC = extend(i_input, mod_RC_a) |> structural_simplify

i_u(t) = 10

tspan=(0,5)
prob = ODEProblem(mod_RC, [],tspan)
sol = solve(prob)
plot(sol, lw=2.5)

This produces the following plot:

image
image724×468 14.3 KB

Next, using MonteCarloMeasurements:

prob = remake(prob; u0=[u => -0.2..0.2], p=[C=> 8±1])
sol = solve(prob)
plot(sol)

produces:

image
image721×477 49.2 KB

[Warning: Using arrays or dicts to store parameters of different types can hurt performance. Consider using tuples instead. – I neglected this advice…]

Question 1: In order for this to work, it is required that at least one parameter and at least one state/initial value simultaneously are specified with uncertainty, right? […to make sure that the quantities are properly promoted to the new type…]

According to the documentations, the above syntax for specifying uncertainty produces 2000 particles.

The following alternative description also works:

Np = 10
prob = remake(prob; u0=[u => 0.2*Particles(Np)], p=[C=> 8 + Particles(Np)])
sol = solve(prob)
plot(sol)

with result:

image
image711×480 35.6 KB

Question 2:: This only works if I specify the same number of particles in every uncertainty, right? [Np above is used in both cases of calling Particles…]

Question 3: Clearly, the blue region in the two plots above are not the “shadows” of plotting multiple parameters – in the last plot, one may discern a “small” number of particle trajectories, probably 10.

  • Does the blue region represent a ribbon plot of the computed mean +/- the computed standard deviation of the particles?
  • Does the thicker blue line in the center of the ribbon plot represent the mean?

I observe that I can change the color of the filled area (“ribbon”?) with setting fill color (fc=...), the line color by specifying lc, and the line thickness by specifying lw.

Question 4: Regarding solutions…

  • Is it possible to specify the line width, line color, line alpha value, and line style individually for the mean and the particles?
2

In DiffEq it’s enough that the initial condition has the correct type, I’m not sure if MTK requires something different here.

Correct

almost, they are quantiles, see plotting for more info.

yes

No. There are also two plot recipes that are interacting here, the recipe for solutions and the one for particles. I know that this interaction is not always working perfectly. If you want full control, you can issue two plot commands, one without sample trajectories and one with.

If you want to use a custom number of particles by default, you can define your own operator, e.g.,

julia> import MonteCarloMeasurements # Don't use using to avoid conflicting import of ..

julia> (..)(a, b) = MonteCarloMeasurements.Particles(10, Uniform(a,b))
.. (generic function with 1 method)

julia> 1 .. 2
1.5 ± 0.303 Particles{Float64, 10}
1 Like
3

Ah. OK – so presumably defaults to 0.025/0.975 quantiles for ribbon?

4

yes