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.

4 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.

3 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.

2 Likes
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
12

This statement has bothered me ever since I read it; after a quick read of George Hart’s “Multidimensonal Analysis”, I have been persuaded that a dimensionful Jacobian can be well defined and should also have a well-defined inverse. Any matrix that represents a set of dimensionally-consistent equations can be given a well-defined structure of its dimensions, as noted here and implemented by one person here. Certainly it is true that many algorithms may not be dimensionally consistent as currently implemented, but I am optimistic that many of these are not unresolvable issues, and maybe with wrappers that carry around the dimensional parts, we can perhaps even still use BLAS, etc.

On the other hand, discussion in Quantity currently incompatible with QuadGK · Issue #40 · SymbolicML/DynamicQuantities.jl · GitHub makes me inclined to believe that encoding dimensions in types hews closer to a mathematically rigorous definition of dimension.

13

Not everything that is in principle feasible is also sensible.

Unit-checking should be a one time activity when checking that a model is dimensionally consistent; there isn’t much sense in keeping to check this at the level of every computation (which is what Chris pointed out already).

And yes, the non-dimensionalisation of models for dimensionality reduction is a diminishing art - the number of times I see people at conferences presenting results where they have squandered numerical effort by simulating completely redundant parameter combinations is mindboggling (not that industry is any better…)

The appropriate scaling of numerical problems is another thing where many people just throw out efficiency while they are busy micro-optimizing code :slight_smile:

14

Unit-checking should be a one time activity

I think by implementing it through the type-system (which I believe is how it’s usually done) this is the case. The compiler should be able to type check the computations “once” and then just run the computation as if it had no units.

1 Like
15

It’s a fair point that encouraging proper nondimensionalization does reduce the number of parameters you have to work with. But in the case I am looking at, we have something like 20 different parameters and physical constants, and a reduction to 16 or so parameters doesn’t help all that much (and may not even produce numbers of order 1).

My main motivation for wanting units carried through lies less in making sure equations are consistent and more in making sure the right units are attached to my final answer. Yes, a single check at the start can answer the consistency question. But it doesn’t prevent me from forgetting whether I put time in seconds or hours at the end, which admittedly is a trivial thing to catch but exemplifies a broader class of the question “what units did we put this in?” A Google search is enough to turn up examples like the Mars orbiter where losing track of which units were used resulted in real, hard-to-diagnose problems.

1 Like
16

My favourite unit conversion problem that could have been disastrous but had a happy ending was the 767 that landed at Gimli.

Personally I like converting to SI units for internal calculations and database use so I know what the units are.

17

I understand the point about controlling units to avoid silly mistakes, but the way I handle this is to have an input checking function.

Essentially like this:

  • It exclusively takes unitful inputs,
  • converts the given units to the desired standard
  • then strips them
  • and then nondimensionalizes and scales.

After running the simulation one just adds units again. Still in full control. :slight_smile:

1 Like
18

That’s precisely the step that makes me nervous. I use a similar kind of workflow, but I’m always worried that I’ll tack the wrong units on because I changed something somewhere else.

3 Likes