Skip to content

Commit 8d247be

Browse files
author
Lucas Morton
committed
Switched to returning units instead of quantities. It makes validation a bit trickier but is more correct. Added screening for unsupported units. Fixed bug with 2nd derivatives. Improved & checked docs. Switched from get_units to get_unit for consistency.
1 parent f4f8302 commit 8d247be

File tree

3 files changed

+215
-77
lines changed

3 files changed

+215
-77
lines changed

docs/src/basics/Validation.md

Lines changed: 116 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -1,18 +1,125 @@
1-
# Model Validation and Units
1+
# [Model Validation and Units](@id units)
22

3-
ModelingToolkit.jl provides extensive functionality for model validation
4-
and unit checking. This is done by providing metadata to the variable
5-
types and then running the validation functions which identify malformed
6-
systems and non-physical equations.
3+
ModelingToolkit.jl provides extensive functionality for model validation and unit checking. This is done by providing metadata to the variable types and then running the validation functions which identify malformed systems and non-physical equations. This approach provides high performance and compatibility with numerical solvers.
74

8-
## Consistency Checking
5+
## Assigning Units
96

10-
```@docs
11-
check_consistency
7+
Units may assigned with the following syntax.
8+
9+
```julia
10+
using ModelingToolkit, Unitful
11+
@variables t [unit = u"s"] x(t) [unit = u"m"] g(t) w(t) [unit = "Hz"]
12+
#Or,
13+
@variables(t, [unit = u"s"], x(t), [unit = u"m"], g(t), w(t), [unit = "Hz"])
14+
#Or,
15+
@variables(begin
16+
t, [unit = u"s"],
17+
x(t), [unit = u"m"],
18+
g(t),
19+
w(t), [unit = "Hz"]
20+
end)
1221
```
1322

14-
## Unit and Type Validation
23+
Do not use `quantities` such as `1u"s"` or `1/u"s"` or `u"1/s"` as these will result in errors; instead use `u"s"` or `u"s^1"`.
24+
25+
## Unit Validation & Inspection
26+
27+
Unit validation of equations happens automatically when creating a system. However, for debugging purposes one may wish to validate the equations directly using `validate`.
1528

1629
```@docs
1730
ModelingToolkit.validate
1831
```
32+
33+
Inside, `validate` uses `get_unit`, which may be directly applied to any term. Note that `validate` will not throw an error in the event of incompatible units, but `get_unit` will. If you would rather receive a warning instead of an error, use `safe_get_unit` which will yield `nothing` in the event of an error. Unit agreement is tested with `ModelingToolkit.equivalent(u1,u2)`.
34+
35+
36+
```@docs
37+
ModelingToolkit.get_unit
38+
```
39+
40+
Example usage below. Note that `ModelingToolkit` does not force unit conversions to preferred units in the event of nonstandard combinations -- it merely checks that the equations are consistent.
41+
42+
```julia
43+
using ModelingToolkit, Unitful
44+
@parameters τ [unit = u"ms"]
45+
@variables t [unit = u"ms"] E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
46+
D = Differential(t)
47+
eqs = eqs = [D(E) ~ P - E/τ,
48+
0 ~ P ]
49+
ModelingToolkit.validate(eqs) #Returns true
50+
ModelingToolkit.validate(eqs[1]) #Returns true
51+
ModelingToolkit.get_unit(eqs[1].rhs) #Returns u"kJ ms^-1"
52+
```
53+
54+
An example of an inconsistent system: at present, `ModelingToolkit` requires that the units of all terms in an equation or sum to be equal-valued (`ModelingToolkit.equivalent(u1,u2)`), rather that simply dimensionally consistent. In the future, the validation stage may be upgraded to support the insertion of conversion factors into the equations.
55+
56+
```julia
57+
using ModelingToolkit, Unitful
58+
@parameters τ [unit = u"ms"]
59+
@variables t [unit = u"ms"] E(t) [unit = u"J"] P(t) [unit = u"MW"]
60+
D = Differential(t)
61+
eqs = eqs = [D(E) ~ P - E/τ,
62+
0 ~ P ]
63+
ModelingToolkit.validate(eqs) #Returns false while displaying a warning message
64+
```
65+
66+
## `Unitful` Literals & User-defined functions
67+
68+
In order for a function to work correctly during both validation & execution, the function must be unit-agnostic. That is, no unitful literals may be used. Any unitful quantity must either be a `parameter` or `variable`. For example, these equations will not validate successfully.
69+
70+
```julia
71+
using ModelingToolkit, Unitful
72+
@variables t [unit = u"ms"] E(t) [unit = u"J"] P(t) [unit = u"MW"]
73+
D = Differential(t)
74+
eqs = [D(E) ~ P - E/1u"ms" ]
75+
ModelingToolkit.validate(eqs) #Returns false while displaying a warning message
76+
77+
myfunc(E) = E/1u"ms"
78+
eqs = [D(E) ~ P - myfunc(E) ]
79+
ModelingToolkit.validate(eqs) #Returns false while displaying a warning message
80+
```
81+
82+
Instead, they should be parameterized:
83+
84+
```julia
85+
using ModelingToolkit, Unitful
86+
@parameters τ [unit = u"ms"]
87+
@variables t [unit = u"ms"] E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
88+
D = Differential(t)
89+
eqs = [D(E) ~ P - E/τ]
90+
ModelingToolkit.validate(eqs) #Returns true
91+
92+
myfunc(E,τ) = E/τ
93+
eqs = [D(E) ~ P - myfunc(E,τ)]
94+
ModelingToolkit.validate(eqs) #Returns true
95+
```
96+
97+
It is recommended *not* to circumvent unit validation by specializing user-defined functions on `Unitful` arguments vs. `Numbers`. This both fails to take advantage of `validate` for ensuring correctness, and may cause in errors in the
98+
future when `ModelingToolkit` is extended to support eliminating `Unitful` literals from functions.
99+
100+
## Other Restrictions
101+
102+
`Unitful` provides non-scalar units such as `dBm`, `°C`, etc. At this time, `ModelingToolkit` only supports scalar quantities. Additionally, angular degrees (`°`) are not supported because trigonometric functions will treat plain numerical values as radians, which would lead systems validated using degrees to behave erroneously when being solved.
103+
104+
## Troubleshooting & Gotchas
105+
106+
If a system fails to validate due to unit issues, at least one warning message will appear, including a line number as well as the unit types and expressions that were in conflict. Some system constructors re-order equations before the unit checking can be done, in which case the equation numbers may be inaccurate. The printed expression that the problem resides in is always correctly shown.
107+
108+
Symbolic exponents for unitful variables *are* supported (ex: `P^γ` in thermodynamics). However, this means that `ModelingToolkit` cannot reduce such expressions to `Unitful.Unitlike` subtypes at validation time because the exponent value is not available. In this case `ModelingToolkit.get_unit` is type-unstable, yielding a symbolic result, which can still be checked for symbolic equality with `ModelingToolkit.equivalent`.
109+
110+
## Parameter & Initial Condition Values
111+
112+
Parameter and initial condition values are supplied to problem constructors as plain numbers, with the understanding that they have been converted to the appropriate units. This is done for simplicity of interfacing with optimization solvers. Some helper function for dealing with value maps:
113+
114+
```julia
115+
remove_units(p::Dict) = Dict(k => Unitful.ustrip(ModelingToolkit.get_unit(k),v) for (k,v) in p)
116+
add_units(p::Dict) = Dict(k => v*ModelingToolkit.get_unit(k) for (k,v) in p)
117+
```
118+
119+
Recommended usage:
120+
121+
```julia
122+
pars = @parameters τ [unit = u"ms"]
123+
p = Dict=> 1u"ms")
124+
ODEProblem(sys,remove_units(u0),tspan,remove_units(p))
125+
```

src/systems/validation.jl

Lines changed: 64 additions & 36 deletions
Original file line numberDiff line numberDiff line change
@@ -1,50 +1,77 @@
11
Base.:*(x::Union{Num,Symbolic},y::Unitful.AbstractQuantity) = x * y
2+
Base.:/(x::Union{Num,Symbolic},y::Unitful.AbstractQuantity) = x / y
23

3-
"Find the units of a symbolic item."
4-
get_units(x::Real) = 1
5-
get_units(x::Unitful.Quantity) = 1 * Unitful.unit(x)
6-
get_units(x::Num) = get_units(value(x))
7-
function get_units(x::Symbolic)
4+
struct ValidationError <: Exception
5+
message::String
6+
end
7+
8+
function screen_unit(result)
9+
result isa Unitful.Unitlike || throw(ValidationError("Unit must be a subtype of Unitful.Unitlike, not $(typeof(result))."))
10+
result isa Unitful.ScalarUnits || throw(ValidationError("Non-scalar units such as $result are not supported. Use a scalar unit instead."))
11+
result == u"°" && throw(ValidationError("Degrees are not supported. Use radians instead."))
12+
end
13+
"Find the unit of a symbolic item."
14+
get_unit(x::Real) = unitless
15+
function get_unit(x::Unitful.Quantity)
16+
result = Unitful.unit(x)
17+
screen_unit(result)
18+
return result
19+
end
20+
equivalent(x,y) = isequal(1*x,1*y)
21+
unitless = Unitful.unit(1)
22+
23+
get_unit(x::Num) = get_unit(value(x))
24+
function get_unit(x::Symbolic)
825
if x isa Sym || operation(x) isa Sym || (operation(x) isa Term && operation(x).f == getindex) || x isa Symbolics.ArrayOp
926
if x.metadata !== nothing
10-
symunits = get(x.metadata, VariableUnit, 1)
27+
symunits = get(x.metadata, VariableUnit, unitless)
28+
screen_unit(symunits)
1129
else
12-
symunits = 1
30+
symunits = unitless
1331
end
14-
return oneunit(1 * symunits)
15-
elseif operation(x) isa Differential || operation(x) isa Difference
16-
return get_units(arguments(x)[1]) / get_units(arguments(arguments(x)[1])[1])
32+
return symunits
33+
elseif operation(x) isa Differential
34+
return get_unit(arguments(x)[1]) / get_unit(operation(x).x)
35+
elseif operation(x) isa Difference
36+
return get_unit(arguments(x)[1]) / get_unit(operation(x).t) #TODO: make this same as Differential
1737
elseif x isa Pow
1838
pargs = arguments(x)
19-
base,expon = get_units.(pargs)
20-
uconvert(NoUnits, expon) # This acts as an assertion
21-
return base == 1 ? 1 : operation(x)(base, pargs[2])
39+
base,expon = get_unit.(pargs)
40+
@assert expon isa Unitful.DimensionlessUnits
41+
if base == unitless
42+
unitless
43+
else
44+
pargs[2] isa Number ? operation(x)(base, pargs[2]) : operation(x)(1*base, pargs[2])
45+
end
2246
elseif x isa Add # Cannot simply add the units b/c they may differ in magnitude (eg, kg vs g)
23-
terms = get_units.(arguments(x))
24-
firstunit = 1*unit(terms[1])
47+
terms = get_unit.(arguments(x))
48+
firstunit = terms[1]
2549
for other in terms[2:end]
26-
isequal(1 * unit(other), firstunit) || throw(Unitful.DimensionError(x, terms))
50+
termlist = join(map(repr,terms),", ")
51+
equivalent(other,firstunit) || throw(ValidationError(", in sum $x, units [$termlist] do not match."))
2752
end
28-
return 1 * firstunit
53+
return firstunit
2954
elseif operation(x) == Symbolics._mapreduce
3055
if x.arguments[2] == +
31-
get_units(x.arguments[3])
56+
get_unit(x.arguments[3])
3257
else
33-
throw(ArgumentError("Unknown array operation $x"))
58+
throw(ValidationError("Unsupported array operation $x"))
3459
end
3560
else
36-
return oneunit(operation(x)(get_units.(arguments(x))...))
61+
return get_unit(operation(x)(1 .* get_unit.(arguments(x))...))
3762
end
3863
end
3964

40-
"Get units of term, returning nothing & showing warning instead of throwing errors."
41-
function safe_get_units(term, info)
65+
"Get unit of term, returning nothing & showing warning instead of throwing errors."
66+
function safe_get_unit(term, info)
4267
side = nothing
4368
try
44-
side = get_units(term)
69+
side = get_unit(term)
4570
catch err
4671
if err isa Unitful.DimensionError
4772
@warn("$info: $(err.x) and $(err.y) are not dimensionally compatible.")
73+
elseif err isa ValidationError
74+
@warn(info*err.message)
4875
elseif err isa MethodError
4976
@warn("$info: no method matching $(err.f) for arguments $(typeof.(err.args)).")
5077
else
@@ -55,7 +82,7 @@ function safe_get_units(term, info)
5582
end
5683

5784
function _validate(terms::Vector, labels::Vector{String}; info::String = "")
58-
equnits = safe_get_units.(terms, info*", ".*labels)
85+
equnits = safe_get_unit.(terms, info*" ".*labels)
5986
allthere = all(map(x -> x!==nothing, equnits))
6087
allmatching = true
6188
first_unit = nothing
@@ -64,9 +91,9 @@ function _validate(terms::Vector, labels::Vector{String}; info::String = "")
6491
if !isequal(terms[idx],0)
6592
if first_unit === nothing
6693
first_unit = equnits[idx]
67-
elseif !isequal(first_unit, equnits[idx])
94+
elseif !equivalent(first_unit, equnits[idx])
6895
allmatching = false
69-
@warn("$info: units $(equnits[1]) for $(labels[1]) and $(equnits[idx]) for $(labels[idx]) do not match.")
96+
@warn("$info: units [$(equnits[1])] for $(labels[1]) and [$(equnits[idx])] for $(labels[idx]) do not match.")
7097

7198
end
7299
end
@@ -76,8 +103,9 @@ function _validate(terms::Vector, labels::Vector{String}; info::String = "")
76103
end
77104

78105
function validate(jump::Union{ModelingToolkit.VariableRateJump, ModelingToolkit.ConstantRateJump}, t::Symbolic; info::String = "")
79-
_validate([jump.rate, 1/t], ["rate", "1/t"], info = info) && # Assuming the rate is per time units
80-
validate(jump.affect!,info = info)
106+
newinfo = replace(info,"eq."=>"jump")
107+
_validate([jump.rate, 1/t], ["rate", "1/t"], info = newinfo) && # Assuming the rate is per time units
108+
validate(jump.affect!,info = newinfo)
81109
end
82110

83111
function validate(jump::ModelingToolkit.MassActionJump, t::Symbolic; info::String = "")
@@ -91,20 +119,20 @@ function validate(jump::ModelingToolkit.MassActionJump, t::Symbolic; info::Strin
91119
end
92120

93121
function validate(jumps::ArrayPartition{<:Union{Any, Vector{<:JumpType}}}, t::Symbolic)
94-
labels = ["in Mass Action Jumps, ", "in Constant Rate Jumps, ", "in Variable Rate Jumps, "]
122+
labels = ["in Mass Action Jumps,", "in Constant Rate Jumps,", "in Variable Rate Jumps,"]
95123
all([validate(jumps.x[idx], t, info = labels[idx]) for idx in 1:3])
96124
end
97125

98-
validate(eq::ModelingToolkit.Reaction; info::String = "") = _validate([oderatelaw(eq)], ["",], info = info)
126+
validate(eq::ModelingToolkit.Reaction; info::String = "") = _validate([oderatelaw(eq)], ["reaction ",], info = info)
99127
validate(eq::ModelingToolkit.Equation; info::String = "") = _validate([eq.lhs, eq.rhs], ["left", "right"], info = info)
100128
validate(eq::ModelingToolkit.Equation, term::Union{Symbolic,Unitful.Quantity,Num}; info::String = "") = _validate([eq.lhs, eq.rhs, term], ["left","right","noise"], info = info)
101-
validate(eq::ModelingToolkit.Equation, terms::Vector; info::String = "") = _validate(vcat([eq.lhs, eq.rhs], terms), vcat(["left", "right"], "noise #".*string.(1:length(terms))), info = info)
129+
validate(eq::ModelingToolkit.Equation, terms::Vector; info::String = "") = _validate(vcat([eq.lhs, eq.rhs], terms), vcat(["left", "right"], "noise #".*string.(1:length(terms))), info = info)
102130

103131
"Returns true iff units of equations are valid."
104-
validate(eqs::Vector; info::String = "") = all([validate(eqs[idx], info = info*"in eq. #$idx") for idx in 1:length(eqs)])
105-
validate(eqs::Vector, noise::Vector; info::String = "") = all([validate(eqs[idx], noise[idx], info = info*"in eq. #$idx") for idx in 1:length(eqs)])
106-
validate(eqs::Vector, noise::Matrix; info::String = "") = all([validate(eqs[idx], noise[idx, :], info = info*"in eq. #$idx") for idx in 1:length(eqs)])
107-
validate(eqs::Vector, term::Symbolic; info::String = "") = all([validate(eqs[idx], term, info = info*"in eq. #$idx") for idx in 1:length(eqs)])
132+
validate(eqs::Vector; info::String = "") = all([validate(eqs[idx], info = info*" in eq. #$idx") for idx in 1:length(eqs)])
133+
validate(eqs::Vector, noise::Vector; info::String = "") = all([validate(eqs[idx], noise[idx], info = info*" in eq. #$idx") for idx in 1:length(eqs)])
134+
validate(eqs::Vector, noise::Matrix; info::String = "") = all([validate(eqs[idx], noise[idx, :], info = info*" in eq. #$idx") for idx in 1:length(eqs)])
135+
validate(eqs::Vector, term::Symbolic; info::String = "") = all([validate(eqs[idx], term, info = info*" in eq. #$idx") for idx in 1:length(eqs)])
108136

109137
"Throws error if units of equations are invalid."
110-
check_units(eqs...) = validate(eqs...) || throw(ArgumentError("Some equations had invalid units. See warnings for details."))
138+
check_units(eqs...) = validate(eqs...) || throw(ValidationError("Some equations had invalid units. See warnings for details."))

test/units.jl

Lines changed: 35 additions & 32 deletions
Original file line numberDiff line numberDiff line change
@@ -5,37 +5,40 @@ MT = ModelingToolkit
55
@variables t [unit = u"ms"] E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
66
D = Differential(t)
77

8-
@test MT.get_units(t) == 1u"ms"
9-
@test MT.get_units(E) == 1u"kJ"
10-
@test MT.get_units(τ) == 1u"ms"
11-
12-
@test MT.get_units(0.5) == 1.0
13-
@test MT.get_units(t) == 1.0u"ms"
14-
@test MT.get_units(P) == 1.0u"MW"
15-
@test MT.get_units(τ) == 1.0u"ms"
16-
17-
@test MT.get_units^-1) == 1/u"ms"
18-
@test MT.get_units(D(E)) == 1.0u"MW"
19-
@test MT.get_units(E/τ) == 1.0u"MW"
20-
@test MT.get_units(2*P) == 1.0u"MW"
21-
@test MT.get_units(t/τ) == 1.0
22-
@test MT.get_units(P - E/τ) == 1.0u"MW"
23-
24-
@test MT.get_units(1.0^(t/τ)) == 1.0
25-
@test MT.get_units(exp(t/τ)) == 1.0
26-
@test MT.get_units(sin(t/τ)) == 1.0
27-
@test MT.get_units(sin(1u"rad")) == 1.0
28-
@test MT.get_units(t^2) == 1.0u"ms"^2
8+
@test MT.get_unit(t) == u"ms"
9+
@test MT.get_unit(E) == u"kJ"
10+
@test MT.get_unit(τ) == u"ms"
11+
12+
@parameters β [unit = u"°"] α [unit = u"°C"] γ [unit = 1u"s"]
13+
@test_throws MT.ValidationError MT.get_unit(β)
14+
@test_throws MT.ValidationError MT.get_unit(α)
15+
@test_throws MT.ValidationError MT.get_unit(γ)
16+
17+
unitless = Unitful.unit(1)
18+
@test MT.get_unit(0.5) == unitless
19+
@test MT.get_unit(t) == u"ms"
20+
@test MT.get_unit(P) == u"MW"
21+
@test MT.get_unit(τ) == u"ms"
22+
23+
@test MT.get_unit^-1) == u"ms^-1"
24+
@test MT.equivalent(MT.get_unit(D(E)),u"MW")
25+
@test MT.equivalent(MT.get_unit(E/τ), u"MW")
26+
@test MT.get_unit(2*P) == u"MW"
27+
@test MT.get_unit(t/τ) == unitless
28+
@test MT.equivalent(MT.get_unit(P - E/τ),u"MW")
29+
@test MT.equivalent(MT.get_unit(D(D(E))),u"MW/ms")
30+
31+
@test MT.get_unit(1.0^(t/τ)) == unitless
32+
@test MT.get_unit(exp(t/τ)) == unitless
33+
@test MT.get_unit(sin(t/τ)) == unitless
34+
@test MT.get_unit(sin(1u"rad")) == unitless
35+
@test MT.get_unit(t^2) == u"ms^2"
2936

3037
@test !MT.validate(E^1.5 ~ E^(t/τ))
3138
@test MT.validate(E^(t/τ) ~ E^(t/τ))
3239

3340
eqs = [D(E) ~ P - E/τ
34-
0.0u"MW" ~ P]
35-
@test MT.get_units(eqs[1].lhs) == 1.0u"MW"
36-
@test MT.get_units(eqs[1].rhs) == 1.0u"MW"
37-
@test MT.validate(eqs[1])
38-
@test MT.validate(eqs[2])
41+
0 ~ P]
3942
@test MT.validate(eqs)
4043
@named sys = ODESystem(eqs)
4144
@named sys = ODESystem(eqs, t, [P, E], [τ])
@@ -164,23 +167,23 @@ eqs = [L ~ v*t,
164167
sys_simple = structural_simplify(sys)
165168

166169
#Jump System
167-
@parameters β [unit = 1/(u"mol"^2*u"s")] γ [unit = 1/(u"mol"*u"s")] t [unit = u"s"]
170+
@parameters β [unit = u"(mol^2*s)^-1"] γ [unit = u"(mol*s)^-1"] t [unit = u"s"] jumpmol [unit = u"mol"]
168171
@variables S(t) [unit = u"mol"] I(t) [unit = u"mol"] R(t) [unit = u"mol"]
169172
rate₁ = β*S*I
170-
affect₁ = [S ~ S - 1u"mol", I ~ I + 1u"mol"]
173+
affect₁ = [S ~ S - 1*jumpmol, I ~ I + 1*jumpmol]
171174
rate₂ = γ*I
172-
affect₂ = [I ~ I - 1u"mol", R ~ R + 1u"mol"]
175+
affect₂ = [I ~ I - 1*jumpmol, R ~ R + 1*jumpmol]
173176
j₁ = ConstantRateJump(rate₁, affect₁)
174177
j₂ = VariableRateJump(rate₂, affect₂)
175178
js = JumpSystem([j₁, j₂], t, [S, I, R], [β, γ],name=:sys)
176179

177-
affect_wrong = [S ~ S - 1u"mol", I ~ I + 1]
180+
affect_wrong = [S ~ S - jumpmol, I ~ I + 1]
178181
j_wrong = ConstantRateJump(rate₁, affect_wrong)
179-
@test_throws ArgumentError JumpSystem([j_wrong, j₂], t, [S, I, R], [β, γ],name=:sys)
182+
@test_throws MT.ValidationError JumpSystem([j_wrong, j₂], t, [S, I, R], [β, γ],name=:sys)
180183

181184
rate_wrong = γ^2*I
182185
j_wrong = ConstantRateJump(rate_wrong, affect₂)
183-
@test_throws ArgumentError JumpSystem([j₁, j_wrong], t, [S, I, R], [β, γ],name=:sys)
186+
@test_throws MT.ValidationError JumpSystem([j₁, j_wrong], t, [S, I, R], [β, γ],name=:sys)
184187

185188
# mass action jump tests for SIR model
186189
maj1 = MassActionJump(2*β/2, [S => 1, I => 1], [S => -1, I => 1])

0 commit comments

Comments
 (0)