I am trying to use DiffEq and SciML sensitivity to build NN controllers for physical systems.

Usually everything works, however when i explicitly call `adjoint_sensitivities`

from SciML-sensitivity to get the gradient of a loss w.r.t. the parameters of the system with an algorithm for stiff equations it fails.

For ODEs I can work around it using a higher order algorithm such as Vern9, not ideal but works.

For DAEs (ODEs with Mass-matrix in my case) I don’t see a solution as such algorithms are not compatible with them.

I checked the suggestions but:

- Autodiff is necessary with NNs in my case.
- I don’t have any allocation for cache or similar
- Not possible.

Error:

```
ERROR: First call to automatic differentiation for the Jacobian
failed. This means that the user `f` function is not compatible
with automatic differentiation. Methods to fix this include:
1. Turn off automatic differentiation (e.g. Rosenbrock23() becomes
Rosenbrock23(autodiff=false)). More details can befound at
https://docs.sciml.ai/DiffEqDocs/stable/features/performance_overloads/
2. Improving the compatibility of `f` with ForwardDiff.jl automatic
differentiation (using tools like PreallocationTools.jl). More details
can be found at https://docs.sciml.ai/DiffEqDocs/stable/basics/faq/#Autodifferentiation-and-Dual-Numbers
3. Defining analytical Jacobians. More details can be
found at https://docs.sciml.ai/DiffEqDocs/stable/types/ode_types/#SciMLBase.ODEFunction
Note: turning off automatic differentiation tends to have a very minimal
performance impact (for this use case, because it's forward mode for a
square Jacobian. This is different from optimization gradient scenarios).
However, one should be careful as some methods are more sensitive to
accurate gradients than others. Specifically, Rodas methods like `Rodas4`
and `Rodas5P` require accurate Jacobians in order to have good convergence,
while many other methods like BDF (`QNDF`, `FBDF`), SDIRK (`KenCarp4`),
and Rosenbrock-W (`Rosenbrock23`) do not. Thus if using an algorithm which
is sensitive to autodiff and solving at a low tolerance, please change the
algorithm as well.
TypeError: in typeassert, expected Float64, got a value of type ForwardDiff.Dual{Nothing, Float64, 12}
```

As a MWE the example from the documentation can be used:

```
using DifferentialEquations
using SciMLSensitivity
function rober(du, u, p, t)
y₁, y₂, y₃ = u
k₁, k₂, k₃ = p
du[1] = -k₁ * y₁ + k₃ * y₂ * y₃
du[2] = k₁ * y₁ - k₃ * y₂ * y₃ - k₂ * y₂^2
du[3] = y₁ + y₂ + y₃ - 1
nothing
end
M = [1.0 0 0
0 1.0 0
0 0 0]
f = ODEFunction(rober, mass_matrix = M)
prob_mm = ODEProblem(f, [1.0, 0.0, 0.0], (0.0, 1e5), [0.04, 3e7, 1e4])
sol = solve(prob_mm, Rodas5(), reltol = 1e-8, abstol = 1e-8)
sensealg = InterpolatingAdjoint(autojacvec=ZygoteVJP(), checkpointing=true)
_, dp = adjoint_sensitivities(sol, Rosenbrock23(); sensealg=sensealg, g=(u, p, t)->sum(u));
```

Is there an algorithm that works here?

Any idea what’s going on? Is sensitivity impossibile for DAE at the moment?

DifferentialEquations version = “7.7.0”

SciMLSensitivity version = “7.20.0”