# DAEs and Turing

**URL:** <https://discourse.julialang.org/t/daes-and-turing/70651>\
**Category:** Modelling & Simulations\
**Tags:** turing\
**Created:** [October 30, 2021, 8:17am UTC](https://discourse.julialang.org/t/daes-and-turing/70651 "2021-10-30T08:17:04Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [October 30, 2021, 8:17am UTC](https://discourse.julialang.org/t/daes-and-turing/70651/1 "2021-10-30T08:17:04Z")

</div>

## 1. Model and DifferentialEquations

I’m trying to estimate initial values and parameters in a DAE using Turing. The structure of the DAE is as follows:

\frac{dx}{dt} = f(x,z,u;\theta) \\ 0 = g(x,z,u;\theta)

I have available known inputs u(t) and measured outputs y where, in general,

y = h(x,z,u;\theta)

The unknowns (“descriptor” in DAE speak) are x,z , or could also be denoted differential variables (x; 3 temperatures) and algebraic variables (z; 3 other temperatures). The DAE is index 1, which means that the system has 3 states. So I need 3 initial values to solve the DAE. One could think of the x variables as states.

The following Julia code works for simulating the DAE:

```julia
function system!(dx,x,p,t)
    Tr,Ts,TFe,Tac,Tad,Tah = x
    mr,ms,mFe = 9260.,6827.,71200.
    Rr,Rs = 0.127e-3,1.95e-6
    chpCu,chpFe,chpa,chpw = 0.385,0.465,1.15,4.18
    mda,mdw = 49.2,54.2
    QdFes,Qdfs = 212.,422.4
    UAr2d,UAs2Fe,UAFe2a,UAx,QdFes,Qdfs = p
    Nstw = UAx/chpw/mdw
    Nsta = UAx/chpa/mda
    Nstd = Nstw - Nsta
    # equations
    dx[1]= (1.1*Rr*If(t)^2-UAr2d*(Tr-Tad))/(mr*chpCu)
    dx[2]= (3*Rs*It(t)^2-UAs2Fe*(Ts-TFe))/(ms*chpCu)
    dx[3]= (UAs2Fe*(Ts-TFe)-UAFe2a*(TFe-Tah)+QdFes)/(mFe*chpFe)
    dx[4] = mda*chpa*(Tac-Tad)+UAr2d*(Tr-Tad)+Qdfs
    dx[5] = mda*chpa*(Tad-Tah)+UAFe2a*(TFe-Tah)
    dx[6] = Tac*(Nstw-Nsta*exp(-Nstd))-Nstd*Tah-Nsta*(1-exp(-Nstd))*Twc(t)
end

```

To set up the model and solve it, I do as follows:

```julia
x0 = [28., 28., 28., 14.,18.,22.]
tspan = (0.0,583*60.)
UAr2d,UAs2Fe,UAFe2a,UAx,QdFes,Qdfs = 2.7,20.,14.3,44.4,212.,422.4
p = [UAr2d,UAs2Fe,UAFe2a,UAx,QdFes,Qdfs]
probfunc = ODEFunction(system!, mass_matrix=Diagonal([1.,1.,1.,0.,0.,0.])) # DAEs to ODEs using Mass matrix reduction technique
prob = ODEProblem(probfunc, x0, tspan, p)
sol = solve(prob,saveat=60.)

```

This set-up works, and produces the solution.

_Question 1_: One thing I’m not 100% sure of is: what happens if there is inconsistency between the initail values of x and z (i.e., `x` in the Julia code)? Will the solver make adjustments so that the algebraic equations are satisfied at every time-step, including at initial time?

## 2. Data fitting with Turing

For data fitting, I have available outputs y consisting of 2 temperatures out of 3 from x (differential), and 2 temperatures out of 3 from z (algebraic).

_Question 2_: Is the above DAE implementation suitable for AD computations with Turing?

I’m assuming that all initial descriptor variables are unknown, i.e., I specify a prior distribution for every element in the descriptor. Of course, if x is known, then z is known from the algebraic equations, but I’m assuming it is possible to postulate uncertainty in both x and z…

_Question 3_: Does it make sense to postulate priors for both differential (x) and algebraic variables (z)? The answer to this question may depend on how the DAE solver handles satisfying the algebraic equations, I guess?

Here is how I set up the Turing model:

```julia
Turing.setadbackend(:forwarddiff) # Automatic Differentiation tool for Jacobian and Hessian
#
@model function fit_system(data,prob,tdata)
    # Measurement variance prior
    V ~ InverseGamma(2,3) # measurement variance https://en.wikipedia.org/wiki/Inverse-gamma_distribution
    # Initial differential variable priors
    Tr0 ~ truncated(Normal(30.,3),25,35) # Initial value prior
    Ts0 ~truncated(Normal(30.,3),25,35) # ""
    TFe0 ~truncated(Normal(30.,2),25,35) # ""
    # Initial algebraic variable priors
    Tac0 ~ truncated(Normal(12.,3),1.,30) # ""
    Tad0 ~truncated(Normal(18.,5),1,30) # ""
    Tah0 ~truncated(Normal(22.,3),1,60) # ""
    # Parameter priors
    UAr2d ~ truncated(Normal(2.7,1),0.2,5) # heat transfer prior
    UAs2Fe ~ truncated(Normal(20,5),1,50) # ""
    UAFe2a ~ truncated(Normal(15,2),0.5,40) # ""
    UAx ~ truncated(Normal(44,5),1,100) # ""
    QdFes ~ truncated(Normal(212.,40),20,400) # ""
    Qdfs ~ truncated(Normal(422.,20),200,500) # ""
    # Initial descriptor x0 and parameter p
    x0 = [Tr0,Ts0,TFe0,Tac0,Tad0,Tah0]
    p = [UAr2d,UAs2Fe,UAFe2a,UAx,QdFes,Qdfs]
    # Problem for differential equation -> solution
    prob_turing = remake(prob,u0=x0,p=p) # u0 and p should be of same type for autodiff
    sol = solve(prob_turing,saveat=60.,save_idxs=[2,3,4,6])
    # Output data model
    for i = 1:Nd
        data[:,i] ~ MvNormal(sol(tdata[i]), V*I)
    end
end
#
model = fit_system(data,prob,tdata); # model

```

Observe that I assume the same data uncertainty (variance `V`) in all 4 measured variables.

_Question 4_: Does this make sense? Perhaps… since all measurements are temperatures of similar size. But in general, this is a problem, I guess (unless data are normalized).

Finally, I run the data fitting:

```julia
Ni = 100
p_est = sample(model,NUTS(0.65),Ni)

```

_Question 5:_ I’ve seen `Ni` referred to as the number of iterations. Is this the same as the number of outcomes/realizations in the _posterior_ distribution?

I’ve also tried to run the code with multiple chains, either:

```julia
Nc = 3 # number of chains
p_est = sample(model,NUTS(.65),MCMCThreads(),Ni,Nc)

```

which seems to run… but forever… so I just have to stop it, or

```julia
p_est = sample(model,NUTS(.65),Ni,Nc)

```

which simply crashes…

_Question 6:_ I _assume_ that when I use multiple chains, this is just running the same sampling for `Nc` chains independently… and that the purpose is to see if I get the same posterior distributions for all chains?? Sort of to check that the found posterior distributions are representative of the true distributions?

Since MCMC is based on Monte Carlo methods, is should be possible to run samples in parallel – or – at least, chains in parallel. I assume that is what `MCMCThreads()` does. I have issued command `Threads.nthreads()` in Jupyter, and the answer is `8`.

_Question 7:_ Does this mean that I have `8` cores available in the current session? Or does it mean that my computer has `8` cores in total, and that I need to specify the number of cores I want to use when I start up Julia?

OK – sorry for many questions… when I run a single chain, typical results are that the posterior distributions are suspiciously narrow:

 ![image](https://global.discourse-cdn.com/julialang/original/3X/6/f/6f092c0faa8addbc7c14b994b24babd7d4d43b58.png)

Although my code runs (for a single chain), it is relatively slow (which is a relative concept). And if I re-run the notebook, many times I don’t get results at all.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [October 30, 2021, 12:57pm UTC](https://discourse.julialang.org/t/daes-and-turing/70651/2 "2021-10-30T12:57:33Z")

</div>

> [@BLI](#):
>
> _Question 1_ : One thing I’m not 100% sure of is: what happens if there is inconsistency between the initail values of x xx and z zz (i.e., `x` in the Julia code)? Will the solver make adjustments so that the algebraic equations are satisfied at every time-step, including at initial time?

Yes. [https://diffeq.sciml.ai/stable/solvers/dae\_solve/#Initialization-Schemes](https://diffeq.sciml.ai/stable/solvers/dae_solve/#Initialization-Schemes)

> [@BLI](#):
>
> _Question 2_ : Is the above DAE implementation suitable for AD computations with Turing?

Yes, both forward and reverse.

> [@BLI](#):
>
> _Question 3_ : Does it make sense to postulate priors for both differential (x xx ) and algebraic variables (z zz )? The answer to this question may depend on how the DAE solver handles satisfying the algebraic equations, I guess?

Yes, though initialization may cause it to be a little odd.

> [@BLI](#):
>
> _Question 4_ : Does this make sense? Perhaps… since all measurements are temperatures of similar size. But in general, this is a problem, I guess (unless data are normalized).

Looks fine.

> [@BLI](#):
>
> _Question 6:_ I _assume_ that when I use multiple chains, this is just running the same sampling for `Nc` chains independently… and that the purpose is to see if I get the same posterior distributions for all chains?? Sort of to check that the found posterior distributions are representative of the true distributions?

Different chains can hit different optima so this can help explore the space. And it can be done in parallel.

> [@BLI](#):
>
> _Question 7:_ Does this mean that I have `8` cores available in the current session? Or does it mean that my computer has `8` cores in total, and that I need to specify the number of cores I want to use when I start up Julia?

Yes.

This is a ton of unrelated questions. You’ll get better answers if you don’t put these in a discussion about DAEs, since I might be the only person to dig through a very long DAE post.

> [@BLI](#):
>
> OK – sorry for many questions… when I run a single chain, typical results are that the posterior distributions are suspiciously narrow:

Run it longer?

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [October 30, 2021, 2:43pm UTC](https://discourse.julialang.org/t/daes-and-turing/70651/3 "2021-10-30T14:43:23Z")

</div>

> [@ChrisRackauckas](#):
>
> since I might be the only person to dig through a very long DAE post.

Nope. But I think so far you have good answers.

The suspiciously narrow output could be a number of things, one might be that it is stuck in initial location. Another might be that the model fits the data extremely well. Are you using real data or simulated? If simulated make sure to add some nontrivial measurement noise.

Try running an optimize to get initial parameters for turing. The docs say to use init\_theta but they are out of date, use keyword init\_params to initialize the sample method.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [October 30, 2021, 2:53pm UTC](https://discourse.julialang.org/t/daes-and-turing/70651/4 "2021-10-30T14:53:46Z")

</div>

I use real data. So I’m a little suspicious of the narrow posteriors. Running data retrodiction, I cannot distinguish between the different draws from the posterior distribution.

What I find strange is that when I get _some_ spread in the posterior distributions, and re-run the code, I may get virtually no variation or results that look meaningsless. Sometimes when I run it with `Ni=1000` iterations (or: rather particles?? or outcomes for the distribution??.. my unanswered Question 5…), I get no reasonable results – other times I get some distribution – but often not big variation.

My Question 7 was poorly phrased… because answer “yes” could mean two things…

Since my code _crashed_ when I tried with 3 chains _without_`MCMCThreads()` – does this mean that I can only run more than one chain when I spread computations on several cores?

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [October 30, 2021, 3:05pm UTC](https://discourse.julialang.org/t/daes-and-turing/70651/5 "2021-10-30T15:05:35Z")

</div>

Start with getting a single chain to run and give reasonable results. When you specify Ni this is the length of the chain, often called the number of iterations or draws or samples. They are samples from the posterior distribution.

> [@BLI](#):
>
> cannot distinguish between the different draws from the posterior distribution.

Most likely it is stuck at the initial position. Likely this position is weird. I suggest to run the optimization step to get a not too weird initial position. After that I sometimes run sample with the MH() method to fuzz that initial position before handing to NUTS

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [October 31, 2021, 9:55am UTC](https://discourse.julialang.org/t/daes-and-turing/70651/6 "2021-10-31T09:55:47Z")

</div>

OK – I’ve read parts of the documentation, and I _think_ I understand parts of your procedure:

1. Run (mode) optimization using Turing, either MLE or MAP, to find initial estimates for the model.
2. Set the estimate from the optimization step as initial value in the `sample` method using kwarg `init_params`.

I can also consider an intermediate step between 1 and 2 above using the `MH()` sampler to fuzz the initial position before handing to NUTS.

This `MH()` step is not clear to me.

- What is the outcome of the `MH()` step? Point estimates that I can use with `init_params` in the `NUTS()` sampler?
- Or is the outcome of the `HM()` step a _distribution_ that I can use as prior in the `NUTS()` step?

One more general question (related to the second bullet point above)… Suppose I run Turing _recursively_, and each time I run the sampler, I want to use the posterior distribution from the _previous_ step as prior distribution in the current step… can this be done?

- If `x_d` is a vector of samples from the distribution of `x`, can I then in the Turing `model` part set… `x ~ x_d` – instead of, say, `x ~ Normal(x_, s)`
- Or would I have to use `Distributions` to _fit_ the distribution sample `x_d` to some theoretical distribution?
- And – how could I take correlations between parameters into account?

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [October 31, 2021, 12:54pm UTC](https://discourse.julialang.org/t/daes-and-turing/70651/7 "2021-10-31T12:54:54Z")

</div>

> [@BLI](#):
>
> What is the outcome of the `MH()` step?

It is samples from the posterior, so you would use the final sample as input to init\_params

> [@BLI](#):
>
> samples from the distribution of `x` , can I then in the Turing `model` part set… `x ~ x_d` – instead of, say, `x ~ Normal(x_, s)`

No, the model uses a distribution represented as a density function, whereas x\_d is samples from the distribution, not a density.

When you try to fit a distribution to the samples, inevitably you lose the dependency structure.

---

<div class="post-metadata">

**Author:** ![BLI](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/bli/32/37206_2.png) [@BLI](https://discourse.julialang.org/u/BLI)\
**Post date:** [October 31, 2021, 2:52pm UTC](https://discourse.julialang.org/t/daes-and-turing/70651/8 "2021-10-31T14:52:02Z")

</div>

Great – then I think I know how to do it.

[I’m a little surprised that I should use the \*final\* sample of the `MH()` samples… that indicates that the samples get progressively better – I thought they were just random samples from the distribution.]

Using `NUTS()`, I observe that with a _low_ number of samples (e.g., 10), the distribution gets “wide”, and as I increase the number of samples, it gets more and more narrow… (with 100 samples, the posterior distribution is narrow, and with 1000 samples, values appear to be constant.) This would be consistent with that as the number of samples get larger and larger, the algorithm gets more and more stuck in a poor initial conditions. But that would imply that the step size in trying to escape a poor initial guess varies inversely with the number of samples??

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [October 31, 2021, 3:24pm UTC](https://discourse.julialang.org/t/daes-and-turing/70651/9 "2021-10-31T15:24:55Z")

</div>

> [@BLI](#):
>
> [I’m a little surprised that I should use the \*final\* sample of the `MH()` samples… that indicates that the samples get progressively better – I thought they were just random samples from the distribution.]

Well, the point here is that you’re trying to find a good starting point. Once you’ve got sampling going right, any random sample is kind of as good as another. But during the initial stages you’re starting _well_ out of equilibrium, and you’re trying to get your point into the high probability set, so generally the last sample will be closest to equilibrium.

Also, in general for MCMC each sample has some correlation with the previous sample, so it takes several steps to get a sample that is “as good as independent”. Hence the output of the ESS (effective sample size) calculation in the summaries isn’t the same as the total number of samples (usually lower, though if the samples are anti-correlated it can be higher)

---

<div class="post-metadata">

**Author:** ![dlakelan](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dlakelan/32/8491_2.png) [@dlakelan](https://discourse.julialang.org/u/dlakelan)\
**Post date:** [October 31, 2021, 4:11pm UTC](https://discourse.julialang.org/t/daes-and-turing/70651/10 "2021-10-31T16:11:28Z")

</div>

> [@BLI](#):
>
> with a _low_ number of samples (e.g., 10), the distribution gets “wide”, and as I increase the number of samples, it gets more and more narrow…

if you start far from the high probability set, you will generally converge towards it. so since you’re moving in the same direction all the time (towards the set) it will appear to be very “wide”. Once you find the high probability set you’ll just move around within it… relative to the convergence period, it could be quite narrow. The problem comes if you have NUTS learn the step size while it’s converging, but then once it’s converged it needs a smaller step size. I usually use NUTS(1000,.75) to tell it to do 1000 initial “learning” steps before doing the real sampling (0.75 = the frequency of acceptance, this lets NUTS be a little looser with its hamiltonian solver).

Then I usually have it sample between 1000 and 5000 “real” samples.