Closed YingboMa closed 8 months ago
Should adding the guess
syntax be considered a breaking change? Because when that's done, the default should take on a stronger meaning of always being satisfied, which would be a breaking change.
Anything that causes a simulation that previously ran and produced a correct answer to not run or not produce a correct answer should be breaking?
What's correct in this case is perhaps a bit fuzzy since if the default previously wasn't satisfied, it couldn't really be considered a correct answer, but it could still have been what the user wanted.
I guess the addition of guess
has to be paired with a change to MTKstdlib to make all default values added to library components guesses rather than fixed.
What's correct in this case is perhaps a bit fuzzy since if the default previously wasn't satisfied, it couldn't really be considered a correct answer, but it could still have been what the user wanted.
The default had a different meaning. It was either an initial condition or a guess depending on the DAE solver's chosen initialization scheme. So if the initial value ended up different from the default (which was actually quite common, it was the behavior on all algebraic equations) then it would not error and that was considered correct behavior.
In v9 we plan to make it so that if you set a default, that has to be satisfied, and only guess
es are allowed to change. That would break code that implicitly was using default to set a guess
(like old examples of Pin
).
Another breaking change, seen in units, is that we want the following code to work:
using ModelingToolkit
using DifferentialEquations
using Plots
using Unitful
@variables t [unit = u"s"]
D = Differential(t)
@connector Pin begin
v(t), [unit = u"V"]
i(t), [connect = Flow, unit = u"A"]
end
@mtkmodel Ground begin
@components begin
g = Pin()
end
@equations begin
g.v ~ 0
end
end
@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 Resistor begin
@extend v, i = twopin = TwoPin()
@parameters begin
R, [unit = u"ฮฉ"]
end
@equations begin
v ~ i * R
end
end
@mtkmodel Capacitor begin
@extend v, i = twopin = TwoPin()
@parameters begin
C, [unit = u"F"]
end
@equations begin
D(v) ~ i / C
end
end
@mtkmodel ConstantVoltage begin
@extend (v,) = twopin = TwoPin()
@parameters begin
V, [unit = u"V"]
end
@equations begin
V ~ v
end
end
@mtkmodel RCModel begin
@components begin
resistor = Resistor(R = 1.0u"ฮฉ")
capacitor = Capacitor(C = 1.0u"F")
source = ConstantVoltage(V = 1.0u"V")
ground = Ground()
end
@equations begin
connect(source.p, resistor.p)
connect(resistor.n, capacitor.p)
connect(capacitor.n, source.n)
connect(capacitor.n, ground.g)
end
end
@mtkbuild rc_model = RCModel(resistor.R = 2.0u"ฮฉ")
u0 = [
rc_model.capacitor.v => 0.0u"V",
]
prob = ODEProblem(sys, u0, (0u"s", 10.0u"s"))
sol = solve(prob)
Right now this code does not work because defaults set to be unitful will cause the simulation process itself to be unitful.
julia> prob = ODEProblem(rc_model, u0, (0u"s", 10.0u"s"))
ODEProblem with uType Vector{Quantity{Float64, ๐ยฒ ๐ ๐โปยน ๐โปยณ, Unitful.FreeUnits{(V,), ๐ยฒ ๐ ๐โปยน ๐โปยณ, nothing}}} and tType Quantity{Float64, ๐, Unitful.FreeUnits{(s,), ๐, nothing}}. In-place: true
timespan: (0.0 s, 10.0 s)
u0: 1-element Vector{Quantity{Float64, ๐ยฒ ๐ ๐โปยน ๐โปยณ, Unitful.FreeUnits{(V,), ๐ยฒ ๐ ๐โปยน ๐โปยณ, nothing}}}:
0.0 V
However, the purpose of the MTK unit handling is to decouple the unit validation process from the solving process. The reason why this fails is technically "just" a solver issue, it would actually work on non-DAE ODEs if you use an explicit RK method but fails on implicit solvers due to linear algebra not supporting unitful. But this is exactly the reason for attempting to perform canonicalization at the MTK level, since making the solvers all support unitful is would not only be painful but also an unnecessary performance hit.
Today, the syntax that will work is:
@mtkbuild rc_model = RCModel(resistor.R = 2.0)
u0 = [
rc_model.capacitor.v => 0.0,
]
prob = ODEProblem(rc_model, u0, (0, 10.0))
sol = solve(prob)
i.e. you build the model with units, but when you actually specify values you don't give them units. This is effectively the canonicalization done by hand, where the model has units (and thus supports the validation tooling) but the solving process does not validate on every function call.
However, the better thing to do would be to require that the user passes in unitful values were are then validated against the units defined in the metadata and then stripped during the XProblem building. This would make there be clearer error reporting back to the user, and allow for automatic converting if different scales are chosen, and we should then error if a unitless value is given as a default or guess for a unitful variable. That said, adding that error would be breaking, given that this is the syntax that works today.
@baggepinnen could you make a PR to remove all the old discrete-time systems code?
Can someone explain exactly what "Default parameter splitting" means?
@baggepinnen could you make a PR to remove all the old discrete-time systems code?
Much of this resides in Symbolics.jl difference.jl
, and it looks like some of it is used for FMUs?
Much of this resides in Symbolics.jl difference.jl, and it looks like some of it is used for FMUs?
@sharanry ?
Much of this resides in Symbolics.jl difference.jl, and it looks like some of it is used for FMUs?
Note that FMUs.jl is no longer in use. We no longer use Symbolics for anything FMU related AFAIK.
Alright just drop it.
states
-> unknowns
: https://github.com/SciML/ModelingToolkit.jl/pull/2432ODAEProblem
: https://github.com/SciML/ModelingToolkit.jl/pull/2435complete
d before creating XProblem
: https://github.com/SciML/ModelingToolkit.jl/pull/2436
states
->unknowns
ODESystem
Unitful
->DynamicQuantities
XProblem
s.ODAEProblem
to justODEProblem