Montecarlo DAE

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

What is the question here?

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

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