# Problem setup for UODE testing

**URL:** <https://discourse.julialang.org/t/problem-setup-for-uode-testing/44422>\
**Category:** Machine Learning\
**Tags:** diffeq\
**Created:** [August 6, 2020, 2:21pm UTC](https://discourse.julialang.org/t/problem-setup-for-uode-testing/44422 "2020-08-06T14:21:51Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![Libbum](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/libbum/32/3935_2.png) [@Libbum](https://discourse.julialang.org/u/Libbum)\
**Post date:** [August 6, 2020, 2:21pm UTC](https://discourse.julialang.org/t/problem-setup-for-uode-testing/44422/1 "2020-08-06T14:21:51Z")

</div>

Hi all,  
I’m following these [Universal ODE](https://github.com/ChrisRackauckas/universal_differential_equations/tree/master/DelayLotkaVolterra) examples, trying to get one of my own systems to work. Having a good deal of issues getting things to solve & I hope someone can see the error I’m making.

Problem is that I see no errors turning up, some of the function calls are just spooling forever. I’ve turned on status bars on the ODE solver calls and I either get  
`ODE 0%| | ETA: N/A`  
or stuck somewhere around 25% with the ETA just continually climbing up, never going down.

My system:

\begin{align} \frac{dB}{dt} &= r\_B\frac{N}{N+H\_1}B-c\_B B^2 - p\_r \frac{PB^2}{H\_4^2+B^2}+i\_B,\\ \frac{dP}{dt} &= p\_r p\_e \frac{PB^2}{H\_4^2+B^2}\frac{V}{H\_2+V} - c\_P P^2 - U\_1(B,P,V) +i\_P, \\ \frac{dV}{dt} &= r\_V V - U\_2(B,P,V) - m\_V \frac{VB^2}{H\_3^2+B^2}. \end{align}

Where U\_{1,2} are the (arbitrary) portions of my system I’d like to try and recover. They are m\_p P and c\_V V^2 respectively in truth.

The `tspan` for my problem is large compared to the L-V example. The system runs on a time unit of days and the dynamics I’m interested in cover 30 years (10,950 days). A normal solve with `Tsit5()` is practically instantaneous over that range.

I’ve set up the training set to cover 3 years, with 30 data points over that span.

Current code:

```julia
using OrdinaryDiffEq
using ModelingToolkit
using DataDrivenDiffEq
using LinearAlgebra, DiffEqSensitivity, Optim
using DiffEqFlux, Flux
using JLD2

println("Generate data")
function lake!(du, u, p, t)
    B, P, V = u
    N, ib, ip, r, H₁, H₂, H₃, H₄, cb, cp, prmax, ce, mp, rv, cv, mv = p

    FR = B^2 / (B^2 + H₄^2)

    du[1] =
        dB = ib + r * (N / (N + H₁)) * B - cb * B^2 - prmax * FR * P
    du[2] =
        dP = ip + ce * prmax * FR * P * (V / (V + H₂)) - mp * P - cp * P^2
    du[3] =
        dV = rv * V - cv * V^2 - mv * (V * B^2 / (H₃^2 + B^2))
end

# Define the experiment
datasize = 30
tspan = (0.0f0, 1095.0f0) # Full data would be 10950: 30 years in days. Currently 3 years with 30 training points
tsteps = range(tspan[1], tspan[2], length = datasize)
u0 = Float32[60.036, 0.738, 11.654] # X1
p_ = Float32[
    0.7,
    2e-5,
    2e-5,
    7.5e-3,
    0.5,
    20,
    20,
    15,
    7.5e-5,
    2.75e-4,
    5e-2,
    0.14,
    2.25e-3,
    7e-3,
    6e-5,
    7e-3,
]
prob = ODEProblem(lake!, u0, tspan, p_)
solution = solve(prob, Vern7(), abstol = 1e-12, reltol = 1e-12, maxiters = 1e7, saveat = tsteps)

# Ideal data
X = Array(solution)
# Add noise to the ideal data
Xₙ = X + Float32(1e-3) * randn(eltype(X), size(X))

# Training implementation

# Input layer (BPV), hidden layer, ouptput layer (u1, u2)
ann_L = FastChain(FastDense(3, 32, tanh), FastDense(32, 32, tanh), FastDense(32, 2))
p = initial_params(ann_L)

function dudt_(u, p, t)
    B, P, V = u
    N, ib, ip, r, H₁, H₂, H₃, H₄, cb, cp, prmax, ce, mp, rv, cv, mv = p_
    # z should be u1, u2 correct?
    z = ann_L(u, p)

    FR = B^2 / (B^2 + H₄^2)

    [ib + r * (N / (N + H₁)) * B - cb * B^2 - prmax * FR * P,
     ip + ce * prmax * FR * P * (V / (V + H₂)) - z[1] - cp * P^2,
     rv * V - z[2] - mv * (V * B^2 / (H₃^2 + B^2))]
end

prob_nn = ODEProblem(dudt_, u0, tspan, p)
# contrete_solve has been depreciated apparently, although line below fails if I just use solve at present
# Depending on initialization, this line is the one that gets a portion of the way then just stops.
# Other times it solves almost instantly.
sol_nn = concrete_solve(prob_nn, Tsit5(), u0, p, saveat = solution.t)

function predict(θ)
    Array(concrete_solve(
        prob_nn,
        Vern7(),
        u0,
        θ,
        saveat = solution.t,
        abstol = 1e-6,
        reltol = 1e-6,
        maxiters = 1e7,
        sensealg = InterpolatingAdjoint(autojacvec = ReverseDiffVJP()),
    ))
end

function loss(θ)
    pred = predict(θ)
    sum(abs2, Xₙ .- pred), pred
end

# Test - this works without issue so long as sol_nn does.
loss(p)

losses = []

callback(θ, l, pred) = begin
    push!(losses, l)
     if !isempty(losses)
        println("Current loss after $(length(losses)) iterations: $(losses[end])")
     end
    false
end

# First train with ADAM for better convergence
res1 = DiffEqFlux.sciml_train(loss, p, ADAM(0.01), cb=callback, maxiters=200)

# Rest of the script is irrelevant for us.

```

And here is where I run into issues. `res1` never starts: the status bar just sits at `Training 0%| | ETA: N/A` which gives me little to go on in terms of troubleshooting - there’s a lot of packages I’m unfamiliar with and the problem could be in any of them, or more likely, my set-up.

The system doesn’t seem to be stiff. I’ve tried with `AutoTsit5(Rosenbrock23())` and `AutoVern7(Rodas5())` just in case - no difference.  
I’ve attempted to solve this on the gpu as well, by

```julia
u0 |> gpu
p |> gpu
p_ |> gpu
res1 = DiffEqFlux.sciml_train(loss, p, ADAM(0.01), cb=callback, maxiters=200)

```

but that also doesn’t do much of anything.

Can you see the issue 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:** [August 6, 2020, 3:59pm UTC](https://discourse.julialang.org/t/problem-setup-for-uode-testing/44422/2 "2020-08-06T15:59:19Z")

</div>

> [@Libbum](#):
>
> `tspan = (0.0f0, 1095.0f0) # Full data would be 10950: 30 years in days. Currently 3 years with 30 training points`

A random neural network probably doesn’t give you something stable over that huge of a time span. You’ll want to play with the initial parameters of the neural network, or use a strategy like demonstrated in [https://diffeqflux.sciml.ai/dev/examples/local\_minima/](https://diffeqflux.sciml.ai/dev/examples/local_minima/)

---

<div class="post-metadata">

**Author:** ![Libbum](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/libbum/32/3935_2.png) [@Libbum](https://discourse.julialang.org/u/Libbum)\
**Post date:** [August 7, 2020, 3:40pm UTC](https://discourse.julialang.org/t/problem-setup-for-uode-testing/44422/3 "2020-08-07T15:40:02Z")

</div>

Thanks for the suggestion Chris, although I don’t even seem to be getting to the point where the local minima strategy can be applied.

I’ve edited some input conditions to give at least some change over a smaller initial time window.

```julia
datasize = 30
tspan = (0.0f0, 10.0f0)                                                                                       
tsteps = range(tspan[1], tspan[2], length = datasize)                                                                                                                                          
u0 = Float32[90.984, 10.0, 80.816]                                                            
p_ = Float32[
    0.1,
    2e-5,
    2e-5,
    7.5e-3,
    0.5,
    20,
    20,
    15,
    7.5e-5,
    2.75e-4,
    5e-2,
    0.14,
    2.25e-3,
    7e-3,
    6e-5,
    7e-3,
]

```

The rest of the setup code is the same.  
If I run `l,pred = loss(p)` before the training, here’s an example output:

 ![Screenshot_20200807_173144](https://global.discourse-cdn.com/julialang/original/3X/8/c/8c59607d4244c2dab762b22eede78bc2be7efbbc.png)  
(lines are the initial solution, markers are predictions based on the initial random parameter values)

The result is still pretty flat, but I don’t see how this could be getting stuck in a local minima?

If I alter the callback function to plot updates:

```julia
callback(θ, l, pred) = begin
    push!(losses, l)
    display(l)
    plt = plot(tsteps[1:size(pred,2)], X[:,1:size(pred,2)]', legend = false)
    scatter!(plt, tsteps[1:size(pred,2)], pred')
    display(plot(plt))
    false
end

```

The `res1` call never manages to even display the loss value `l`, let alone update the plot. Something else must be amiss yeah?

**Edit:**  
Parameter set that isn’t as flat over the training window, still the same issue:

```julia
datasize = 30
tspan = (0.0f0, 10.0f0)
tsteps = range(tspan[1], tspan[2], length = datasize)
u0 = Float32[60.036, 0.738, 11.654]
p_ = Float32[
    0.7,
    2e-5,
    2e-5,
    7.5e-3,
    0.5,
    20,
    20,
    15,
    7.5e-3, # 7.5e-5,
    2.75e-6, # 2.75e-4,
    5e-2,
    14.0, # 0.14
    2.25e-4, #mp
    5e-2, #7e-3
    6e-7, #cv
    7e-3,
]

```

Initial parameters guess before `sciml_train`:

![Screenshot_20200807_182833](https://global.discourse-cdn.com/julialang/original/3X/2/5/256496e7647a77e7f5b9f330fe86d4ef71ce3380.png)

**Second Edit:** Actually, the parameter set above worked! It just took around 28 minutes to show the first step, then converged from there in 2 minutes.

So I can move on from here - thanks!

---

<div class="post-metadata">

**Author:** ![Libbum](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/libbum/32/3935_2.png) [@Libbum](https://discourse.julialang.org/u/Libbum)\
**Post date:** [August 11, 2020, 8:45am UTC](https://discourse.julialang.org/t/problem-setup-for-uode-testing/44422/4 "2020-08-11T08:45:52Z")

</div>

With the help of @ChrisRackauckas’ previous answer, I’ve managed to compress my timeline to a more functional range. Using the same equations as above, my new parameter set is

```julia
# Artificially collapse time scale to an appropreate range for the NN.
τ = 300
p_ = Float32[
    0.7,
    2e-5 * τ,
    2e-5 * τ,  
    7.5e-3 * τ,
    0.5,
    20,
    20,
    15,
    7.5e-5 * τ,
    2.75e-4 * τ,
    5e-2 * τ,   
    0.14,
    2.25e-3 * τ, #mp * τ, will need to unscale later
    7e-3 * τ,
    6e-5 * τ, #cv * τ, will need to unscale later
    7e-3 * τ,
]

```

Which enables my training data to be a more respectable

```julia
datasize = 30
tspan = (0.0f0, 3.0f0)
tsteps = range(tspan[1], tspan[2], length = datasize)

```

Thus, solving the NN portion is now a straightforward step, Yielding a good fit for my unknowns:

```julia
# Ideal data
# m_p P, c_v V^2
L̄ = [p_[13] * X[2,:]';p_[15] * (X[3,:] .^ 2)']
# Neural network guess
L̂ = ann_L(Xₙ, res2.minimizer)

```

![Screenshot_20200811_101307](https://global.discourse-cdn.com/julialang/original/3X/a/d/ad484a66e7ca711f1c727375401108304d5c6eca.png)

Now, I’m trying to recover m\_p P and c\_v V^2 using [SInDy](https://datadriven.sciml.ai/dev/quickstart/#Nonlinear-Systems-Sparse-Identification-of-Nonlinear-Dynamics-1), but don’t seem to be having much luck.

Here are the details you’ll need to follow along:  
Pretty simple `basis`:

```julia
@variables u[1:3]
h = Operation[u; u.^2; u.^3; sin.(u); cos.(u); 1]
basis = Basis(h, u)

```

My ‘real’ data:

```julia
# Ideal data, (B, P, V)
X = Array(solution)
# Add noise to the ideal data
Xₙ = X + Float32(1e-3) * randn(eltype(X), size(X))

```

> Float32[60.03501 58.813198 57.746056 56.817295 56.002934 55.284653 54.649925 54.08962 53.590942 53.149014 52.751167 52.39717 52.081867 51.79673 51.542362 51.31143 51.102936 50.915524 50.74641 50.58813 50.444656 50.31366 50.191795 50.080654 49.977543 49.879482 49.787895 49.698246 49.618042 49.544212; 0.7365888 0.73702586 0.7382404 0.73592335 0.73757845 0.73639107 0.73894846 0.73881716 0.73736215 0.7374129 0.7376731 0.7376398 0.73752075 0.73964417 0.74060017 0.7400208 0.74057287 0.74170345 0.74299425 0.74339664 0.7477569 0.7492879 0.7483499 0.7524558 0.7527696 0.7543931 0.7592906 0.7618538 0.76314175 0.76721525; 11.653926 11.654253 11.671077 11.692297 11.721293 11.757507 11.797537 11.843286 11.894215 11.949576 12.006296 12.067214 12.133057 12.1960335 12.265695 12.333571 12.405001 12.475596 12.548772 12.623673 12.696516 12.772361 12.844946 12.921922 12.996918 13.073674 13.146856 13.222248 13.296648 13.370284]

and the neural network guess `L̂ ` of my two unknowns:

> Float32[0.49594843 0.49605167 0.49617392 0.4962982 0.49644905 0.49659967 0.49678326 0.49696237 0.49714184 0.49734688 0.4975723 0.4978165 0.4980893 0.49843413 0.4988172 0.49922848 0.4997449 0.50036114 0.5011006 0.50196916 0.5030898 0.50432986 0.5056533 0.50740856 0.50924486 0.51142704 0.51392204 0.5166917 0.5196348 0.52292186; 2.4320884 2.4411485 2.4520636 2.464492 2.4785047 2.4943013 2.5114942 2.5303664 2.5509937 2.5731423 2.5966578 2.6217601 2.648635 2.6758575 2.7052617 2.7353282 2.7669077 2.7988124 2.8319352 2.8665643 2.9005435 2.9360812 2.9712234 3.0073447 3.043308 3.079943 3.1147962 3.150848 3.1859117 3.2197812]

Following the documentation of [DataDrivenDiffEq](https://datadriven.sciml.ai/dev/), I attempt

```julia
opt = SR3(1e-3, 1.0)
Ψ = SInDy(Xₙ[:, 2:end], L̂[:, 2:end], basis, opt=opt, maxiter=100000, normalize=true, denoise=true)

```

Which fails to identify a result for Equation 1:

```juliarepl
julia> println(Ψ)
Sparse Identification Result
No. of Parameters : 2
Active terms : 2
   Equation 1 : 0
   Equation 2 : 2
Overall error (L2-Norm) : 5.511965
   Equation 1 : 2.7072322
   Equation 2 : 2.8047326
AICC :
   Equation 1 : -3.8355591
   Equation 2 : 0.8347673

SR3{Float64,UnionAll}(0.3, 1.0, ProximalOperators.NormL1) converged after 30400 iterations.

julia> print_equations(Ψ)
1 dimensional basis in ["u₁", "u₂", "u₃"]
f_1 = u₃ ^ 2 * p₁ + u₃ ^ 3 * p₂

julia> ODESystem(Ψ)
ERROR: AssertionError: length(b) == length(variables(b))

```

The documentation is currently a little thin on what to do at this point. I’ve (blindly) tried different optimisers and parameters, altered the ranges in my input arrays, etc. Nothing sane in the results. In [this example](https://github.com/ChrisRackauckas/universal_differential_equations/blob/cdae0b41d5d7489f1e5e35104c828542d003c4dd/DelayLotkaVolterra/VolterraExp.jl#L164), Chris uses two extra parameters:

1. `λ = exp10.(-7:0.1:3)`, which verbatim doesn’t help my situation and I’m unsure what makes sense for my system anyhow.
2. `f_target(x, w) = iszero(x[1]) ? Inf : norm(w.*x, 2)`, but the `f_target` keyword seems to have been dropped in the current version v0.3.2 of the SInDy function. Should I be using an older version?

My questions now would be:

- How can I help the solver identify a 2 dimensional basis rather than dropping the first term in my result?
- How do I identify sane optimiser defaults for any given system? (should I not be using SR3 in this instance for example)?
- What range of \lambda makes sense for my problem and how do I identify that?
- Is `f_target` still relevant?

---

<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:** [August 11, 2020, 11:21am UTC](https://discourse.julialang.org/t/problem-setup-for-uode-testing/44422/5 "2020-08-11T11:21:02Z")

</div>

You might want to ask @Julius_Martensen how the latest Pareto stuff is going, because that should be a bit better at recovering a sparse form.

But that said, how good was the sparse form it gave you? It might be interesting to have a regularization that prefers some basis functions over others as well, biasing towards lower order or something.

---

<div class="post-metadata">

**Author:** ![Julius\_Martensen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/julius_martensen/32/17077_2.png) [@Julius\_Martensen](https://discourse.julialang.org/u/Julius_Martensen)\
**Post date:** [August 11, 2020, 11:50am UTC](https://discourse.julialang.org/t/problem-setup-for-uode-testing/44422/6 "2020-08-11T11:50:11Z")

</div>

Its getting there. We might fall back onto a mix of `0.3.1` and `0.3.2` given the [open PR](https://github.com/SciML/DataDrivenDiffEq.jl/pull/133) where the user is able to provide a feature space and an evaluation metric. I found some bugs / dumb mistakes today, so it might be better of in the future.

Choosing `λ` is indeed hard, given that no one knows a priori the appropriate size of the parameters. So this is try and error. Additionally, this is influenced by the `normalize` option, which scales the parameters. `SR3` tends to perform better than `STRRidge` if properly tuned ( you can also allow it to diverge more from the ridge regression ). I would try to perform a hyperparameter optimization around these and the optimizers.

---

<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:** [August 11, 2020, 11:58am UTC](https://discourse.julialang.org/t/problem-setup-for-uode-testing/44422/7 "2020-08-11T11:58:37Z")

</div>

It looks like it did grab a sparse form though, just a different one than intended, so it could just be unidentifiable given the data.

@Libbum you might want to sample more neural network input/outputs

---

<div class="post-metadata">

**Author:** ![Libbum](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/libbum/32/3935_2.png) [@Libbum](https://discourse.julialang.org/u/Libbum)\
**Post date:** [August 11, 2020, 12:03pm UTC](https://discourse.julialang.org/t/problem-setup-for-uode-testing/44422/8 "2020-08-11T12:03:41Z")

</div>

Thanks for the replies guys - will keep an eye on the open issues and play around with this system further.

This was an arbitrary test anyway and I’m not completely sure my system is sound in the first place!

Good suggestions as to where to go from here, and nice to know there’s nothing silly with my approach so far at least.

---

<div class="post-metadata">

**Author:** ![Libbum](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/libbum/32/3935_2.png) [@Libbum](https://discourse.julialang.org/u/Libbum)\
**Post date:** [August 11, 2020, 12:22pm UTC](https://discourse.julialang.org/t/problem-setup-for-uode-testing/44422/9 "2020-08-11T12:22:55Z")

</div>

> [@ChrisRackauckas](#):
>
> how good was the sparse form it gave you?

I couldn’t get that far to test due to this `length` assertion error when attempting to create the `ODESystem`…

---

<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:** [August 11, 2020, 12:42pm UTC](https://discourse.julialang.org/t/problem-setup-for-uode-testing/44422/10 "2020-08-11T12:42:47Z")

</div>

Hmm, the updated tested version should be here BTW: [https://github.com/SciML/DataDrivenDiffEq.jl/blob/master/test/applications/partial\_lotka\_volterra.jl](https://github.com/SciML/DataDrivenDiffEq.jl/blob/master/test/applications/partial_lotka_volterra.jl) . I think you might want to fallback to v0.3.2 until that open PR is finalized though.

---

<div class="post-metadata">

**Author:** ![Julius\_Martensen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/julius_martensen/32/17077_2.png) [@Julius\_Martensen](https://discourse.julialang.org/u/Julius_Martensen)\
**Post date:** [August 13, 2020, 8:32am UTC](https://discourse.julialang.org/t/problem-setup-for-uode-testing/44422/11 "2020-08-13T08:32:11Z")

</div>

The PR is stable and on master. I will publish a new version either today or tomorrow.

---

<div class="post-metadata">

**Author:** ![Libbum](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/libbum/32/3935_2.png) [@Libbum](https://discourse.julialang.org/u/Libbum)\
**Post date:** [August 13, 2020, 8:38am UTC](https://discourse.julialang.org/t/problem-setup-for-uode-testing/44422/12 "2020-08-13T08:38:33Z")

</div>

Awesome, thanks!