Index reduction to model DAE as ODE

[taylormcd]

DAEs can be converted to ODEs using index reduction. It would be nice to try that for this package so that more solvers/features from DifferentialEquations can be used. The sensitivity analysis in DifferentialEquations doesn't cover DAEs, so that would be an example of at least one area in which this would be useful.

[ChrisRackauckas]

The sensitivity analysis does cover index-1 DAEs in mass-matrix form. But not fully implicit DAEs at this moment (though that should change). Index reduction can be done with ModelingToolkit.jl

[taylormcd]

@ChrisRackauckas That's good to know. Does the sensitivity analysis also cover non-constant mass matrices created by using DiffEqArrayOperators? I think trying to combine sensitivity analyses and non-constant mass matrices might have caused me to come to a false conclusion about DifferentialEquations' support for DAEs.

On a related note, we're planning on implementing the unsteady adjoint method for DAEs like the one modeled by this package (fully implicit, or at least a non-constant mass matrix) for some of the work we're doing in our lab. So we might be able to help get this capability into DifferentialEquations. We'll probably implement it outside of DifferentialEquations first, but if you have any tips for getting this going let me know.

[ChrisRackauckas]

On a related note, we're planning on implementing the unsteady adjoint method for DAEs like the one modeled by this package (fully implicit, or at least a non-constant mass matrix) for some of the work we're doing in our lab

You mean backsolve adjoint? That is unlikely to be stable for any DAEs: we know that from tests on the mass matrix case. I think the more important thing would be to define the InterpolatingAdjoint and QuadratureAdjoint cases. They would just need the relevant DAEFunction dispatches and the code for handling the reinitialization. Basically just this part:

https://github.com/SciML/DiffEqSensitivity.jl/blob/master/src/interpolating_adjoint.jl#L101-L123

The reinitialization piece of course is the really tricky part of the general DAE adjoint. It might be good to do it first without the reinitialization schemes and then add those in. I wonder if @frankschae is up for the challenge here.

Does the sensitivity analysis also cover non-constant mass matrices created by using DiffEqArrayOperators? I think trying to combine sensitivity analyses and non-constant mass matrices might have caused me to come to a false conclusion about DifferentialEquations' support for DAEs.

It does not. If you look at the supplemental of https://arxiv.org/pdf/2001.04385.pdf, equations 35 and 36 relies on the mass matrix being constant in order for the consistent initialization to result in a linear solve. Otherwise you will need to have the full DAE reinitialization procedure done there, which cannot use the standard DAE reinit phase. Adding the more general reinitialization fallbacks though would be what's needed to fix that (but of course it's a lot more expensive)

[frankschae]

Sounds interesting. I'd be willing to help here -- I don't have much experience with DAEs and the associated adjoints though. So it's probably more of a medium-to-long-term thing I think.

[taylormcd]

You mean backsolve adjoint? That is unlikely to be stable for any DAEs: we know that from tests on the mass matrix case. I think the more important thing would be to define the InterpolatingAdjoint and QuadratureAdjoint cases.

Any of these would probably work. We haven't settled on a specific method, though InterpolatingAdjoint probably works for our use case. I think it would be nice to test out a few options.

The reinitialization piece of course is the really tricky part of the general DAE adjoint. It might be good to do it first without the reinitialization schemes and then add those in.

I'll keep that in mind. I still have a bit of learning to do on the adjoint method, but my time frame is more short-to-medium term so the learning should happen quickly.