I am trying to build a simulation driven by tabulated data and use unit validation. This requires interpolation of the table, and it seems DataInterpolation.jl is the intended tool (Interpolations.jl could work, but I think it needs some more dispatches defined). Below is a (not so minimal) example that illustrates the problems.
In the Source component I tried interpolating the Unitful values, which passes validation but fails during initialization in the generated function attempting to add a unitless value to the output of the interpolation which has units of voltage.
ERROR: DimensionError: -0.10400756471751627 V and 0.0 are not dimensionally compatible.
If instead, I strip units before interpolating, the system fails unit validation with
Warning: in eq. #1: units [V] for left and [] for right do not match.
In my actual code I get a third problem in the event handling code, assigning the result of the interpolation to a value on the event fails since a Unitful value cannot be stored in a unitless vector. I have not yet been able to make a minimal example of this one.
What is the right way around this? Is there some dispatch I can define to strip the units before evaluation?
using ModelingToolkit
using DataInterpolations
using Unitful
using Random
using DifferentialEquations
@variables t, [unit=u"s"]
function generate_data()
rng = MersenneTwister(54321)
time = 0.0u"s":10.0u"minute":24.0u"hr"
N = length(time)
A = 280.0u"K" .+ rand(rng, N)*u"K"
B = rand(rng, N)*u"V"
data = (; time, A, B)
end
@connector function Pin(;name)
sts = @variables begin
V(t), [unit=u"V"]
I(t), [unit=u"A", connect=Flow]
end
ODESystem(Equation[], t, sts, []; name=name)
end
function Source(data; name)
Aitp = LinearInterpolation(data.A, data.time)
# Unit mismatch on verification
# Bitp = LinearInterpolation(ustrip.(data.B), data.time)
# Unit mismatch on evaluation
Bitp = LinearInterpolation(data.B, data.time)
@parameters T=300.0, [unit=u"K"]
T = GlobalScope(T)
@named out = Pin()
@named gnd = Pin()
evts = [data.time => [T ~ Aitp(t)]]
eqs = [
out.V - gnd.V ~ Bitp(t)
0 ~ out.I + gnd.I
]
compose(ODESystem(eqs, t, [], [T]; name=name), out, gnd)
end
function Ground(;name)
@named a=Pin()
eqs = [
a.V ~ 0
]
compose(ODESystem(eqs, t, [], []; name=name), a)
end
function Resistor(R;name)
@named a=Pin()
@named b=Pin()
@parameters begin
R₀=R, [unit=u"Ω"]
γ=0.1, [unit=u"Ω/K"]
T₀=300.0, [unit=u"K"]
T, [unit=u"K"]
end
T = GlobalScope(T)
R = R₀ + (T-T₀)*γ
eqs = [
a.V - b.V ~ R*a.I
0 ~ a.I + b.I
]
compose(ODESystem(eqs, t, [], [T₀, T, R₀, γ]; name=name), a, b)
end
function Capacitor(C;name)
@named a=Pin()
@named b=Pin()
@parameters begin
C₀=C, [unit=u"F"]
γ=0.1, [unit=u"F/K"]
T₀=300.0, [unit=u"K"]
T, [unit=u"K"]
end
@variables V(t), [unit=u"V"]
D = Differential(t)
T = GlobalScope(T)
C = C₀ + (T - T₀)*γ
eqs = [
V ~ a.V - b.V
a.I + b.I ~ 0
a.I ~ C*D(V)
]
compose(ODESystem(eqs, t, [V], [T₀, T, C₀, γ]; name=name), a, b)
end
function example()
data = generate_data()
@named src = Source(data)
@named R = Resistor(1.0)
@named C = Capacitor(1.0)
@named gnd = Ground()
eqs = [
connect(src.out, R.a)
connect(R.b, C.a)
connect(C.b, gnd.a, src.gnd)
]
sys = compose(ODESystem(eqs, t, [], []; name=:example), src, R, C, gnd)
sys = structural_simplify(sys)
u0 = [ C.V => 0.0 ]
prob = ODAEProblem(sys, u0, extrema(data.time))
soln = solve(prob)
end