Unit inconsistency leads to more confusing errors

[mtiller-jh]

This model has a legitimate issue which is that the voltage (which has units of V) is being assigned to a constant value. So let's start there:

using ModelingToolkit
using ExampleLibrary
using DifferentialEquations
using Plots
using DataInterpolations
using Unitful

@variables t [unit = u"s"]
D = Differential(t)

@mtkmodel TwoPin begin
    @components begin
        p = Pin()
        n = Pin()
    end
    @variables begin
        v(t), [unit = u"V"]
        i(t), [unit = u"A"]
    end
    @equations begin
        v ~ p.v - n.v
        0 ~ p.i + n.i
        i ~ p.i
    end
end

@mtkmodel TimeVaryingVoltage begin
    @structural_parameters begin
        f
    end
    @parameters begin
        u = 1.0, [unit = u"V"]
    end
    @extend (v,) = twopin = TwoPin()
    @equations begin
        v ~ 5.0
    end
end

time_points = [0, 0.5, 0.8, 100.0]
voltage_levels = [0, 0, 20, 20]

f = LinearInterpolation(voltage_levels, time_points)

@mtkmodel RLCModel begin
    @components begin
        source = TimeVaryingVoltage(f=f)
        ground = Ground()
    end
    @equations begin
        connect(source.n, ground.g)
    end
end

@mtkbuild rlc_model = RLCModel()
u0 = []
prob = ODEProblem(rlc_model, u0, (0, 2))
sol = solve(prob)
display(plot(sol, idxs=[rlc_model.source.v]))

This leads to the predictable:

┌ Warning:  in eq. #1: units [V] for left and [] for right do not match.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/validation.jl:176
ERROR: ModelingToolkit.ValidationError("Some equations had invalid units. See warnings for details.")
Stacktrace:
  [1] check_units(eqs::Vector{Equation})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/validation.jl:277
  [2] ODESystem(tag::UInt64, deqs::Vector{Equation}, iv::SymbolicUtils.BasicSymbolic{Real}, dvs::Vector{Any}, ps::Vector{SymbolicUtils.BasicSymbolic{Real}}, tspan::Nothing, var_to_name::Dict{Any, Any}, ctrls::Vector{Any}, observed::Vector{Equation}, tgrad::Base.RefValue{Vector{Num}}, jac::Base.RefValue{Any}, ctrl_jac::Base.RefValue{Any}, Wfact::Base.RefValue{Matrix{Num}}, Wfact_t::Base.RefValue{Matrix{Num}}, name::Symbol, systems::Vector{Any}, defaults::Dict{Any, Any}, torn_matching::Nothing, connector_type::Nothing, preface::Nothing, cevents::Vector{ModelingToolkit.SymbolicContinuousCallback}, devents::Vector{ModelingToolkit.SymbolicDiscreteCallback}, metadata::Nothing, gui_metadata::ModelingToolkit.GUIMetadata, tearing_state::Nothing, substitutions::Nothing, complete::Bool, discrete_subsystems::Nothing, unknown_states::Nothing, split_idxs::Nothing, parent::Nothing; checks::Bool)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/diffeqs/odesystem.jl:166
  [3] ODESystem
    @ ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/diffeqs/odesystem.jl:151 [inlined]
  [4] ODESystem(deqs::Vector{Equation}, iv::Num, dvs::Vector{Any}, ps::Vector{Num}; controls::Vector{Num}, observed::Vector{Equation}, systems::Vector{Any}, tspan::Nothing, name::Symbol, default_u0::Dict{Any, Any}, default_p::Dict{Any, Any}, defaults::Dict{Any, Any}, connector_type::Nothing, preface::Nothing, continuous_events::Nothing, discrete_events::Nothing, checks::Bool, metadata::Nothing, gui_metadata::ModelingToolkit.GUIMetadata)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/diffeqs/odesystem.jl:226
  [5] __TimeVaryingVoltage__(; name::Symbol, f::LinearInterpolation{Vector{Int64}, Vector{Float64}, true, Int64}, u::Nothing)
    @ Main ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/model_parsing.jl:82
  [6] __TimeVaryingVoltage__
    @ ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/model_parsing.jl:82 [inlined]
  [7] #_#172
    @ ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/model_parsing.jl:25 [inlined]
  [8] macro expansion
    @ ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/abstractsystem.jl:1003 [inlined]
  [9] __RLCModel__(; name::Symbol, source__f::Nothing)
    @ Main ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/model_parsing.jl:82
 [10] __RLCModel__
    @ ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/model_parsing.jl:82 [inlined]
 [11] #_#172
    @ ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/model_parsing.jl:25 [inlined]
 [12] top-level scope
    @ ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/abstractsystem.jl:1003

So far this makes sense and one fix would be to change the TimeVaryingVoltage model as follows:

@mtkmodel TimeVaryingVoltage begin
    @structural_parameters begin
        f
    end
    @parameters begin
        u = 1.0, [unit = u"V"]
    end
    @extend (v,) = twopin = TwoPin()
    @equations begin
        v ~ 5.0*u
    end
end

This would seem to give the system of equations the correct units. But this then leads to:

ERROR: MethodError: no method matching cleansyms(::Vector{Any})

Closest candidates are:
  cleansyms(::Tuple)
   @ SciMLBase ~/.julia/packages/SciMLBase/eK30d/src/symbolic_utils.jl:81
  cleansyms(::LinearIndices)
   @ SciMLBase ~/.julia/packages/SciMLBase/eK30d/src/symbolic_utils.jl:83
  cleansyms(::CartesianIndices)
   @ SciMLBase ~/.julia/packages/SciMLBase/eK30d/src/symbolic_utils.jl:84
  ...

Stacktrace:
 [1] macro expansion
   @ ~/.julia/packages/SciMLBase/eK30d/src/solutions/solution_interface.jl:262 [inlined]
 [2] apply_recipe(plotattributes::AbstractDict{Symbol, Any}, sol::SciMLBase.AbstractTimeseriesSolution)
   @ SciMLBase ~/.julia/packages/RecipesBase/BRe07/src/RecipesBase.jl:300
 [3] _process_userrecipes!(plt::Any, plotattributes::Any, args::Any)
   @ RecipesPipeline ~/.julia/packages/RecipesPipeline/BGM3l/src/user_recipe.jl:38
 [4] recipe_pipeline!(plt::Any, plotattributes::Any, args::Any)
   @ RecipesPipeline ~/.julia/packages/RecipesPipeline/BGM3l/src/RecipesPipeline.jl:72
 [5] _plot!(plt::Plots.Plot, plotattributes::Any, args::Any)
   @ Plots ~/.julia/packages/Plots/sxUvK/src/plot.jl:223
 [6] plot(args::Any; kw::Base.Pairs{Symbol, V, Tuple{Vararg{Symbol, N}}, NamedTuple{names, T}} where {V, N, names, T<:Tuple{Vararg{Any, N}}})
   @ Plots ~/.julia/packages/Plots/sxUvK/src/plot.jl:102
 [7] top-level scope
   @ ~/Source/modeling-toolkit-introduction/srcs/error_cases/tmp2.jl:59

However, what I really want is to be able to use the function, f, being passed in as a @structural_parameters. But that also leads to odd results. One might expect that changing 5.0 in the original model to f(t) should lead to the same units issue (since f is just a function with no particular units associated with it), but this leads to an entirely different error:

@mtkmodel TimeVaryingVoltage begin
    @structural_parameters begin
        f
    end
    @parameters begin
        u = 1.0, [unit = u"V"]
    end
    @extend (v,) = twopin = TwoPin()
    @equations begin
        v ~ f(t)
    end
end

┌ Warning:  in eq. #1right: 1.0 s and 0.0 are not dimensionally compatible.
└ @ ModelingToolkit ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/validation.jl:150
ERROR: ModelingToolkit.ValidationError("Some equations had invalid units. See warnings for details.")
Stacktrace:
  [1] check_units(eqs::Vector{Equation})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/validation.jl:277
  [2] ODESystem(tag::UInt64, deqs::Vector{Equation}, iv::SymbolicUtils.BasicSymbolic{Real}, dvs::Vector{Any}, ps::Vector{SymbolicUtils.BasicSymbolic{Real}}, tspan::Nothing, var_to_name::Dict{Any, Any}, ctrls::Vector{Any}, observed::Vector{Equation}, tgrad::Base.RefValue{Vector{Num}}, jac::Base.RefValue{Any}, ctrl_jac::Base.RefValue{Any}, Wfact::Base.RefValue{Matrix{Num}}, Wfact_t::Base.RefValue{Matrix{Num}}, name::Symbol, systems::Vector{Any}, defaults::Dict{Any, Any}, torn_matching::Nothing, connector_type::Nothing, preface::Nothing, cevents::Vector{ModelingToolkit.SymbolicContinuousCallback}, devents::Vector{ModelingToolkit.SymbolicDiscreteCallback}, metadata::Nothing, gui_metadata::ModelingToolkit.GUIMetadata, tearing_state::Nothing, substitutions::Nothing, complete::Bool, discrete_subsystems::Nothing, unknown_states::Nothing, split_idxs::Nothing, parent::Nothing; checks::Bool)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/diffeqs/odesystem.jl:166
  [3] ODESystem
    @ ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/diffeqs/odesystem.jl:151 [inlined]
  [4] ODESystem(deqs::Vector{Equation}, iv::Num, dvs::Vector{Any}, ps::Vector{Num}; controls::Vector{Num}, observed::Vector{Equation}, systems::Vector{Any}, tspan::Nothing, name::Symbol, default_u0::Dict{Any, Any}, default_p::Dict{Any, Any}, defaults::Dict{Any, Any}, connector_type::Nothing, preface::Nothing, continuous_events::Nothing, discrete_events::Nothing, checks::Bool, metadata::Nothing, gui_metadata::ModelingToolkit.GUIMetadata)
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/diffeqs/odesystem.jl:226
  [5] __TimeVaryingVoltage__(; name::Symbol, f::LinearInterpolation{Vector{Int64}, Vector{Float64}, true, Int64}, u::Nothing)
    @ Main ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/model_parsing.jl:82
  [6] __TimeVaryingVoltage__
    @ ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/model_parsing.jl:82 [inlined]
  [7] #_#172
    @ ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/model_parsing.jl:25 [inlined]
  [8] macro expansion
    @ ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/abstractsystem.jl:1003 [inlined]
  [9] __RLCModel__(; name::Symbol, source__f::Nothing)
    @ Main ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/model_parsing.jl:82
 [10] __RLCModel__
    @ ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/model_parsing.jl:82 [inlined]
 [11] #_#172
    @ ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/model_parsing.jl:25 [inlined]
 [12] top-level scope
    @ ~/.julia/packages/ModelingToolkit/7CI1l/src/systems/abstractsystem.jl:1003

The issue must be in TimeVaryingVoltage, but how is time coming into this?

[ChrisRackauckas]

@AayushSabharwal this seems like it was a SII thing. Can you check if this is fixed?

[AayushSabharwal]

It's not an SII thing. The problem is that sampling f (a LinearInterpolation) with a unit throws an error:

using ModelingToolkit: t
time_points = [0, 0.5, 0.8, 100.0]
voltage_levels = [0, 0, 20, 20]

f = LinearInterpolation(voltage_levels, time_points)
f(t)
# works, traces symbolically
f(1u"s")
# errors

[ChrisRackauckas]

I meant the cleansyms one. That should be fixed?

[AayushSabharwal]

Yep that's fixed

[bradcarman]

I think I'm getting essentially the same problem, MWE...

using ModelingToolkit
using ModelingToolkit: t, D
using DynamicQuantities

vars = @variables begin
    ẋ(t), [unit=u"m/s"]
end

eqs = [
    D(ẋ) ~ 1u"m/s^2"
]

ModelingToolkit.validate(eqs) #true
@mtkbuild sys = ODESystem(eqs, t, vars, []) #ERROR: DimensionError: 0 and 1.0 m s⁻² have incompatible dimensions

I would expect if ModelingToolkit.validate(eqs) passes then the system should build.

[ChrisRackauckas]

That's https://github.com/SciML/ModelingToolkit.jl/issues/2340 which would take a bit more work. Unit literals actually have a lot more quirks to get them right and will take quite lot more work.

[bradcarman]

OK, so I need to do the following for now...

using ModelingToolkit
using ModelingToolkit: t, D
using DynamicQuantities

vars = @variables begin
    ẋ(t), [unit=u"m/s"]
end

pars = @parameters begin
    force=1, [unit=u"m*s^-2"]
end

eqs = [
    D(ẋ) ~ force
]
@mtkbuild sys = ODESystem(eqs, t, vars, pars) # OK

OK, making progress. But, now how do I actually solve the problem?

julia> prob = ODEProblem(sys, [10], (0,1), []) 
┌ Warning:  in eq. #1right, in sum 10 - ẋ(t), units [1.0 , 1.0 m s⁻¹] do not match.
└ @ ModelingToolkit C:\Users\bradl\.julia\packages\ModelingToolkit\kByuD\src\systems\unit_check.jl:176
ERROR: ModelingToolkit.ValidationError("Some equations had invalid units. See warnings for details.")

And if I try with units I get...

julia> prob = ODEProblem(sys, [10u"m/s"], (0,1), []) 
ERROR: DimensionError: 0 and 10.0 m s⁻¹ have incompatible dimensions

[bradcarman]

OK, I see, the unit literals issue pops up with initialization...

julia> initsys = ModelingToolkit.generate_initializesystem(sys; u0map=[ẋ=>10u"m/s"])
       equations(initsys)
1-element Vector{Equation}:
 0 ~ 10.0 m s⁻¹ - ẋ(t)

[ChrisRackauckas]

There are many other issues with unit literals that would not be documented yet. There's a fundamental issue with it so it's easy to construct such edge cases.

[ChrisRackauckas]

OK, making progress. But, now how do I actually solve the problem?

For now it's assumed that you build the numerical problem with just numbers. The numerical constructors cannot handle re-validating and converting. The validation is currently too early in the process and it's not reused in the later stages yet.

[bradcarman]

OK, I see now I was defining the problem incorrectly, I needed to use a map like...

prob = ODEProblem(sys, [sys.ẋ => 10], (0, 10))

[wang890]

That's #2340 which would take a bit more work. Unit literals actually have a lot more quirks to get them right and will take quite lot more work.

The unit is useful. I am currently trying the following solution

using ModelingToolkit: t, D # with DynamicQuantities Unit

step 1：construct a standard unit dictionary SI like the type class of Modelica (ctr+F to search the words: package SI) SI = Dict( "Length" => Dict("quantity"=>"Length", "unit"=>"m", "u"=>"u\"m\"", "prefixes" => Dict("quantity"=>["final"], "unit"=>["final"], )),

"Area" => Dict("quantity"=>"Area", "unit"=>"m2", "u"=>"u\"m^2\"", "prefixes" => Dict("quantity"=>["final"], "unit"=>["final"], )) ) The type class of Modelica is similar with Julia Struct, I think, but the Metadata style in ModelingToolkit is also very good.

Step2: function u as following u(u_name) = SI[u_name]["u"], then u("Area") returns the u\"m^2\" which can be recognized by DynamicQuantities.jl

Step3: using as following @mtkmodel VoltageSource begin @extend i,v,p,n = onePort = OnePort() @components begin
ramp = Ramp()
end @constants begin unit_balance = 1.0, [unit = u("ElectricPotential"), description = "Balance out the units of the equation"] end
@equations begin
v ~ ramp.output.u * unit_balance # v is a ElectricPotential quantity
end end

However, some models have many equations, many of which have the unitless variables from Blocks. It is very troublesome to balance the units, So I commented out the body of function check_units, because the param checks can not be assigned during @mtkbuild sys = Model() like ODESystem(eqs, t, checks=MT.CheckComponents)

In addition, standard unit must be used at present. In future, it would be considered that customizing units of different scales, such as nm for microscopic, like Unitful.preferunits(u"nm,s,A,K,cd,mg,mol"...), also like MTK docs say "In the future, the validation stage may be upgraded to support the insertion of conversion factors into the equations"

Unit is only auxiliary to modeling and solving, and could be improved in future. I am a new user of Julia and MTK, and what I said above may not be correct. Thanks a lot for your great work. @ChrisRackauckas @YingboMa

[ChrisRackauckas]

That was all pretty correct. Yeah we need to improve a few things. We don't have the modelica wildcard units, which seem to be under-defined, but the bigger difficulty is we will need to build a new unit system or modify DynamicQuantities.jl if we want to add that. So it'll take a bit to get this right.