Montecarlo DAE

Specific Domains Numerics
1

Hi everyone -

I am having a little problem defining a montecarlo DAE problem with varying parameters -
Any help would be appreciated.

Here is toy example, comparinng standard correct solution and montecarlo results.


#Satndard DAE

tspan = (0.0,1.0)
@everywhere function ode_f(du,u,p,t)
    du[1]=p[1]*u[1]
    du[2]=p[2]*u[2]
    du[3]=u[2]+u[1]-u[3]
end
mm=zeros(3,3)
mm[1,1]=mm[2,2]=1

system1=ODEFunction(ode_f;mass_matrix=mm)
prob1 = ODEProblem(system1,[10.0,10.0, 20.0],tspan,x[1])
sol1 = solve(prob1,Rodas5())

Note that sol1[3]=sol1[1]+sol1[2] as expected from the DAe deffinition

For my toy problem this also holds when running montecarlo mode


tspan = (0.0,1.0)
@everywhere function ode_f(du,u,p,t)
    du[1]=p[1]*u[1]
    du[2]=p[2]*u[2]
    du[3]=u[2]+u[1]-u[3]
end
mm=zeros(3,3)
mm[1,1]=mm[2,2]=1

system=ODEFunction(ode_f;mass_matrix=mm)
prob = ODEProblem(system,[10.0,10.0, 20.0],tspan)

@everywhere x=[[1.0,2.0],[3.0,4.0]]

@everywhere function prob_func(prob,i,repeat) # x is the new parameter vector
    p_new = x[i]
    x0 = prob.u0
    tspan = prob.tspan
    f = prob.f
    return ODEProblem(f,x0,tspan,p_new)
end
monte_prob=MonteCarloProblem(prob,prob_func=prob_func)
sol = solve(monte_prob, Rodas5(); num_monte=2)

However when I move to my complex system

@everywhere MCparams=fill(zeros(36),10)
@everywhere subparams=[]
for i in 1:size(PARAMS)[1]
     [MCparams[i][j-2]=PARAMS[i,j] for j in 3:(size(PARAMS)[2])]
end
MCparams[10]
print(MCparams[10])
print(IC)

[0.025, 2.52, 3.78, 7.0, 0.072, 0.183756, 334.0, 2.2e6, 6.22, 12.5, 77.4267, 0.551291, 0.635899, 0.1, 0.465208, 0.44, 65.0, 13.1, 0.558397, 21000.0, 2.5e7, 65.1, 212.1, 88.2, 55.93, 0.524, 0.25, 0.035, 0.00022, 25.0, 0.761, 241.126, 8.0, 81.6699, 0.314806, 8.0]
[457.516, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 4.0, 0.00419, 0.0, 0.0, 0.0, 0.00419, 0.0347196, 0.0, 10001.0, 0.131975, 1.49971]

Mass=zeros(21,21)
[Mass[i,i]=1.0 for i in 1:15]

ODESystem=ODEFunction(S.PKPD_Model;mass_matrix=Mass)
Problem1=ODEProblem(ODESystem,IC,tspan,MCparams[2])
solve(Problem1,Rodas5())

works without a glitch - I get a solution

but when trying to use montecarlo

Problem = ODEProblem(ODESystem,IC,tspan)

@everywhere function Problem_func(prob,i,repeat) # x is the new parameter vector


    p_new = MCparams[i]
    x0 = prob.u0
    tspan = prob.tspan
    f = prob.f

    return ODEProblem(f,x0,tspan,p_new)


end

monte_problem=MonteCarloProblem(Problem,prob_func=Problem_func)
Sol = solve(monte_problem, Rodas5(); num_monte=10)

I get inestability warning and don’t get a solution -

Why I dont understand is why I can solve the problem in a standard DAE format, but when I do montecarlo simulation with 10 param sets I get inestability warning.

PS- Sorry for the mistake I pressed to publish before time

2

What is the question here?

3

just edited -
I pressed publish by mistake - I don’t understand why montecarlo problems generate inestabilies in problems where standard ODEsolver solves without problem.

Thanks,

A

4

So I figure it out →
When defining MCparams

@everywhere MCparams=fill(zeros(36),10) # kept the params as 0 in the montecarlo simulation



# Correcting the for loop and defining PARAMS with @everywhere made the parameters available to montecarlo solver and solved inestability issue


@everywhere PARAMS=CSV.read("PARAMS.csv")

@everyhwhere for i in 1:size(PARAMS)[1]
     [MCparams[i][j-2]=PARAMS[i,j] for j in 3:(size(PARAMS)[2])]
end

It took me long to figure this, but hope this will help others :slight_smile:

A

5

PS: mark it as solved