Wrap my function in an MTK component

Specific Domains Modelling & Simulations
1

I am trying to wrap an arbitrary function in an MTK component and I am a bit lost in all the possibilities available these days.

Right now I have:

seconds_to_datetime(t, t0::DateTime) = t0 + Millisecond(round(Int, t * 1e3))

function SolarPosition.SolarPositionComponent(;
    name,
    t0::DateTime,
    observer::Observer,
    algorithm::SolarAlgorithm = PSA(),
    refraction::RefractionAlgorithm = NoRefraction(),
)
    pos(t) = solar_position(observer, seconds_to_datetime(t, t0), algorithm, refraction)

    # and now define the component?
end

where solar_position will return a struct like:

struct SolPos{T} <: AbstractSolPos where {T<:AbstractFloat}
    "Azimuth (degrees, 0=N, +clockwise, range [-180, 180])"
    azimuth::T
    "Elevation (degrees, range [-90, 90])"
    elevation::T
    "Zenith = 90 - elevation (degrees, range [0, 180])"
    zenith::T
end

I’m trying to build a component, where I can just refer to positions like:

@named sun = SolarPositionComponent(t0 = t0, observer = obs)
sun.azimuth
sun.elevation

Some pointing to the right docs and examples would be appreciated.

1 Like
2

An example would be ModelingToolkitNeuralNets.jl where NeuralNetworkBlock wraps a function call.

What you need to do it to create a System with a pos ~ solar_position(observer, ...) equation (or whatever relations would define the system) and then you would define the parameters of your system as the non-numeric constants that your function needs.

Currently using the struct symbolically might be a bit tricky, but if you have just have 3 separate variables, then the last code block would just work since azimuth and elevation would just be @variables in the system and you would have some equations for them.

2 Likes
3

do you mean in an @parameters block inside the SolarPositionComponent?

Sorry I can’t really follow. Shouldn’t azimuth and elevation already be @variables inside my component?

4

I tried the following approach:


seconds_to_datetime(t_sec, t0::DateTime) = t0 + Millisecond(round(Int, t_sec * 1e3))

@register_symbolic seconds_to_datetime(t_sec::Num, t0)
@register_symbolic solar_position(observer, time, algorithm, refraction)

function SolarPosition.SolarPositionBlock(; name)

    @parameters t0::DateTime [tunable = false] observer::Observer [tunable = false]
    @parameters algorithm::SolarAlgorithm = PSA() [tunable = false]
    @parameters refraction::RefractionAlgorithm = NoRefraction() [tunable = false]

    @variables pos(t) [output = true]

    eqs =
        [pos ~ solar_position(observer, seconds_to_datetime(t, t0), algorithm, refraction)]

    return System(eqs, t; name = name)
end

Which then fails with

  Got exception outside of a @test
  MethodError: no method matching seconds_to_datetime(::Symbolics.Num, ::SymbolicUtils.BasicSymbolic{DateTime})
  The function `seconds_to_datetime` exists, but no method is defined for this combination of argument types.

I’m really not sure what I’m supposed to do here.

5

For the function registration, try to register with

@register_symbolic seconds_to_datetime(t_sec, t0::DateTime)
@register_symbolic solar_position(observer::Observer, time::DateTime, algorithm::SolarAlgorithm, refraction::RefractionAlgorithm)

i.e. specify all the types that are not symbolic

1 Like
6

thank you!

I came to the following:

using SolarPosition: Observer, SolarAlgorithm, RefractionAlgorithm, PSA, NoRefraction
using ModelingToolkit: @parameters, @variables, System, @register_symbolic, t_nounits
using Symbolics
using Symbolics: term
using Dates: DateTime, Millisecond

import SolarPosition: SolarPositionBlock, solar_position


seconds_to_datetime(t_sec, t0::DateTime) = t0 + Millisecond(round(Int, t_sec * 1e3))

@register_symbolic seconds_to_datetime(t_sec, t0::DateTime)
@register_symbolic solar_position(
    observer::Observer,
    time::DateTime,
    algorithm::SolarAlgorithm,
    refraction::RefractionAlgorithm,
)

function SolarPositionBlock(; name)

    @parameters t0::DateTime [tunable = false] observer::Observer [tunable = false]
    @parameters algorithm::SolarAlgorithm = PSA() [tunable = false]
    @parameters refraction::RefractionAlgorithm = NoRefraction() [tunable = false]

    @variables pos(t_nounits) [output = true]

    time_expr = term(seconds_to_datetime, t_nounits, t0; type = DateTime)

    eqs = [pos ~ solar_position(observer, time_expr, algorithm, refraction)]

    return System(eqs, t_nounits; name = name)
end

Now @named sun = SolarPositionBlock() works and gives:

Model sun:
Equations (1):
  1 standard: see equations(sun)
Unknowns (1): see unknowns(sun)
  pos(t)
Parameters (4): see parameters(sun)
  observer
  algorithm [defaults to PSA(2020)]
  t0
  refraction [defaults to NoRefraction()]Variables: SymbolicUtils.BasicSymbolic[]

I don’t see how to get to azimuth and elevation from here.

7

After struggling a little bit more this seems to work beautifully now. Some critical feedback to the efficiency of this method would be greatly appreciated.


using SolarPosition: Observer, SolarAlgorithm, RefractionAlgorithm, PSA, NoRefraction
using ModelingToolkit: @parameters, @variables, System, @register_symbolic, t_nounits
using Symbolics
using Symbolics: term
using Dates: DateTime, Millisecond

import SolarPosition: SolarPositionBlock, solar_position


seconds_to_datetime(t_sec, t0::DateTime) = t0 + Millisecond(round(Int, t_sec * 1e3))

# Helper functions to extract fields from solar position
get_azimuth(pos) = pos.azimuth
get_elevation(pos) = pos.elevation
get_zenith(pos) = pos.zenith

@register_symbolic seconds_to_datetime(t_sec, t0::DateTime)
@register_symbolic solar_position(
    observer::Observer,
    time::DateTime,
    algorithm::SolarAlgorithm,
    refraction::RefractionAlgorithm,
)
@register_symbolic get_azimuth(pos)
@register_symbolic get_elevation(pos)
@register_symbolic get_zenith(pos)

function SolarPositionBlock(; name)

    @parameters t0::DateTime [tunable = false] observer::Observer [tunable = false]
    @parameters algorithm::SolarAlgorithm = PSA() [tunable = false]
    @parameters refraction::RefractionAlgorithm = NoRefraction() [tunable = false]

    @variables azimuth(t_nounits) [output = true]
    @variables elevation(t_nounits) [output = true]
    @variables zenith(t_nounits) [output = true]

    time_expr = term(seconds_to_datetime, t_nounits, t0; type = DateTime)
    pos = solar_position(observer, time_expr, algorithm, refraction)

    eqs = [
        azimuth ~ get_azimuth(pos),
        elevation ~ get_elevation(pos),
        zenith ~ get_zenith(pos),
    ]

    return System(eqs, t_nounits; name = name)
end

which I can then put in an ODEProblem like this:

@named sun = SolarPositionBlock()

sys = mtkcompile(sun)

pmap = [
    sys.observer => obs,
    sys.t0 => t0,
    sys.algorithm => PSA(),
    sys.refraction => NoRefraction(),
]

# create ODEProblem for 24 hours (86400 seconds)
tspan = (0.0, 86400.0)
prob = ODEProblem(sys, pmap, tspan)

sol = solve(prob; saveat = 60.0)

image
image2400×1600 236 KB

8

I’m not sure if we can avoid the explicit term. I think this should be fine.

1 Like