Using interpolation with unit verification in ModelingToolkit

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

I believe I have something working correctly.

GitHub - wnoise/mtkunitinterpolation: Julia ModelingToolkit, Unitful Unit Validation, and interpolation modifies your example (primarily in promoting T to a variable from a parameter), and the file “struct-wrapped.jl” appears to work. At least, it solves, and modifying the units in the interpolation construction make it complain about the units. No comment on whether the system it implements actually makes any sense after my modifications.

See README.md for the full explanation.

TODO: file a bug on Symbolics for @register_symbolic to work on callable structs.