Fix Symbolics derivative wrapping for structural simplification

[ChrisRackauckas]

Should fix SciML/ModelingToolkit.jl#2190

@bradcarman can you run your test case with this?

[bradcarman]

Unfortunately this still fails. If I run the following at the end of derivative_tests.jl I still get ERROR: Differentiation with array expressions is not yet supported

f = LinearInterpolation(u, t)
expand_derivatives(D(f))

If I understand correctly this is what ModelingToolkit will be doing

[bradcarman]

MWE

    using ModelingToolkit
    using ModelingToolkitStandardLibrary.Blocks: TimeVaryingFunction
    using DataInterpolations

    dt = 4e-4
    t_end = 10.0
    time = 0:dt:t_end
    x = @. time^2 + 1.0

    vars = @variables y(t)=1 dy(t)=0 ddy(t)=0
    @named src = TimeVaryingFunction(LinearInterpolation(x, time))
    @named int = Integrator()
    @named iosys = ODESystem([y ~ src.output.u
            D(y) ~ dy
            D(dy) ~ ddy
            connect(src.output, int.input)],
        t,
        systems = [int, src])
    sys = structural_simplify(iosys)

This gives: ERROR: Differentiation with array expressions is not yet supported

However...

julia> Symbolics.expand_derivatives(D(equations(iosys)[end].rhs))
DataInterpolations.derivative([1.0, 1.00000016, 1.00000064, 1.00000144, 1.00000256, 1.000004, 1.00000576, 1.00000784, 1.00001024, 1.00001296  …  9.9964, 9.9968, 9.9972, 9.9976, 9.998, 9.9984, 9.9988, 9.9992, 9.9996, 10.0], 
t)

So it seems somehow structural_simplify is calling expand_derivatives incorrectly?

[bradcarman]

OK, I think I understand now what's happening. In this case structural_simplify is attempting to expand the 2nd derivative. Therefore, the following model will work when D(dy) ~ ddy is commented out, and error when it's included...

using ModelingToolkit
using ModelingToolkitStandardLibrary.Blocks: TimeVaryingFunction, Integrator
using DataInterpolations

@parameters t
D = Differential(t)

dt = 4e-4
t_end = 10.0
time = 0:dt:t_end
x = @. time^2 + 1.0

vars = @variables y(t)=1 dy(t)=0 ddy(t)=0
@named src = TimeVaryingFunction(LinearInterpolation(x, time))
@named int = Integrator()
@named sys = ODESystem([y ~ src.output.u
        D(y) ~ dy
        # D(dy) ~ ddy # <-- including this gives the error
        connect(src.output, int.input)],
    t,
    systems = [int, src])
sys = structural_simplify(sys)

This error can be recreated simply with...

u = [0.0, 1.5, 0.0]
t = [0.0, 0.5, 1.0]
@variables τ
D = Symbolics.Differential(τ)
f = LinearInterpolation(u, t)

df = expand_derivatives(D(f(τ))) # this is OK
ddf = expand_derivatives(D(D(f(τ)))) # ERROR: Differentiation with array expressions is not yet supported

[ChrisRackauckas]

Okay so at least this work with first derivatives though?

[SebastianM-C]

@bradcarman I just tried

using ModelingToolkit
using ModelingToolkitStandardLibrary.Blocks: TimeVaryingFunction, Integrator
using DataInterpolations

@parameters t
D = Differential(t)

dt = 4e-4
t_end = 10.0
time = 0:dt:t_end
x = @. time^2 + 1.0

vars = @variables y(t) = 1 dy(t) = 0 ddy(t) = 0
@named src = TimeVaryingFunction(LinearInterpolation(x, time))
@named int = Integrator()
@named sys = ODESystem([y ~ src.output.u,
              D(y) ~ dy,
              connect(src.output, int.input)],
       t,
       systems=[int, src])
sys = structural_simplify(sys)

with the patch from this PR applied and I'm still getting the

ERROR: Differentiation with array expressions is not yet supported

even if

g = LinearInterpolation(x, time)
expand_derivatives(D(g(t)))

works. Am I missing something?

[bradcarman]

Hi @SebastianM-C , actually I can confirm the same error as you. I suspect structural_simplify is doing something different now. Maybe a second dummy derivative is being generated?

[ChrisRackauckas]

@sathvikbhagavan can you put the final touches on this so we can get it merged?

[sathvikbhagavan]

@ChrisRackauckas I was looking into this and was having some issues in defining second order derivatives in symbolics. I want to define a function like https://github.com/SciML/DataInterpolations.jl/blob/ChrisRackauckas-patch-1/ext/DataInterpolationsSymbolicsExt.jl#L22 for second order which then can be lowered to computing the values using ForwardDiff numerically.
What is the right way to do it?

[ChrisRackauckas]

Let's just get first derivatives merged first.