Still looking for a way to use units meaningfully with SciML ecosystem

General Usage
,
1

I’ve spent frustrating hours now working on an optimization problem across what should be a relatively simple set of ODEs.

I am working in a field where some very commonly used units are frankly unpleasant (e.g. cm^2\cdot Torr \cdot hr/g , for one). I feel like working with these kinds of units is the whole point of packages like Unitful.jl. I would like to apply sensitivity analysis, as in the tutorial here, for parameter optimization.

As I learned here and read here and here, it’s not as simple as one would hope to do unit-aware calculations once there are solvers (and linear algebra, and automatic differentiation) involved. I understand there are clear technical difficulties, some of which could maybe be resolved by using a package like DynamicQuantities.jl where units are not encoded in the type.

I spent some time today trying to follow the ModelingToolkit.jl tutorial for units, but ultimately was frustrated by the fact that–because I want to work with quantities of compatible dimensions but different units–I still have to convert everything to the same base units before I start, which means that my parameter definitions can’t be nice and clean (and it’s not clear to me, from the documentation, what a clean approach would be if I need to do unit conversions before I define my @constants).

Looking for possible workarounds for making my ODE function compatible with ForwardDiff, I followed the old suggestion of wrapping the function in such a way that units are added at the entrance and stripped at the end. But as I learned the hard way, it is possible to break even that if you add plain numbers to quantities where units cancel (e.g. 1 + 5u"mK/K" breaks ForwardDiff in the following MWE, which I plan to submit somewhere as an issue):

using DifferentialEquations
using Unitful

function ode_system!(du, u, p, t)
    R0, τ, Tref = p
    
    T = u[1]*u"K"

    dTdt =  -T / (1 + R0*(Tref - T)) / τ

    du[1] = ustrip(u"K/hr", dTdt)
end

params = (3.0u"1/mK", 100.0u"s", 260u"K")
u0_ndim = [250.0, 100]
tspan = (0.0, 100.0)

prob = ODEProblem(ode_system!, u0_ndim, tspan, params)
sol = solve(prob, Rosenbrock23())

I want to write Julia code where these parameters are easy to read and compare to literature, and the answer that “you just have to strip all the units off yourself” feels really unsatisfying. When I am trying to persuade colleagues to use Julia, I feel like this would be a really nice cherry on top. Is there a way to make progress here? Are there roadblocks (with Unitful or DynamicQuantities or ForwardDiff) I could contribute to fixing (as an engineering PhD student determined to use Julia)?

3 Likes
2

I think in theory DynamicQuantities.jl could solve these issues (while Unitful in general cannot), but we just need a lot more overloads around. It would be good to document which functions exactly that it runs into.

3 Likes
3

So, to put this suggestion more concretely–you are suggesting trying out DynamicQuantities.jl with these problems (e.g. optimization, ODEs), trace the errors, and identify what functions need new methods to handle the Quantity type?

That is something I think I have the skills to start doing (if not always the time). The next question, then, is: should these errors get posted as GitHub issues in DynamicQuantities.jl, SciML libraries, or some mix of the two? Would there need to be an extension making DifferentialEquations.jl compatible with DynamicQuantities?

4

Yes exactly.

Either one, I’ll follow both.

I presume the differential equation solvers won’t need any changes, but things like linear solvers might run into issues and require extension libraries.

1 Like