Interpolation within ModelingToolkit framework

Hello,

Is there a way to interpolate within a ModelingToolkit’s (MTK) component model using the parameters/variables of the model? Interpolations.jl seems to struggle with Num types, and from the example in the documentation (Composing Ordinary Differential Equations · ModelingToolkit.jl) I understand that it is only possible to use pre-interpolated data that is then fed to the model as a time-variable function (?).

To illustrate this idea, I include a simple example using t as coordinate to obtain the interpolated variable.

using ModelingToolkit
using DifferentialEquations
using Interpolations


# Function to interpolate a time series
function interpolateTimeSeries(t, time, values)

	interpolator = CubicSplineInterpolation(time, values)
	interpolatedValue = interpolator(t)

end

@register interpolateTimeSeries(t, time, values)

# Source to produce the output of a timeseries table
function lookupTable(;name, time, values)
	
	@parameters t
	@variables output(t)

	equations = [output ~ interpolateTimeSeries(t, time, values)]

	ODESystem(equations, t, output, []; name=name)

end

############################## lookupTable source model #################
# Data
time = LinRange(1, 100, 100)
values = randn(100)

# Instantiate lookup-table source model
@named sourceLookupTable = lookupTable(time = time, values = values)
simpSource = structural_simplify(sourceLookupTable)
u0 = [sourceLookupTable => ]

# Define the problem
prob = ODEProblem(simpSource, u0, (0.0, 100.0))
sol = solve(prob, Tsit5())

This produces the following error:

ERROR: LoadError: TypeError: non-boolean (Num) used in boolean context
Stacktrace:
 [1] inbounds_index at $PATH\.julia\packages\Interpolations\qHlUr\src\extrapolation\extrapolation.jl:108 [inlined]
 [2] inbounds_position at $PATH\.julia\packages\Interpolations\qHlUr\src\extrapolation\extrapolation.jl:99 [inlined]
 [3] Extrapolation at $PATH\.julia\packages\Interpolations\qHlUr\src\extrapolation\extrapolation.jl:48 [inlined]
 [4] interpolateTimeSeries(::Num, ::LinRange{Float64}, ::Array{Float64,1}) at $PATH\SourcesSinks.jl:10
 [5] lookupTable(; name::Symbol, time::LinRange{Float64}, values::Array{Float64,1}) at $PATH\SourcesSinks.jl:22
 [6] top-level scope at $PATH\SourcesSinks.jl:34
 [7] include(::String) at .\client.jl:457
 [8] top-level scope at REPL[7]:1
in expression starting at $PATH\SourcesSinks.jl:34

That works IIRC? But even if it doesn’t, you can just @register your function.

Thanks for the quick reply @ChrisRackauckas. I did @register my function in the previous example (maybe not correctly?).

I also tried DataInterpolations.jl as you suggested, both with and without an external registered function, but I get the same error as in the example included in the first post.

Without external function

using ModelingToolkit
using DifferentialEquations
using Interpolations
using DataInterpolations


function lookupTable(;name, time, values)
	
	@parameters t
	@variables output(t)
 	interpolator = CubicSpline(values, time)

	equations = [output - interpolator(t)]

	ODESystem(equations, t, [output], []; name=name)

end

# Data
time = LinRange(1, 100, 100)
values = randn(100)

# Instantiate lookup-table source model
@named sourceLookupTable = lookupTable(time = time, values = values)
simpSource = structural_simplify(sourceLookupTable)

With external registered function

# Function to interpolate a time series
function interpolateTimeSeries(t, time, values)

	interpolator = CubicSpline(values, time)
	interpolatedValue = interpolator(t)

end

@register interpolateTimeSeries(t, time, values)

# Source to produce the output of a timeseries table
function lookupTable(;name, time, values)
	
	@parameters t
	@variables output(t)

	equations = [output ~ interpolateTimeSeries(t, time, values)]

	ODESystem(equations, t, [output], []; name=name)

end

############################## lookupTable source model #################
# Data
time = LinRange(1, 100, 100)
values = randn(100)

# Instantiate lookup-table source model
@named sourceLookupTable = lookupTable(time = time, values = values)
#@named sourceLookupTable = lookupTable()
simpSource = structural_simplify(sourceLookupTable)

In both cases, and the first example using Interpolations.jl, I get the same error:

ERROR: LoadError: TypeError: non-boolean (Num) used in boolean context

These errors originate from the use of @parameters as an input during the interpolation, e.g. interpolator(t) or interpolateTimeSeries(t, time, values) .

You just forgot to specify the non-symbolic parts arguments. The fix is:

@register interpolateTimeSeries(t, time::AbstractVector, values::AbstractVector)

Many thanks for your help @ChrisRackauckas. The work of the MKT team in this package is fantastic.

Just some constructive feedback. I feel that the @register functionality might be useful to many people. However, a quick search of @register in the search docs bar only leads to Composing Ordinary Differential Equations - Specifying a time-variable forcing function, where it says:

" MTK allows to “register” arbitrary Julia functions, which are excluded from symbolic transformations but are just used as-is. So, you could, for example, interpolate a given time series using DataInterpolations.jl."

From this description, at least for a new user of MTK like me, it is a bit difficult to infer that I can specify some variables as non-symbolic and the remaining as symbolic. I took the liberty of using the example included in the documentation and include a function that interpolates some time-series data with a mix of non-symbolic and symbolic inputs. Here is the code:

using ModelingToolkit
using DifferentialEquations
using DataInterpolations 


@parameters tau
@variables t, x(t), f(t)
D = Differential(t)

function f_interpolate(t, table_t, table_u)

	interpolator = QuadraticInterpolation(table_u, table_t)
	output = interpolator(t)

end

@register f_interpolate(t, table_t::AbstractVector, table_u::AbstractVector)

table_t, table_u = LinRange(0.5, 10.5, 10), randn(10)

@named fol_external_f = ODESystem([f ~ f_interpolate(t, table_t, table_u), D(x) ~ (f - x) / tau])
prob = ODEProblem(structural_simplify(fol_external_f), [x => 0.0], (0.0, 10.0), [tau => 0.75])

sol = solve(prob)
plot(sol, vars=[x, f])
scatter!(table_t, table_u)

Capture3

If you feel this may complement the existing documentation, feel free to add it. Otherwise it is left here so it may help future new users.

That would be a good tutorial: Mixing Symbolic Models with Data, or something.

2 Likes

This is just a performance tip but I would have f_interpolate return the interpolating function. Right now every time f_interpolate is called it constructs a interpolation object. I would rewrite it as follows

f_interpolate = QuadraticInterpolation(table_u, table_t)
@register f_interpolate(t)

Assuming @register accepts DataInterpolations

It accepts any Julia function.

Agree. Thanks for the remark. The code above can become inefficient when one considers arrays.

I tried some of the above answers, but I too am running into errors trying to include CubicSpline interpolation from DataInterpolations.jl in my loss function. In the example below I am interpolating the model onto the same grid just as a test. Thank you in advance!

julia> using ModelingToolkit, DataInterpolations, Polynomials

julia> @variables c0 c1
2-element Vector{Num}:
 c0
 c1

julia> @parameters x[1:10] data[1:10]
2-element Vector{Symbolics.Arr{Num, 1}}:
 x[1:10]
 data[1:10]

julia> function f_interpolate(t, table_t, table_u)

               interpolator = CubicSpline(table_u, table_t)
               output = interpolator.(t)

       end
f_interpolate (generic function with 1 method)

julia> rms(x,y) = sqrt(sum((x .- y).^2) / length(x))
rms (generic function with 1 method)

julia> typeof(x)
Symbolics.Arr{Num, 1}

julia> @register f_interpolate(t::Vector{Symbolics.Arr{Num, 1}}, table_t::Vector{Symbolics.Arr{Num, 1}}, table_u::Vector{Symbolics.Arr{Num, 1}})

julia> loss = rms(data, f_interpolate(x, x, Polynomial([c0, c1]).(x)))
ERROR: MethodError: no constructors have been defined for Any
Stacktrace:
 [1] _broadcast_getindex_evalf
   @ ./broadcast.jl:670 [inlined]
 [2] _broadcast_getindex
   @ ./broadcast.jl:643 [inlined]
 [3] getindex
   @ ./broadcast.jl:597 [inlined]
 [4] copy
   @ ./broadcast.jl:899 [inlined]
 [5] materialize(bc::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{1}, Nothing, Type{Any}, Tuple{Vector{Term{Any, Nothing}}}})
   @ Base.Broadcast ./broadcast.jl:860
 [6] munge_data(u::Symbolics.Arr{Any, 1}, t::Symbolics.Arr{Num, 1})
   @ DataInterpolations ~/.julia/packages/DataInterpolations/HJsu4/src/interpolation_utils.jl:44
 [7] CubicSpline(u::Symbolics.Arr{Any, 1}, t::Symbolics.Arr{Num, 1})
   @ DataInterpolations ~/.julia/packages/DataInterpolations/HJsu4/src/interpolation_caches.jl:161
 [8] f_interpolate(t::Symbolics.Arr{Num, 1}, table_t::Symbolics.Arr{Num, 1}, table_u::Symbolics.Arr{Any, 1})
   @ Main ./REPL[4]:3
 [9] top-level scope
   @ REPL[9]:1

julia>