Expose a non dummy derivative index lowering

[ChrisRackauckas]

Fixes #2513

[baggepinnen]

Is this supposed to be user facing?

[ChrisRackauckas]

To an extent, because it has different trade-offs in the system that is generated. However, the non-DD version seems to be a bit buggy so I'm not sure if it's ready.

[ChrisRackauckas]

MWE:

using ModelingToolkit, OrdinaryDiffEq, Test
using ModelingToolkit: t_nounits as t, D_nounits as D
@parameters g
@variables x(t) [state_priority = 10] y(t) λ(t)

eqs = [
       D(D(x)) ~ λ * x
       D(D(y)) ~ λ * y - g
       x^2 + y^2 ~ 1
      ]
@named pend = ODESystem(eqs,t)
pend = complete(pend)
ss = structural_simplify(pend; dummy_derivative = false)

unknowns(ss)
equations(ss)

#=
julia> unknowns(ss)
7-element Vector{SymbolicUtils.BasicSymbolic{Real}}:
 x(t)
 y(t)
 xˍtt(t)
 yˍtt(t)
 λ(t)
 xˍt(t)
 yˍt(t)

julia> equations(ss)
5-element Vector{Equation}:
 Differential(t)(x(t)) ~ xˍt(t)
 Differential(t)(y(t)) ~ yˍt(t)
 0 ~ -xˍtt(t) + x(t)*λ(t)
 0 ~ -g - yˍtt(t) + y(t)*λ(t)
 0 ~ 1 - (x(t)^2) - (y(t)^2)
=#

@named pend = ODESystem(eqs,t)
ss = complete(dae_index_lowering(expand_connections(pend)))

#=
julia> equations(ss)
3-element Vector{Equation}:
 Differential(t)(Differential(t)(x(t))) ~ x(t)*λ(t)
 Differential(t)(Differential(t)(y(t))) ~ -g + y(t)*λ(t)
 0 ~ -2(Differential(t)(y(t))^2) - 2(Differential(t)(x(t))^2) - 2(x(t)^2)*λ(t) - 2y(t)*(-g + y(t)*λ(t))
 =#

 ss = complete(ode_order_lowering(dae_index_lowering(expand_connections(pend))))

 #=
 julia> equations(ss)
5-element Vector{Equation}:
 Differential(t)(xˍt(t)) ~ x(t)*λ(t)
 Differential(t)(yˍt(t)) ~ -g + y(t)*λ(t)
 Differential(t)(x(t)) ~ xˍt(t)
 Differential(t)(y(t)) ~ yˍt(t)
 0 ~ -2(Differential(t)(y(t))^2) - 2(Differential(t)(x(t))^2) - 2(x(t)^2)*λ(t) - 2y(t)*(-g + y(t)*λ(t))
 =#

It seems just doing dummy_derivative=false isn't what I expected. I do think as a debugging tool we should expose the last to the user somehow in a documented way. @YingboMa how would you prefer that be done?

[YingboMa]

We can make structural_simplify(sys, dummy_derivative = false) behave like structural_simplify(dae_index_lowering(expand_connections(sys))).

[ChrisRackauckas]

The manual order lowering is required? At least from my tests there. Is there a way to make pantelides do that part?

[YingboMa]

I am not sure if ode_order_lowering(dae_index_lowering(expand_connections(pend))) always return a valid semi-implicit ODE in general.

[YingboMa]

The manual order lowering is required? At least from my tests there. Is there a way to make pantelides do that part?

I think the only reliable way of lowering to first order is structural_simplify.

[ChrisRackauckas]

That's what I tried first, but the structural_simplify way (this PR) gives an invalid system (the code in this PR), so something went wrong 😅 .

[ChrisRackauckas]

Silly is fine if it's correct 😅. It's a bit odd to implement in practice because I assume it needs to be after the validity checks, but then there's state in the tearing state so I reconstruct it. Let me know if this state reconstruction should have another step in there.

[YingboMa]

I don’t think it’s valid to keep the state. I can take a look tomorrow.

[ChrisRackauckas]

I don't keep the state, I generate a new one.

[YingboMa]

Ah, I meant mm.