Add support for modeling with SI units

[chriscoey]

Something that has been mentioned before in passing and might be worth discussing down the track: allow modeling with SI units using https://github.com/ajkeller34/Unitful.jl. It may be fairly trivial, I'm not sure.

DifferentialEquations.jl seems to allow this, for example here.

GPKit is a geometric programming toolkit in Python and it supports specification of units and automatically does dimensional analysis to convert or verify units in expressions/constraints. From the GPKit docs, an example of a variable declaration:

# Declare a variable, z, with units of meters, and a description
z = Variable("z", "m", "A variable called z with units of meters")

and a simple model:

# consider a block with dimensions x, y, z less than 1
# constrain surface area less than 1.0 m^2
x = Variable("x", "m")
y = Variable("y", "m")
z = Variable("z", "m")
S = Variable("S", 1.0, "m^2")
c = (2*x*y + 2*x*z + 2*y*z <= S)

(There are other ideas we could borrow from GPKit to make JuMP even more user friendly. I like the little solution tables it can print after a solve.)

[blegat]

You can do this by extending the @variable macro as follows

using JuMP
using Unitful
m = Model()
function unitvar(vref::JuMP.VariableRef, u::Unitful.FreeUnits)
    JuMP.GenericAffExpr(0.0 * u, vref => 1.0 * u)
end
struct Unit
    u::Unitful.FreeUnits
end
struct UnitVariable <: JuMP.AbstractVariable
    v::JuMP.ScalarVariable
    u::Unitful.FreeUnits
end
function JuMP.buildvariable(_error::Function, info::JuMP.VariableInfo, u::Unit)
    UnitVariable(JuMP.ScalarVariable(info), u.u)
end
function JuMP.addvariable(m::Model, v::UnitVariable, name::String="")
    vref = JuMP.addvariable(m, v.v, name)
    # add bridge to `m.moibackend` if its does not support constraints with Unitful coefficients
    unitvar(vref, v.u)
end
@variable m v Unit(u"m")
@show v + 2v + 1u"m" # -> 3.0 m v + 1.0 m

To support adding constraints such as @constraint m 2v + 1 <= 2 , you can simply create a constraint bridge that transforms constraints with functions of Unitful coefficient type to constraints where the coefficients are the values of the unitful quantities (note that you need jump-dev/MathOptInterface.jl#397 to be is fixed for it to work with JuMP in non-Direct mode though).

Both the new constraint bridge and the @variable macro extension would fit nicely into a JuMP extension.
The extension would work nicely with other extension since you don't need to define a new model constructor or define a solvehook.

By the way this could be a nice example for the extensions doc don't you think ?

[odow]

By the way this could be a nice example for the extensions doc don't you think ?

Yes. This is ridiculous. Hooray for a good abstraction.

[IainNZ]

Could be fun to try with https://github.com/IainNZ/RationalSimplex.jl too.

[odow]

Pyomo has this: https://pyomo.readthedocs.io/en/stable/advanced_topics/units_container.html

It came up yesterday in a discussion with @zolanaj. It's particularly useful if you have a large model with lots of different physical units, and you want to avoid issues like adding meters and millimeters.

[trulsf]

I have been playing around with the above ideas to combine JuMP and Unitful. My primary goal was to use units in combination with the @variable and @constraint macros for linear constraints. I started out with the approach suggested by @blegat, but I found it hard to play nicely with the @rewrite macro. Instead I have tried an approach where constraints are handled similar to variables and are bundled with a unit. Being a fairly newcomer to both Julia and JuMP, I have not made a deep dive into the details of MutableArithmetics, and have only implemented limited support for the interface for expressions involving units.

I have made a proof-of-concept package: https://github.com/trulsf/UnitJuMP.jl

It runs on simple examples, but is in no way robust or complete. I would be grateful if anyone with a better grasp of JuMP and its inner workings, could comment upon the overall approach. Is it a viable way of approaching the problem? Or would it be better approach to change JuMP to better accomodate units in GenericAffExpr?

Here is a simple example of its usage:

using UnitJuMP, GLPK

m = Model(GLPK.Optimizer)

@variable(m, x >= 0, u"m/s")
@variable(m, y >= 0, u"ft/s")

max_speed = 60u"km/hr"

@constraint(m, x + y <= max_speed)
@constraint(m, x <= 0.5y)
obj = @objective(m, Max, x + y)

optimize!(m)

println("x = $(value(x))  y = $(value(y))")

#output x = 5.555555555555556 m s^-1  y = 36.45377661125693 ft s^-1

[odow]

Nice! I'm traveling for the next week, but I will take a look when I get a chance. If it's been 2 weeks and I haven't replied, please ping me as a reminder.

…

On Thu, 13 Jan 2022, 8:52 AM Truls Flatberg, ***@***.***> wrote: I have been playing around with the above ideas to combine JuMP and Unitful. My primary goal was to use units in combination with the @variable and @constraint macros for linear constraints. I started out with the approach suggested by @blegat <https://github.com/blegat>, but I found it hard to play nicely with the @rewrite macro. Instead I have tried an approach where constraints are handled similar to variables and are bundled with a unit. Being a fairly newcomer to both Julia and JuMP, I have not made a deep dive into the details of MutableArithmetics, and have only implemented limited support for the interface for expressions involving units. I have made a proof-of-concept package: https://github.com/trulsf/UnitJuMP.jl It runs on simple examples, but is in no way robust or complete. I would be grateful if anyone with a better grasp of JuMP and its inner workings, could comment upon the overall approach. Is it a viable way of approaching the problem? Or would it be better approach to change JuMP to better accomodate units in GenericAffExpr? Here is a simple example of its usage: using UnitJuMP, GLPK m = Model(GLPK.Optimizer) @variable(m, x >= 0, ***@***.***(m, y >= 0, u"ft/s") max_speed = 60u"km/hr" @constraint(m, x + y <= ***@***.***(m, x <= 0.5y) obj = @objective(m, Max, x + y) optimize!(m) println("x = $(value(x)) y = $(value(y))") #output x = 5.555555555555556 m s^-1 y = 36.45377661125693 ft s^-1 — Reply to this email directly, view it on GitHub <#1350 (comment)>, or unsubscribe <https://github.com/notifications/unsubscribe-auth/AB6MQJK4YPNENDDDFQNSZKTUVXLZLANCNFSM4FFTAHUQ> . Triage notifications on the go with GitHub Mobile for iOS <https://apps.apple.com/app/apple-store/id1477376905?ct=notification-email&mt=8&pt=524675> or Android <https://play.google.com/store/apps/details?id=com.github.android&referrer=utm_campaign%3Dnotification-email%26utm_medium%3Demail%26utm_source%3Dgithub>. You are receiving this because you commented.Message ID: <jump-dev/JuMP. ***@***.***>

[blegat]

I checked the implementation of UnitJuMP.jl, it looks very good.
About MutableArithmetics, implementing its interface is optional. It should work without it but be slower. So in the long term, it's best to make it work with MA but you don't need it if you just want to have things work. Once you implement: _MA.mutability(::Type{UnitAffExpr}) = _MA.IsMutable() then you'll indeed need to implement the interface to make it work. For the expressions to be used outside the macros, you'll still need to implement Base.:+, ... in addition to the MA interface.
That is quite a lot of work so I would recommend using GenericAffExpr{Quantity{U},VariableRef} instead of UnitAffExpr with unit U. Then you'll need to implement operations between GenericAffExpr and UnitVariableRef but operations with VariableRef might just work. I'd recommend try to make GenericAffExpr{Quantity{U},VariableRef} work without MutableArithmetics first, you can even test that additions etc... work outside the macros.
Let's me know if that's clear :)

One small comment on performance: the types of the fields of structure should be concrete, ScalarVariable, Units are not concrete types.

[trulsf]

Thanks for the good feedback.
My initial attempts was focused on using GenericAffExpr{Quantity{U},VariableRef}, but I quickly ran into problems with the operators (like Base.:+()) only being defined for arguments of the same type GenericAffExpr{C,V}. Since I want to allow for using different units (of same dimension) in the same constraint, I end up with errors in the code.
I guess for this to work, I will have to implement operators along the line

function Base.:+(
    lhs::GenericAffExpr{Quantity(U),V},
    rhs::GenericAffExpr{Quantity(W),V},
) where {U,W<:Units,V<:_JuMPTypes}

which should be doable, but some work. In addition, there are also problems with other operators like Base.:*(lhs::_Constant, rhs::_GenericAffOrQuadExpr) when involving units. So there seems to be quite a lot of work involved in going down this road as well. It may be that I am doing something wrong and there is a simple solution to this, so any tips?

[blegat]

I want to allow for using different units (of same dimension) in the same constraint

What do you mean ? How would it make sense to add affine expressions of different units ?
Yes, it might be some work but I would expect probably less methods needed than the approach with UnitAffExpr. The downside is that you need to play with the internal of GenericAffExpr.

[trulsf]

I want e..g. to combine variables of different units in the same constraint (there may e.g. be more natural units for different variables, typicially for our use is different scaling like kW and MW used for different technologies) . Returning to your first implementation above, I would want to support

@variable m v Unit(u"m")
@variable m w Unit(u"km")
@show v + w

that will now error due to v and w having different units. On the other hand, the following works (using the automatic promotion of Unitful)

@show v + 2u"km"

I can make some progress by allowing adding GenericAffExpr{C,V} with coefficients of different type for + and - operators

function Base.:+(
    lhs::GenericAffExpr{C,V},
    rhs::GenericAffExpr{D,V},
) where {C,D,V<:_JuMPTypes}

and continue with the same for the various add_to_expression!.
Even assuming using only the same unit I run into various problems with the @constraint-macro, e.g. in the operator_to_set function '<=` is mapped to a MOI.LessThan(0.0) causing problems when shifting a constant rhs with units. There are also problems when multiplying a quantity with a GenericAffExpr.

It is probably possible for one more familiar with the internals of JuMP and MathOptInterface to fix things, but I think it will be too time consuming for me at this point.

Having a separate UnitAffExpr to work with makes it easier for me to experiment with different solutions and I think I will continue a bit along that road just to learn a bit more before potentially diving into the details of JuMP/MathOptInterface/MutableArithmetics.

[blegat]

Sounds good, just wanted to give some advice but feel free to experiment in any direction!

[odow]

This is very cool! I like the use of a JuMP extension.

The printing of the units beside each constraint is very slick.

julia> m = Model(GLPK.Optimizer)
A JuMP Model
Feasibility problem with:
Variables: 0
Model mode: AUTOMATIC
CachingOptimizer state: EMPTY_OPTIMIZER
Solver name: GLPK

julia> @variable(m, x >= 0, u"m/s")
x m s⁻¹

julia> @variable(m, y >= 0, u"ft/s")
y ft s⁻¹

julia> max_speed = 60u"km/hr"
60 km hr⁻¹

julia> @constraint(m, x + y <= max_speed)
x + 0.3048 y ≤ 16.666666666666668 [m s⁻¹]

julia> @constraint(m, x <= 0.5y)
x - 0.1524 y ≤ 0.0 [m s⁻¹]

julia> obj = @objective(m, Max, x + y)
x + 0.3048 y m s⁻¹

julia> optimize!(m)

julia> println("x = $(value(x))  y = $(value(y))")
x = 5.555555555555556 m s⁻¹  y = 36.45377661125693 ft s⁻¹

I want e..g. to combine variables of different units in the same constraint

One way to simplify everything (and probably reduce the likelihood of bugs) is to not allow this. Instead, force the user to convert their units manually.

So instead of

julia> @expression(m, x + y)
x + 0.3048 y m s⁻¹

you would go

julia> y_expr = @expression(m, 1y)  # Could be avoided with some additional methods
y ft s⁻¹

julia> y_ms = UnitJuMP.convert(y_expr, u"m/s")
0.3048 y m s⁻¹

julia> @expression(m, x + y_ms)
x + 0.3048 y m s⁻¹

Automatic promotion seems useful, but it's likely to just complicate things. (It's too easy to accidentally type something weird and have the units all match, but differ from what was intended.) We should throw an error if the user provides mis-matching types. I think this is what the Pyomo package does as well.

[trulsf]

Thanks for the feedback!
In my view I think the introduction of units should both be a correctness and convenience measure. While the Pyomo solution seems to go for only the correctness issue by ensuring consistent units in a constraint, I kind of like the ability to combine different units in the same constraint with automatic conversions. Inconsistent use of dimensions would still produce an error.

julia> @variable(m, x>=0, u"m")
x m

julia> @variable(m, y>=0, u"s")
y s

julia> x + y
ERROR: DimensionError: m and s are not dimensionally compatible.

In some cases it is more convenient to have some variables in their natural units, e.g. in an energy model it can be that energy input is given in kcal while the energy requirements is based on mechanical work in kJ.

I have pushed an update that implements most operators for handling linear expressions with units, e.g.

julia> expr = x - 3u"km/hr"*y + 1.5u"km"
x - 0.8333333333333334 y + 1500 [m]

[odow]

Perhaps mixing units is okay. I wonder if we can use the tagging feature of constraints as well, which would allow:

max_speed = 1u"km/h"
@constraint(model, 2x + y <= max_speed, u"m/s")

That would give a little more control over what the final units are!

[trulsf]

There is already optional support for this through using unit as a keyword argument.

max_speed = 1u"km/h"
@constraint(model, 2x + y <= max_speed, unit=u"m/s")

I was not aware it was possible to manage it by tagging it along with the syntax you described. Any tips on how that could be done?

[odow]

You can pass positional arguments to build_constraint: https://jump.dev/JuMP.jl/stable/developers/extensions/#Build

Instead of the kwarg stuff you do here:
https://github.com/trulsf/UnitJuMP.jl/blob/17ff4b56c6f677d67a228f666d110c02575b070d/src/units.jl#L83-L87

[trulsf]

Thanks for the tip! Worked as a charm. I'll update with the new syntax option.

[mlubin]

I was convinced by @trulsf's JuMP-dev talk (https://www.youtube.com/watch?v=JQ6_LZfYRqg) that this functionality should ideally be incorporated into JuMP.

[odow]

this functionality should ideally be incorporated into JuMP.

An open question is whether we should vendor some form of UnitJuMP in as a JuMPExtension within JuMP, or whether we should extend VariableRef to include a unit field.

The former is probably easier, although it might make it hard for other JuMP extensions to use units.

The latter introduces a whole slew of problems, including the need to have a non-concrete field in VariableRef.

[odow]

This is potentially a good candidate for the new weakdeps extension packages that are landing in Julia 1.9: https://pkgdocs.julialang.org/dev/creating-packages/#Conditional-loading-of-code-in-packages-(Extensions)

We could move UnitJuMP to /ext/UnitfulExt.jl, and then if someone goes using JuMP, Unitful, it would load the extension.

[trulsf]

Seems like a good solution.

I am currently trying to add support for quadratic expressions into UnitJuMP and a rudimentary first version is available in the 'quadratic' branch. The approach used led to a lot of cases needed to be covered (quadratic expressions with units could arise from a lot of different combinations of variables, unit variables and quantities) and is not yet complete. I think there is room for a lot of improvement and any tips for possible refactoring/design changes are welcome.