Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Unit inconsistency leads to more confusing errors #2342

Open
mtiller-jh opened this issue Nov 1, 2023 · 13 comments
Open

Unit inconsistency leads to more confusing errors #2342

mtiller-jh opened this issue Nov 1, 2023 · 13 comments
Assignees

Comments

@mtiller-jh
Copy link

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
Copy link
Member

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

@AayushSabharwal
Copy link
Member

AayushSabharwal commented Feb 23, 2024

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
Copy link
Member

I meant the cleansyms one. That should be fixed?

@AayushSabharwal
Copy link
Member

Yep that's fixed

@bradcarman
Copy link
Contributor

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
Copy link
Member

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.

@bradcarman
Copy link
Contributor

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
Copy link
Contributor

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
Copy link
Member

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
Copy link
Member

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
Copy link
Contributor

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
Copy link

wang890 commented Jun 27, 2024

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
Copy link
Member

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.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

5 participants