Adding uncertainty to initial conditions and parameter value

[lakshaya17]

Hey,

I am trying to put uncertainty in both initial condition and parameter value for the bouncing ball example. It is giving me errors.

Is it possible for you to show how that can be done on the existing code for Bouncing ball example?

Thanks in advance.

[ChrisRackauckas]

Can you share your MWE? I don't think anything special would need to be done here.

[lakshaya17]

using Distributions

cor_dist = truncated(Normal(0.9, 0.02), 0.9 - 3 * 0.02, 1.0) # For parameter

cor_dist_1 = truncated(Normal(mean_x, std_x), mean_x- 3 std_x, mean_x+3 std_x) # For initial conditions for u[1] cor_dist_2 = truncated(Normal(mean_z, std_z), mean_z- 3 std_z, mean_z+3 std_z) # For initial conditions for u[3]

trajectories = 100

prob_func(prob, i, repeat) = remake(prob, u0= [rand(cor_dist_1), u[2], rand(cor_dist_2), u[4]],p = [p[1], rand(cor_dist)])

........ Will the distribution be like this ?

using SciMLExpectations gd = GenericDistribution(cor_dist, cor_dist_1, cor_dis_2) h(x, u, p) = [u[2];u[4]] [p[1]; x[1]] sm = SystemMap(prob, Tsit5(), callback = cbs) exprob = ExpectationProblem(sm, obs, h, gd; nout = 1) sol = solve(exprob, Koopman(), ireltol = 1e-5) sol.u

[ArnoStrouwen]

I think the following will get you started. It is a bit tricky how the h function works. I think in your code, the issue is that the first element h returns is a vector of only length 2.

#adapted from https://docs.sciml.ai/SciMLExpectations/stable/tutorials/optimization_under_uncertainty/
using DifferentialEquations
using SciMLExpectations
using Distributions

function ball!(du, u, p, t)
    du[1] = u[2]
    du[2] = 0.0
    du[3] = u[4]
    du[4] = -p[1]
end

ground_condition(u, t, integrator) = u[3]
ground_affect!(integrator) = integrator.u[4] = -integrator.p[2] * integrator.u[4]
ground_cb = ContinuousCallback(ground_condition, ground_affect!)

u0 = [0.0, 2.0, 50.0, 0.0]
tspan = (0.0, 50.0)
p = [9.807, 0.9]

prob = ODEProblem(ball!, u0, tspan, p)

stop_condition(u, t, integrator) = u[1] - 25.0
stop_cb = ContinuousCallback(stop_condition, terminate!)
cbs = CallbackSet(ground_cb, stop_cb)

sol = solve(prob, Tsit5(), callback = cbs)

obs(sol, p) = abs2(sol[3, end] - 25)

cor_dist = truncated(Normal(0.9, 0.02), 0.9 - 3 * 0.02, 1.0)
u0_randomness = truncated(Normal(50.0, 1.0), 45.0, 55.0) # THIS LINE IS NEW; will represent randomness in u[3]

gd = GenericDistribution(cor_dist,u0_randomness) # this now has 2 sources of randomness
h(x, u, p) = [u[1], u[2], x[2], u[4]], [p[1]; x[1]] # compare how this line is different from the original tutorials
# h is a function which has inputs, randomly drawn numbers from gd (x) and then the initial conditions (u) and parameters (p) from the original ODEProblem
# h outputs the initial conditions and parameters taking into account the randomness from x

sm = SystemMap(prob, Tsit5(), callback = cbs)
exprob = ExpectationProblem(sm, obs, h, gd; nout = 1)
sol = solve(exprob, Koopman(), ireltol = 1e-5)
sol.u