# Montecarlo DAE

**URL:** <https://discourse.julialang.org/t/montecarlo-dae/25497>\
**Category:** Numerics\
**Created:** [June 20, 2019, 7:52pm UTC](https://discourse.julialang.org/t/montecarlo-dae/25497 "2019-06-20T19:52:50Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![alewolf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alewolf/32/9000_2.png) [@alewolf](https://discourse.julialang.org/u/alewolf)\
**Post date:** [June 20, 2019, 7:52pm UTC](https://discourse.julialang.org/t/montecarlo-dae/25497/1 "2019-06-20T19:52:50Z")

</div>

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.

```julia

#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

```julia

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

```julia
@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]

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

```

```julia
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

```julia
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

---

<div class="post-metadata">

**Author:** ![dpsanders](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dpsanders/32/3573_2.png) [@dpsanders](https://discourse.julialang.org/u/dpsanders)\
**Post date:** [June 20, 2019, 7:55pm UTC](https://discourse.julialang.org/t/montecarlo-dae/25497/2 "2019-06-20T19:55:31Z")

</div>

What is the question here?

---

<div class="post-metadata">

**Author:** ![alewolf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alewolf/32/9000_2.png) [@alewolf](https://discourse.julialang.org/u/alewolf)\
**Post date:** [June 20, 2019, 8:01pm UTC](https://discourse.julialang.org/t/montecarlo-dae/25497/3 "2019-06-20T20:01:51Z")

</div>

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

---

<div class="post-metadata">

**Author:** ![alewolf](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/alewolf/32/9000_2.png) [@alewolf](https://discourse.julialang.org/u/alewolf)\
**Post date:** [June 20, 2019, 8:12pm UTC](https://discourse.julialang.org/t/montecarlo-dae/25497/4 "2019-06-20T20:12:48Z")

</div>

So I figure it out →   
When defining MCparams

```julia
@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 🙂

A

---

<div class="post-metadata">

**Author:** ![mauro3](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mauro3/32/292_2.png) [@mauro3](https://discourse.julialang.org/u/mauro3)\
**Post date:** [June 21, 2019, 7:21am UTC](https://discourse.julialang.org/t/montecarlo-dae/25497/5 "2019-06-21T07:21:54Z")

</div>

PS: mark it as solved