Distinguishing between unobserved parameters and control variables in DifferentialEquations

General Usage
, , ,
1

Hello, I’m fresh off the boat to using Julia (coming from Python and Stan), and I am looking for the preferred generic solution to specifying an ODE with both unknown parameters and control variables/system forcings.

As a simple case, I have a vector response function, Y = f(u,t,\theta;X), where X are control variables/system characteristics (like the volume of a tank, say, or inflow concentrations) potentially to be optimized to meet some constraints, and \theta are parameters to be learned (like reaction rates, settling velocities) from observed concentration data (u) via calibration (in my case, preferably Bayesian to support uncertainty quantification in the decision domain due to the unknown parameters).

I can easily make an ODEProblem work to separate user-specified “inputs” (i.e., X values) and parameters in the derivative function (thanks to the excellent SciMLTutorials), e.g.

    function derivatives!(du,u,p,t)
    theta = p[1:n]
    x1 = p[n+1]    # p[n+1] a tuple
    ...
    xn = p[last]    # p[last] a vector
    
    # now compute derivatives
    end

but this solution throws method errors I don’t yet understand how to interpret when passing the same problem to Turing for parameter estimation. Any tips/pointers would be much appreciated - especially so if there is a settled way of handling this pattern, generally. I am working on a simple example to see if the SciML ecosystem is right for my research focus (optimization of wasteload allocations in aquatic systems). So while this is currently just an ODEProblem, it will quickly become a PDE, SDE and perhaps UDE problem, as well.

2

For greater clarity, I did not include the errors I’m getting in the question above. But here is an example when passing a problem formed as I described to remake() as shown in the Turing examples:

MethodError: MethodError: no method matching remake(::ODEProblem{Vector{Float64}, Tuple{Float64, Float64}, true, Tuple{Float64, Float64, Float64, Float64, Float64, Float64, Float64, Int64, Float64, Float64, Float64, Float64, Float64, Float64, Float64, NTuple{9, Float64}, Vector{Float64}}, ODEFunction{true, typeof(derivatives!), LinearAlgebra.UniformScaling{Bool}, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, Nothing, typeof(SciMLBase.DEFAULT_OBSERVED), Nothing}, Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, SciMLBase.StandardODEProblem}, ::Vector{Float64})
Closest candidates are:
  remake(::ODEProblem; f, u0, tspan, p, kwargs...) at C:\Users\...
3

We probably can not reproduce your problem without a MWE.

4

Unfortunately I don’t have time right now to post a working example. I will ASAP, but I was hoping this was a general enough requirement that someone could just point me to the right docs or examples.

5

Here is a working example, but I am still unsure if this is the preferred way to pass unobserved parameters and control variables for Bayesian inference. The example is simply a tank of fixed volume, inflow and outflow rates, and constant loading at the inflow. A mass of organic nitrogen enters at the inflow and undergoes hydrolysis to ammonium as well as settling to the bottom of the tank, where it is considered to have left the system permanently. Two parameters control the hydrolysis rate, and one controls the settling rate. Several other parameters define the system characteristics as control variables, which should be supplied by the user (so I do not want to hard-code them). Inflow concentrations are the system forcing and will be used in downstream workflows as a decision variable to determine the maximum loading given some constraints (so they also cannot be hard-coded). Here is the code:

using Turing, Distributions, DifferentialEquations
using MCMCChains, Plots, StatsPlots
using Random
Random.seed!(14);

function derivs!(dc,c,p,t)
    # uncertain parameters (theta)
    khn20,khn_th,von = p[1]
    # system characteristics (simplified to constant values)
    A01, A12, elev, H, T, V, Qin, Qout = p[2]
    # system forcings
    c_inflow = p[3]
    # compute current reaction rates
    khn = khn20*khn_th^(T-20)
    # organic n terms
    on_hyd = khn*c[1]
    on_stl = von/V*A12*c[1]
    # compute net source terms for each state
    sink = -on_hyd-on_stl
    # compute derivatives
    dc = Qin/V*c_inflow[1]-Qout/V*c[1]+sink
end

# define inputs
A01 = 2.5e4 #surface area (m^2)
A12 = 2.5e4 #bottom area (m^2)
H = 10.0 #mean depth
V = H/3*(A01+A12+(A01*A12)^0.5)
elev = 100.0 #surface elevation (meters above sea level)
Qin = 7500.0 #inflow rate (m^3/day)
Qout = Qin
T = 21.0 #water temperature (Celsius)
w = [0.50] # inflow concentrations, g/m3 (cin)
x = [A01, A12, elev, H, T, V, Qin, Qout] # system characteristics
p = [0.10185547, #khn_20
    1.06767578, #khn_th
    0.06699219] #von
p = (p,x,w);

# simulate noisy data for parameter inference
c0 = [1.00]
tspan = (0.0,40.0)
prob = ODEProblem(derivs!,c0,tspan,p)
sol = solve(prob,Tsit5());
ts = collect(range(0.4,stop=38.0,length=20));
# interpolate solution at each time in t, add normal noise, convert to array
randomized = VectorOfArray([sol(ts[i]) .+ rand(Normal(0,0.05),1) for i in 1:length(ts)])
data = convert(Array,randomized);

Turing.setadbackend(:forwarddiff)

@model function fit(ts, data, prob, x, w)
    #priors
    khn20 ~ Uniform(0.05,0.15) #khn20
    khn_th ~ Uniform(1.0,1.1) #khn_th
    von ~ Uniform(0.05,0.25) #von
    σ ~ InverseGamma(2, 3) # standard deviation from functional mean

    p = [khn20,khn_th,von]
    # tuple of parameters, system constants, and forcings
    p = (p,x,w)
    prob = remake(prob, p=p)
    predicted = solve(prob,Tsit5(),saveat=ts)

    for i = 1:length(predicted)
        data[:,i] ~ MvNormal(predicted[i], σ)
    end
end

model = fit(ts, data, prob, x, w);
chain = sample(model, NUTS(0.95), 1000, progress=false)

Again, any constructive feedback is much appreciated.

1 Like
6

You took the W in MWE literally, I believe: the code is working on Julia 1.6.3.

7

It is - thanks for verifying. Again, I’m not looking so much for help de-bugging as I am in understanding the most robust/performant way to construct a problem of this basic type using Julia’s inherent design, given that this problem is a toy problem and the focus of my work is on much higher dimensional problems with many (potentially hundreds) more uncertain parameters. That said, it doesn’t seem to be of much interest to others, so I’ll just leave it open a few more days and close out if there are no further comments from the community. Thanks for your feedback, @goerch .