ModelingToolkit and Noise Processes for uniform noise definition

NinaOmejc · February 25, 2025, 12:44pm

Hi,
I am using MTK to set up a hierarchical dynamical system. As the input to the (otherwise ODE) system can also be random uniform noise, I need to set up SDESystem (instead of ODESystem). However, I cannot find out how to define uniformly distributed noise with MTK, only the Brownian noise, by using the macro @brownian.

TLTR; Is it possible to do it and if someone has an example of how to do it, I would really greatly appreciate it. Thank you.

A longer description of what I have tried:
i) Simulation with DifferentialEquations.jl to have the ground truth. (Although, interestingly, the simulation also works with the ODEProblem definition here, not necessarily the SDEProblem)

Code and output for simulating with DifferentialEquations.jl


using DifferentialEquations
using Plots
using Interpolations

# Sigmoid function
function sigmoid(V, e0, r, θ)
    return 2 * e0 / (1 + exp(r * (θ - V)))
end

# Jansen-Rit equations
function jansen_rit!(du, u, p, t)
    A0, a0, A1, a1, A2, a2, C1, C2, C3, C4, e0, r, v0, I_ext = p
    
    du[1] = u[4]
    du[2] = u[5]
    du[3] = u[6]
    du[4] = A0 * a0 * sigmoid(u[2] - u[3], e0, r, v0) - 2 * a0 * u[4] - a0^2 * u[1]
    du[5] = A1 * a1 * (I_ext(t) + C2 * sigmoid(C1 * u[1], e0, r, v0)) - 2 * a1 * u[5] - a1^2 * u[2]
    du[6] = A2 * a2 * (C4 * sigmoid(C3 * u[1], e0, r, v0)) - 2 * a2 * u[6] - a2^2 * u[3]
end

function noise!(du, u, p, t)
    du[1] = 0.
    du[2] = 0.
    du[3] = 0.
    du[4] = 0.
    du[5] = 0.
    du[6] = 0.
end

# Simulation settings
A0, a0, A1, a1, A2, a2 = 3.25, 100.0, 3.25, 100.0, 22., 50.0
C1, e0, r, v0 = 135.0, 5., 0.56, 6.0
C2, C3, C4 = 0.8 * C1, 0.25 * C1, 0.25 * C1

u0 = [0.39057, -204.28131, -3.88133, -0.11387, -7.68439, -253.32597]
tspan = (0.0, 10.0)
dts = 1e-2
dt = 1e-4

# noise settings 1
time_vals = 0:dts:tspan[2] |> collect 
noise_vals = rand(length(time_vals)) .* 320
I_ext = LinearInterpolation(time_vals, noise_vals, extrapolation_bc=Line())

p = [A0, a0, A1, a1, A2, a2, C1, C2, C3, C4, e0, r, v0, I_ext]

# Solve the ODE problem
prob = SDEProblem(jansen_rit!, noise!, u0, tspan, p)
sol = solve(prob, EM(), dt=dt, saveat=dts)

plot(sol[1, :])

ii) Simulation with MTK by having a placeholder for noise inside SDESystem, but then putting the true uniform noise (via NoiseGrid) only inside SDEProblem. This doesn’t work - the noise argument inside SDEProblem is completely disregarded, regardless of the values I put in.

Code and output for ii)

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using DifferentialEquations
using Plots

# Model
function sigmoid(V, e0, r, θ)
    return 2 * e0 / (1 + exp(r * (θ - V)))
end

params = @parameters A0 a0 A1 a1 A2 a2 C1 C2 C3 C4 e0 r v0
vars = @variables (u(t))[1:6]

eqs = [D(u[1]) ~ u[4],
       D(u[2]) ~ u[5],
       D(u[3]) ~ u[6],
       D(u[4]) ~ A0 * a0 * sigmoid(u[2] - u[3], e0, r, v0) - 2 * a0 * u[4] - a0^2 * u[1],
       D(u[5]) ~ A1 * a1 * (C2 * sigmoid(C1 * u[1], e0, r, v0)) - 2 * a1 * u[5] - a1^2 * u[2],
       D(u[6]) ~ A2 * a2 * (C4 * sigmoid(C3 * u[1], e0, r, v0)) - 2 * a2 * u[6] - a2^2 * u[3]
       ]

noiseeqs = [0., 0., 0., 0., 1. , 0.]

# Simulation settings
params_map = Dict(A0 => 3.25, a0 => 100.0, A1 => 3.25, a1 => 100.0, 
                  A2 => 22.0, a2 => 50.0, C1 => 135.0, C2 => 33.75,
                  C3 => 33.75, C4 => 108.0, e0 => 5.0, r => 0.56, v0 => 6.0)

inits = [0.39057, -204.28131, -3.88133, -0.11387, -7.68439, -253.32597]
tspan = (0.0, 10.0)
dts = 1e-2
dt = 1e-4

sys = structural_simplify(SDESystem(eqs, noiseeqs, t, vars, params; name=Symbol("JR")))

# Noise definition
time_vals = 0:dts:tspan[2] |> collect
noise_vals = rand(6, length(time_vals)) .* 320
noisegrid = NoiseGrid(time_vals, [noise_vals[:,i] for i in 1:size(noise_vals)[2]])

prob = SDEProblem(sys, inits, tspan, params_map; noise=noisegrid) 
sol = solve(prob, EM(), dt=dt, saveat=dts)
plot(sol[1, :])

iii) Simulation with MTK using Brownian distribution. It somehow works, but I cannot get the same output as with ground truth, it is only an approximation. If there is a better way, it would be great.

Code and output for iii)

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using DifferentialEquations
using Plots

function sigmoid(V, e0, r, θ)
    return 2 * e0 / (1 + exp(r * (θ - V)))
end

params = @parameters A0 a0 A1 a1 A2 a2 C1 C2 C3 C4 e0 r v0
vars = @variables (u(t))[1:6]

eqs = [D(u[1]) ~ u[4],
       D(u[2]) ~ u[5],
       D(u[3]) ~ u[6],
       D(u[4]) ~ A0 * a0 * sigmoid(u[2] - u[3], e0, r, v0) - 2 * a0 * u[4] - a0^2 * u[1],
       D(u[5]) ~ A1 * a1 * (C2 * sigmoid(C1 * u[1], e0, r, v0)) - 2 * a1 * u[5] - a1^2 * u[2],
       D(u[6]) ~ A2 * a2 * (C4 * sigmoid(C3 * u[1], e0, r, v0)) - 2 * a2 * u[6] - a2^2 * u[3]
       ]

noiseeqs = [0., 0., 0., 0., 102400. , 0.]

# Parameters
params_map = Dict(A0 => 3.25, a0 => 100.0, A1 => 3.25, a1 => 100.0, 
                  A2 => 22.0, a2 => 50.0, C1 => 135.0, C2 => 33.75,
                  C3 => 33.75, C4 => 108.0, e0 => 5.0, r => 0.56, v0 => 6.0)

inits = [0.39057, -204.28131, -3.88133, -0.11387, -7.68439, -253.32597]
tspan = (0.0, 10.0)
dts = 1e-2
dt = 1e-4

sys = structural_simplify(SDESystem(eqs, noiseeqs, t, vars, params; name=Symbol("JR")))
prob = SDEProblem(sys, inits, tspan, params_map) # noise=noisegrid) 
sol = solve(prob, EM(), dt=dt, saveat=dts)

plot(sol[1, :])

I anyone has some ideas and would be willing to help, would be great. Thank you!

NinaOmejc · February 25, 2025, 1:07pm

I tried another approach with symbolically registering rand() function and “artificially” adding noise equations to the ODESystem to create SDESystem, and seems this is a solution for my problem:

using ModelingToolkit: t_nounits as t, D_nounits as D
using DifferentialEquations
using Plots
using Interpolations

# Sigmoid function
function sigmoid(V, e0, r, θ)
    return 2 * e0 / (1 + exp(r * (θ - V)))
end

# Define a time-varying random function
function random_input(t)
    return 320.0 * rand()
end
@register_symbolic (random_input)(t)

params = @parameters A0=3.25 a0=100.0 A1=3.25 a1=100.0 A2=22.0 a2=50.0 C1=135.0 C2=0.8 * C1 C3=0.25 * C1 C4=0.25 * C1 e0=5.0 r=0.56 v0=6.0
vars = @variables (u(t))[1:6]

eqs = [D(u[1]) ~ u[4],
       D(u[2]) ~ u[5],
       D(u[3]) ~ u[6],
       D(u[4]) ~ A0 * a0 * sigmoid(u[2] - u[3], e0, r, v0) - 2 * a0 * u[4] - a0^2 * u[1],
       D(u[5]) ~ A1 * a1 * (random_input(t) + C2 * sigmoid(C1 * u[1], e0, r, v0)) - 2 * a1 * u[5] - a1^2 * u[2],
       D(u[6]) ~ A2 * a2 * (C4 * sigmoid(C3 * u[1], e0, r, v0)) - 2 * a2 * u[6] - a2^2 * u[3]
       ]

noiseeqs = [0., 0., 0., 0., 1., 0.]

# Parameters
params_map = Dict(A0 => 3.25, a0 => 100.0, A1 => 3.25, a1 => 100.0, 
                  A2 => 22.0, a2 => 50.0, C1 => 135.0, C2 => 0.8 * 135.0, 
                  C3 => 0.25 * 135.0, C4 => 0.25 * 135.0, e0 => 5.0, r => 0.56, v0 => 6.0)

inits = [0.39057, -204.28131, -3.88133, -0.11387, -7.68439, -253.32597]
tspan = (0.0, 10.0)
dts = 1e-2
dt = 1e-4

sys = structural_simplify(SDESystem(eqs, noiseeqs, t, vars, params; name=Symbol("JR")))
prob = SDEProblem(sys, inits, tspan, params_map) 
sol = solve(prob, EM(), dt=dt, saveat=dts)
plot(sol[1, :])

Mason · February 25, 2025, 1:33pm

Note that this is not a stochastic differential equation.

You don’t need to use an SDE here at all, you’re just solving an ODE with a pre-defined source term that happened to be randomly generated. You could have just written your DifferentialEquations.jl version as

using Plots, Interpolations, OrdinaryDiffEq

function sigmoid(V, e0, r, θ)
    return 2 * e0 / (1 + exp(r * (θ - V)))
end

function jansen_rit!(du, u, p, t)
    (;A0, a0, A1, a1, A2, a2, C1, C2, C3, C4, e0, r, v0, I_ext) = p
    
    du[1] = u[4]
    du[2] = u[5]
    du[3] = u[6]
    du[4] = A0 * a0 * sigmoid(u[2] - u[3], e0, r, v0) - 2 * a0 * u[4] - a0^2 * u[1]
    du[5] = A1 * a1 * (I_ext(t) + C2 * sigmoid(C1 * u[1], e0, r, v0)) - 2 * a1 * u[5] - a1^2 * u[2]
    du[6] = A2 * a2 * (C4 * sigmoid(C3 * u[1], e0, r, v0)) - 2 * a2 * u[6] - a2^2 * u[3]
end

A0, a0, A1, a1, A2, a2 = 3.25, 100.0, 3.25, 100.0, 22., 50.0
C1, e0, r, v0 = 135.0, 5., 0.56, 6.0
C2, C3, C4 = 0.8 * C1, 0.25 * C1, 0.25 * C1

u0 = [0.39057, -204.28131, -3.88133, -0.11387, -7.68439, -253.32597]
tspan = (0.0, 10.0)
dts = 1e-2
dt = 1e-4

# noise settings 1
time_vals = 0:dts:tspan[2] |> collect 
noise_vals = rand(length(time_vals)) .* 320
I_ext = LinearInterpolation(time_vals, noise_vals, extrapolation_bc=Line())

p = (;A0, a0, A1, a1, A2, a2, C1, C2, C3, C4, e0, r, v0, I_ext);
prob = ODEProblem(jansen_rit!, u0, tspan, p);
sol = solve(prob, Tsit5(), dt=dt, saveat=dts)
plot(sol[1, :])

NinaOmejc · February 25, 2025, 2:00pm

Thank you for the response!
Yes, I know, I just wanted to make it comparable between the DE and MTK versions that I was trying.

What is interesting though is that the same concept doesn’t work in the MTK world. When I tried to solve it via ODESystem, the ODE solver wasn’t able to solve it, there was always an error during simulation.

ChrisRackauckas · February 25, 2025, 2:05pm

Not quite. A randomly generated source does not necessarily have the regularity properties required for most ODE solvers to do this appropriately. It may be some kind of definitional / mathematical thing but if the noise source is not differentiable, for example Brownian motion is only Holder continuous for every epsilon < 0.5, and thus it does not match the assumptions required for ODE solvers to properly converge to the correct answer. If you put a linear interpolation between the points, you end up with something that is a piecewise smooth f, and thus Tsit5 ends up working on the problem for the specific pre-defined source term, though it can have some adaptivity issues with the distance between points is rather small because it needs an f with 5 derivatives defined and so as you sample more and more points this because less and less true, and the adaptivity effectively needs to take steps smaller than the noise sample size.

So this somewhat gets to a definitional thing. Yes, if you have a noisy process that you don’t sample very densely, and you only want to use this one noisy process, then you can put it into the ODE solver. But as you sample that more densely, it will not necessarily converge to the correct answer. So if you’re looking for something that truly simulates uniform noise, and believe that the model truly has this noise property, then this approach is not correct. But this is also making the implicit assumption that there is no maximum frequency to the randomness, that every possible speed of random fluctuations is possible up to infinity, and as dt->0 you see more and more of these frequencies, which is possibly a non-physical assumption.

But then for:

NinaOmejc:

I tried another approach with symbolically registering rand() function and “artificially” adding noise equations to the ODESystem to create SDESystem, and seems this is a solution for my problem:

Note that this kind of “accidentally” works. This will “work” for non-adaptive solvers only, because it’s not going to check the error which implicitly assumes a form of Brownian motion properties.

Effectively what you have here is a random ODE with a uniform noise process

Random ODEs are more general than Stochastic ODEs because Stochastic ODEs assume a Brownian process, or colored noise which is then represented via Brownian processes.

There is some theory on this. For the Euler method on random ODEs, you get convergence of O(dt^alpha) where alpha is the maximal Holder continuity of the process. For Brownian motion that is 1/2. For uniform noise, I don’t know? So this is only convergent if alpha > 0. But if you’re willing to just relax on the mathematics a bit, I’d say using the Random ODE formulation with uniform noise would be the most mathematically rigorous, where knowing the convergence rate (and if it converges) would require some mathematical proof.

That said, 1. ModelingToolkit cannot target the Random ODE form at this time. And 2…

This does not have the property that as dt->0 the random perturbations go to zero, so I cannot see this as a convergent process, I think alpha < x for every x can be proven because of the lack of this property. So I would assume this needs to be tempered like:

function random_input(t)
    return 320.0 * dt * rand()
end

to have any possibility of working. Basically, without this you don’t have the property that as dt->0 that the average change → 0, instead it goes to a constant, and so the average integral must go to infinity.

Note that the Random ODE formulation effectively does this automatically.

Mason · February 25, 2025, 2:12pm

Oh yeah, I’m not saying that that this is a good way to approach this problem. I was just saying that the problem as laid out wasn’t actually touching any of the SDE machinery because it didn’t actually have any noise (from the solver’s point of view).

It was just solving an ODE with jagged randomly generated source term. I definitely wasn’t recommending this approach though

NinaOmejc · February 25, 2025, 2:37pm

Thank you for the response!

Yeah, it felt like cheating a bit, I tried all combinations to make it work with MTK… but then if I stick to non-adaptive solvers like EM(), with a fixed dt, I can keep everything as it is? In a sense that here: 320.0 * dt * rand(), dt is just a constant and then it wouldn’t affect convergence? Or is this approach wrong anyhow and I shouldn’t continue doing this?

And secondly, this comment below was really insightful, I have to think about this part… thanks!

Also, actually the returned solution is RODESolution, which seems in order.

ChrisRackauckas · February 25, 2025, 3:15pm

Yeah that should be fine. I mean, to really prove it’s convergent you’d have to sit down with pen and paper and prove something quite difficult, but I would put this in the territory of “yeah that looks like it would be reasonably safe to assume it converges at some. speed, though that speed might be slow”.

SDEs in DifferentialEquations.jl use an RODESolution because there’s no difference in the solution of the two, both are stochastic. The difference really is that an SDEProblem has slightly more information for the solver than the RODEProblem, which the solver then specializes on. The solutions have nothing really to specialize on though.

NinaOmejc · February 25, 2025, 3:50pm

Great, thank you for all the input!

Topic		Replies	Views
Simulating SDEs with shared noise in `DifferentialEquations.jl` Modelling & Simulations	2	60	February 23, 2025
StochasticDiffEq with Acausal modeling Modelling & Simulations	4	61	November 3, 2024
Stochastic model building: assistance with the DiffEq ecosystem Numerics diffeq , sde	0	394	September 14, 2020
Is it possible to simulate SDDE with modeling toolkit? Modelling & Simulations modelingtoolkit	3	79	February 6, 2025
Combining ODE and SDE General Usage sde , differentialequation	4	516	March 29, 2022

ModelingToolkit and Noise Processes for uniform noise definition

Related topics