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)?

5 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?

1 Like
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.

2 Likes
5

I am glad to see progress towards this, and I am willing to help.

Unitful is one of the things that drew me into Julia in the first place. That combined with this talk by @StefanKarpinski where he showed @ChrisRackauckas’s example of using Measurements + DifferentialEquaitons and it just working!

I personally believe that for many applications (mine included) using units is one of the best ways to “test” your code. You are effectively checking dimensional consistency on every expression you write! No more bugs because of accidentally using the continuous spectrum instead of the discrete because the units of the final solution would come out wrong.

So I was a bit disillusioned when I found out Unitful does not work as well with some of the other reasons I want to use Julia: automatic diff and the SciML ecosystem.

I now understand why this is. Linear algebra (etc) does not work with Unitful because of the units used as parameters in the type. But I also understand why that is the correct choice for Unitful (e.g. dispatch on dimensions).

I am excited about DynamicQuantities and am glad to see the seamless conversion between it and Unitful. Hopefully, this means that when SciML works with DynamicQuantities it would be straightforward to make it work with Unitful as well (by converting back and forth).

2 Likes
6

Making things work with DynamicQuantities.jl is difficult but possible (there’s some things ongoing in that direction). There’s certain things which just aren’t possible with Unitful without just stripping out the units in certain parts.

Many numerical methods are just not unit compatible though. Comparing it to Measurements.jl doesn’t even make sense. What are the units of a Jacobian inverse? The inverse and other matrix operations are linear combinations of all elements, and so if your u has different units then these operations are almost surely defining a Jacobian inverse that doesn’t make sense. And if linear algebraic operations don’t make sense, then Newton’s method and all sorts of numerical operations will fail. The problem there isn’t the units library, it’s that the algorithms themselves do not respect units. The only way to make them work is to strip the units, which is somewhat counter-productive to having units in the first place.

Also, units are just misplaced as number types. You want to unit check at the modeling level, not at the numerical implementation level. ModelingToolkit’s unit checking and verification interface is thus a much better fit because that canonicalizes before the numerical algorithm.

2 Likes
7

It’s really best to create a dimensionless version of your engineering problem as the first step after writing your equations. There are a number of reasons for this but one really important one is dimension reduction by 3 for most problems (Buckingham Pi theorem and mass,length, and time as the fundamental dimensions).

This needs to be taught more! People are often just unaware of this incredibly powerful mathematical technique.

1 Like
8

@ChrisRackauckas thanks for the quick and thoughtful reply. I am always amazed by this community’s responsiveness. :slight_smile:

Ah this all makes sense, thanks for explaining it this way! Sounds like it won’t be an easy problem for DynamicQuantities either, but glad there are efforts underway.

I’ll have to take a look at this to understand what you mean… I’ll keep learning more about the issue with units in computation. But this points me in the right direction.

9

Meant to link to the wiki page… But wasn’t quick enough

1 Like
10

Actually, Newton’s method does get physical dimensions/units right. The Adadelta paper has a short discussion (for unit less loss functions). Similarly, the units the ij-th entry of the Jacobian are just given as [\frac{\partial f_i}{\partial x_j}] = \frac{[f_i]}{[x_j]}.

I would view dimensions, e.g., time, distance etc, as akin to the “physical type” of a quantity. In turn, dimensional analysis could (should?) be done independently of the actual computation – like type checking or symbolic evaluation. Seems that Unitful and DynamicQuantities are missing some overloads to get matrix inversion etc working.

2 Likes
11

I just came across UnifulLinearAlgebra.jl. It creates a UnitfulMatrix type that does not have the units/dimensions in the type parameter. I am wondering what are the differences/advantages between this approach and that of DynamicQuantities. In any case using Unitful requires converting back and forth. E.g. these simple functions work for me:

inv(M::Matrix{<:Quantity}) = Matrix(inv(UnitfulMatrix(M)))
\(A::Matrix{<:Quantity}, B::Matrix{<:Quantity}) = Matrix(UnitfulMatrix(A) \ UnitfulMatrix(B))
2 Likes