Skip to content

Commit 88fe45e

Browse files
Merge pull request #1156 from lamorton/uniting
Validation of Unitful equations
2 parents dd64720 + 3e32fd2 commit 88fe45e

File tree

15 files changed

+519
-85
lines changed

15 files changed

+519
-85
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/ModelingToolkit.jl

Lines changed: 2 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -122,7 +122,6 @@ include("systems/diffeqs/sdesystem.jl")
122122
include("systems/diffeqs/abstractodesystem.jl")
123123
include("systems/diffeqs/first_order_transform.jl")
124124
include("systems/diffeqs/modelingtoolkitize.jl")
125-
include("systems/diffeqs/validation.jl")
126125
include("systems/diffeqs/basic_transformations.jl")
127126

128127
include("systems/jumps/jumpsystem.jl")
@@ -135,10 +134,9 @@ include("systems/control/controlsystem.jl")
135134

136135
include("systems/pde/pdesystem.jl")
137136

138-
include("systems/dependency_graphs.jl")
139-
140137
include("systems/discrete_system/discrete_system.jl")
141-
138+
include("systems/validation.jl")
139+
include("systems/dependency_graphs.jl")
142140
include("systems/systemstructure.jl")
143141
using .SystemStructures
144142

src/systems/abstractsystem.jl

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -206,9 +206,10 @@ Setfield.get(obj::AbstractSystem, ::Setfield.PropertyLens{field}) where {field}
206206
:(getfield(obj, $(Meta.quot(fn))))
207207
end
208208
end
209+
kwarg = :($(Expr(:kw, :checks, false))) # Inputs should already be checked
209210
return Expr(:block,
210211
Expr(:meta, :inline),
211-
Expr(:call,:(constructorof($obj)), args...)
212+
Expr(:call, :(constructorof($obj)), args..., kwarg)
212213
)
213214
else
214215
error("This should never happen. Trying to set $(typeof(obj)) with $patch.")

src/systems/control/controlsystem.jl

Lines changed: 12 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -72,11 +72,14 @@ struct ControlSystem <: AbstractControlSystem
7272
parameters are not supplied in `ODEProblem`.
7373
"""
7474
defaults::Dict
75-
function ControlSystem(loss, deqs, iv, dvs, controls, ps, observed, name, systems, defaults)
76-
check_variables(dvs, iv)
77-
check_parameters(ps, iv)
78-
check_equations(deqs, iv)
79-
check_equations(observed, iv)
75+
function ControlSystem(loss, deqs, iv, dvs, controls, ps, observed, name, systems, defaults; checks::Bool = true)
76+
if checks
77+
check_variables(dvs, iv)
78+
check_parameters(ps, iv)
79+
check_equations(deqs, iv)
80+
check_equations(observed, iv)
81+
check_units(deqs)
82+
end
8083
new(loss, deqs, iv, dvs, controls, ps, observed, name, systems, defaults)
8184
end
8285
end
@@ -87,7 +90,8 @@ function ControlSystem(loss, deqs::AbstractVector{<:Equation}, iv, dvs, controls
8790
default_u0=Dict(),
8891
default_p=Dict(),
8992
defaults=_merge(Dict(default_u0), Dict(default_p)),
90-
name=nothing)
93+
name=nothing,
94+
kwargs...)
9195
name === nothing && throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro"))
9296
if !(isempty(default_u0) && isempty(default_p))
9397
Base.depwarn("`default_u0` and `default_p` are deprecated. Use `defaults` instead.", :ControlSystem, force=true)
@@ -105,7 +109,7 @@ function ControlSystem(loss, deqs::AbstractVector{<:Equation}, iv, dvs, controls
105109
collect_defaults!(defaults, dvs′)
106110
collect_defaults!(defaults, ps′)
107111
ControlSystem(value(loss), deqs, iv′, dvs′, controls′,
108-
ps′, observed, name, systems, defaults)
112+
ps′, observed, name, systems, defaults, kwargs...)
109113
end
110114

111115
struct ControlToExpr
@@ -197,5 +201,5 @@ function runge_kutta_discretize(sys::ControlSystem,dt,tspan;
197201
equalities = vcat(stages,updates,control_equality)
198202
opt_states = vcat(reduce(vcat,reduce(vcat,states_timeseries)),reduce(vcat,reduce(vcat,k_timeseries)),reduce(vcat,reduce(vcat,control_timeseries)))
199203

200-
OptimizationSystem(reduce(+,losses, init=0),opt_states,ps,equality_constraints = equalities, name=nameof(sys))
204+
OptimizationSystem(reduce(+,losses, init=0),opt_states,ps,equality_constraints = equalities, name=nameof(sys), checks = false)
201205
end

src/systems/diffeqs/odesystem.jl

Lines changed: 9 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -88,10 +88,13 @@ struct ODESystem <: AbstractODESystem
8888
"""
8989
preface::Any
9090

91-
function ODESystem(deqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, structure, connection_type, preface)
92-
check_variables(dvs,iv)
93-
check_parameters(ps,iv)
94-
check_equations(deqs,iv)
91+
function ODESystem(deqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, structure, connection_type, preface; checks::Bool = true)
92+
if checks
93+
check_variables(dvs,iv)
94+
check_parameters(ps,iv)
95+
check_equations(deqs,iv)
96+
check_units(deqs)
97+
end
9598
new(deqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, structure, connection_type, preface)
9699
end
97100
end
@@ -107,6 +110,7 @@ function ODESystem(
107110
defaults=_merge(Dict(default_u0), Dict(default_p)),
108111
connection_type=nothing,
109112
preface=nothing,
113+
checks = true,
110114
)
111115
name === nothing && throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro"))
112116
deqs = collect(deqs)
@@ -140,7 +144,7 @@ function ODESystem(
140144
if length(unique(sysnames)) != length(sysnames)
141145
throw(ArgumentError("System names must be unique."))
142146
end
143-
ODESystem(deqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, nothing, connection_type, preface)
147+
ODESystem(deqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, nothing, connection_type, preface, checks = checks)
144148
end
145149

146150
function ODESystem(eqs, iv=nothing; kwargs...)

src/systems/diffeqs/sdesystem.jl

Lines changed: 9 additions & 5 deletions
Original file line numberDiff line numberDiff line change
@@ -86,10 +86,13 @@ struct SDESystem <: AbstractODESystem
8686
"""
8787
connection_type::Any
8888

89-
function SDESystem(deqs, neqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, connection_type)
90-
check_variables(dvs,iv)
91-
check_parameters(ps,iv)
92-
check_equations(deqs,iv)
89+
function SDESystem(deqs, neqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, connection_type; checks::Bool = true)
90+
if checks
91+
check_variables(dvs,iv)
92+
check_parameters(ps,iv)
93+
check_equations(deqs,iv)
94+
check_units(deqs,neqs)
95+
end
9396
new(deqs, neqs, iv, dvs, ps, var_to_name, ctrls, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, connection_type)
9497
end
9598
end
@@ -103,6 +106,7 @@ function SDESystem(deqs::AbstractVector{<:Equation}, neqs, iv, dvs, ps;
103106
defaults=_merge(Dict(default_u0), Dict(default_p)),
104107
name=nothing,
105108
connection_type=nothing,
109+
checks = true,
106110
)
107111
name === nothing && throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro"))
108112
deqs = collect(deqs)
@@ -130,7 +134,7 @@ function SDESystem(deqs::AbstractVector{<:Equation}, neqs, iv, dvs, ps;
130134
ctrl_jac = RefValue{Any}(Matrix{Num}(undef, 0, 0))
131135
Wfact = RefValue(Matrix{Num}(undef, 0, 0))
132136
Wfact_t = RefValue(Matrix{Num}(undef, 0, 0))
133-
SDESystem(deqs, neqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, connection_type)
137+
SDESystem(deqs, neqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, tgrad, jac, ctrl_jac, Wfact, Wfact_t, name, systems, defaults, connection_type, checks = checks)
134138
end
135139

136140
SDESystem(sys::ODESystem, neqs; kwargs...) = SDESystem(equations(sys), neqs, get_iv(sys), states(sys), parameters(sys); kwargs...)

src/systems/diffeqs/validation.jl

Lines changed: 0 additions & 26 deletions
This file was deleted.

src/systems/discrete_system/discrete_system.jl

Lines changed: 8 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -54,9 +54,12 @@ struct DiscreteSystem <: AbstractTimeDependentSystem
5454
in `DiscreteSystem`.
5555
"""
5656
default_p::Dict
57-
function DiscreteSystem(discreteEqs, iv, dvs, ps, var_to_name, ctrls, observed, name, systems, default_u0, default_p)
58-
check_variables(dvs,iv)
59-
check_parameters(ps,iv)
57+
function DiscreteSystem(discreteEqs, iv, dvs, ps, var_to_name, ctrls, observed, name, systems, default_u0, default_p; checks::Bool = true)
58+
if checks
59+
check_variables(dvs, iv)
60+
check_parameters(ps, iv)
61+
check_units(discreteEqs)
62+
end
6063
new(discreteEqs, iv, dvs, ps, var_to_name, ctrls, observed, name, systems, default_u0, default_p)
6164
end
6265
end
@@ -75,6 +78,7 @@ function DiscreteSystem(
7578
default_u0=Dict(),
7679
default_p=Dict(),
7780
defaults=_merge(Dict(default_u0), Dict(default_p)),
81+
kwargs...,
7882
)
7983
name === nothing && throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro"))
8084
eqs = collect(eqs)
@@ -97,7 +101,7 @@ function DiscreteSystem(
97101
if length(unique(sysnames)) != length(sysnames)
98102
throw(ArgumentError("System names must be unique."))
99103
end
100-
DiscreteSystem(eqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, name, systems, default_u0, default_p)
104+
DiscreteSystem(eqs, iv′, dvs′, ps′, var_to_name, ctrl′, observed, name, systems, default_u0, default_p, kwargs...)
101105
end
102106

103107
"""

src/systems/jumps/jumpsystem.jl

Lines changed: 8 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -53,9 +53,12 @@ struct JumpSystem{U <: ArrayPartition} <: AbstractTimeDependentSystem
5353
type: type of the system
5454
"""
5555
connection_type::Any
56-
function JumpSystem{U}(ap::U, iv, states, ps, var_to_name, observed, name, systems, defaults, connection_type) where U <: ArrayPartition
57-
check_variables(states, iv)
58-
check_parameters(ps, iv)
56+
function JumpSystem{U}(ap::U, iv, states, ps, var_to_name, observed, name, systems, defaults, connection_type; checks::Bool = true) where U <: ArrayPartition
57+
if checks
58+
check_variables(states, iv)
59+
check_parameters(ps, iv)
60+
check_units(ap,iv)
61+
end
5962
new{U}(ap, iv, states, ps, var_to_name, observed, name, systems, defaults, connection_type)
6063
end
6164
end
@@ -68,6 +71,7 @@ function JumpSystem(eqs, iv, states, ps;
6871
defaults=_merge(Dict(default_u0), Dict(default_p)),
6972
name=nothing,
7073
connection_type=nothing,
74+
checks = true,
7175
kwargs...)
7276
name === nothing && throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro"))
7377
eqs = collect(eqs)
@@ -98,7 +102,7 @@ function JumpSystem(eqs, iv, states, ps;
98102
process_variables!(var_to_name, defaults, states)
99103
process_variables!(var_to_name, defaults, ps)
100104

101-
JumpSystem{typeof(ap)}(ap, value(iv), states, ps, var_to_name, observed, name, systems, defaults, connection_type)
105+
JumpSystem{typeof(ap)}(ap, value(iv), states, ps, var_to_name, observed, name, systems, defaults, connection_type, checks = checks)
102106
end
103107

104108
function generate_rate_function(js::JumpSystem, rate)

0 commit comments

Comments
 (0)