Input Differential in Observables

[bradcarman]

Take the following...

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D

pars = @parameters begin
    k = 1
end

vars = @variables begin
    x(t) = 0
    dx(t) = 0
    ddx(t) = 0

    y(t)
    dy(t)
    ddy(t)

    input(t) = 0, [input=true]
end

eqs = [
    D(x) ~ dx
    D(dx) ~ ddx

    D(y) ~ dy
    D(dy) ~ ddy

    ddy ~ k*(y - x)

    D(ddy) ~ input
]

@named sys = ODESystem(eqs, t, vars, pars)
sys, = structural_simplify(sys, ([input], Num[]))

This gives the following equations and observables...

julia> full_equations(sys)
3-element Vector{Equation}:
 Differential(t)(ddy(t)) ~ input(t)
 Differential(t)(dy(t)) ~ ddy(t)
 Differential(t)(y(t)) ~ dy(t)

julia> observed(sys)
8-element Vector{Equation}:
 x(t) ~ (-ddy(t) + k*y(t)) / k
 ddyˍtt(t) ~ Differential(t)(input(t))
 yˍtt(t) ~ ddy(t)
 xˍt(t) ~ (-input(t) + k*dy(t)) / k
 xˍtt(t) ~ (-ddyˍtt(t) + k*yˍtt(t)) / k
 dx(t) ~ xˍt(t)
 dxˍt(t) ~ xˍtt(t)
 ddx(t) ~ dxˍt(t)

Note: the term ddyˍtt(t) ~ Differential(t)(input(t)). Inputs should be considered piecewise constant and therefore the derivative should be assumed zero. Or if that assumption is not acceptable, there should be a way to then provide the definition of D(input).

This issue expands on the requirements laid out in #3823

[ChrisRackauckas]

It should be handled through the same mechanism as #3493

[ChrisRackauckas]

Basically finishing #2529, i.e. lowering inputs to discrete variables / time-dependent parameters

[AayushSabharwal]

Just from looking at the result, it looks like MTK differentiates ddy ~ k*(y - x) twice to be able to find a value for ddx.

The first differentiation yields D(ddy) ~ k * (D(y) - D(x)) which is just input ~ k * (dy - dx) which it rearranges to (essentially) dx ~ (k * dy - input) / k. The second differentiation yields D(D(ddy)) ~ k * (D(dy) - D(dx)) which is ddy_t ~ k * (ddy - ddx) which it rearranges to ddx ~ (ddy * k - ddy_t) / k. ddy_t is just D(D(ddy)) which is D(input). It does not know D(input). I'm not sure if we can provide D(input) because that processing happens before the inputs are handled.

[ChrisRackauckas]

Such inputs are assumed piecewise constant in cosimulation

[bradcarman]

Moving to v10 I can confirm this issue is not resolved...

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using Test

pars = @parameters begin
    k = 1
end

vars = @variables begin
    x(t) = 0
    dx(t) = 0
    ddx(t) = 0

    y(t)
    dy(t)
    ddy(t)

    input(t) = 0, [input=true]
end

eqs = [
    D(x) ~ dx
    D(dx) ~ ddx

    D(y) ~ dy
    D(dy) ~ ddy

    ddy ~ k*(y - x)

    D(ddy) ~ input
]

@mtkcompile sys = System(eqs, t, vars, pars) inputs=[input]

vs = Set()
for eq in equations(sys)
    ModelingToolkit.vars!(vs, eq)
end
for eq in observed(sys)
    ModelingToolkit.vars!(vs, eq)
end

@test !(D(sys.input) in vs)

Gives...

Test Failed at c:\Work\Packages\CatapultModels\test\mwe1.jl:43
  Expression: !(D(sys.input) in vs)
   Evaluated: !(Differential(t)(input(t)) in Set(Any[Differential(t)(ddy(t)), yˍtt(t), ddx(t), input(t), ddy(t), k, xˍtt(t), y(t), ddyˍtt(t), dx(t), x(t), Differential(t)(input(t)), dxˍt(t), Differential(t)(y(t)), dy(t), xˍt(t), Differential(t)(dy(t))]))

[ChrisRackauckas]

Yes #2529 didn't merge. I know what needs to be done here but it just didn't complete because it needed MTK v10 so we can swap states and parameters in the compiler without having it metadata.