I have an ODE system defined through modelingtoolkit. The system is formed of multiple ODESystems and simplified with system = structural_simplify(all).
One of the ODESystem is defined by the following function
function create_ref(func,lookahead_times; name = :ref)
@parameters t
@variables o[1:length(lookahead_times)+1](t)
eqs = vcat(
o[1] ~ func(t),
[o[i+1] ~ func(t+lookahead_times[i]) for i in 1:length(lookahead_times)]
)
ODESystem(eqs; name=name)
end
where func is a function passed as an argument.
Initially, I had func be essentially a sum of sinusoidal (see code below).
However, I now want this function to be a random signal band-limited signal. In particular, I want it to be a filtered ornestein-ulhnebeck process. The issue is that I can’t define such a system as an abstract function as the sum of sinusoidalas above. I need to interpolate the solution of and OU process (see code below).
So essentially I just changed the func in the ODESystem above from a sum of sins to an interpolated Ornestein-Ulhenbeck process.
This change slows down the solution of my overall system immensely.
I ran a profiler and found that most of the run time is spent with the check_gridded function of the interpolate (see picture below for part of the flamegraph).
I was wondering if there is a better interpolator I can use to define the function that will slow down less the solution of the ODE. Or even, if there might be a way of defining the ODEsystem without the need of the interpolator.
Thank for the help in advance.
Definition of sum of sinusoidals
function sinFunDef(fc,numSin=4,dt=0.001)
omegas = rand(fc/(numSin+1):dt:fc,numSin)
phaseDif = rand(0:dt:pi/2,numSin)
amp = rand(0.1:dt:1,numSin)
function refSin(t) # sum of sinusoidals
sum(1/numSin*amp.*sin.(omegas.*t.+phaseDif))
end
return refSin
end
Definition of interpolated O-U process
function interpolateSeries(t,time,values)
# interpolator = CubicSplineInterpolation(time, values)
interpolator = LinearInterpolation(time, values)
interpolator(t)
end
@register interpolateSeries(t, time::AbstractVector, values::AbstractVector)
"""
create random reference from O-U process filtered and interpolatedValue
"""
function ouFunDef(θ,sigma,μ,T,dt,fc,w0S=0.5)
# ornstein uhlenbeck process
t0 = 0.0
w0 = w0S*(rand()-0.5)
W = OrnsteinUhlenbeckProcess(θ,μ,sigma,t0,w0)
prob = NoiseProblem(W,(t0,T))
sol = DifferentialEquations.solve(prob;dt=dt)
# filter
responsetype = Lowpass(fc) # type of filter
designmethod = Butterworth(3) # window type
filteredX = filt(digitalfilter(responsetype, designmethod), sol.u)
function refF(t)
interpolateSeries(t,sol.t,filteredX)
end
return refF
end