SingularException with Variable constrained not to change.

[strangeli]

If we constraint a variable not to change, the mass matrix based solvers fail with a SingularException. This came up in a design for composing open loop and controller dynamics. Example:

function tf!(du, u, p, t)
  du[1] = 0.
end

tof = ODEFunction(tf!, mass_matrix=zeros(1,1))
tp = ODEProblem(tof, [2.], (0., 5.))
solve(tp, Rodas4())

[strangeli]

Error is

SingularException(1)
checknonsingular at factorization.jl:19 [inlined]
#lu!#3(::Bool, ::Int64, ::typeof(RecursiveFactorization.lu!), ::Array{Float64,2}, ::Array{Int64,1}, ::Val{true}) at lu.jl:32
#lu!#2 at none:0 [inlined]
lu! at lu.jl:10 [inlined]
lu! at lu.jl:10 [inlined]
#_#407(::Nothing, ::Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}}, ::DefaultLinSolve, ::Array{Float64,1}, ::Array{Float64,2}, ::Array{Float64,1}, ::Bool) at linear_nonlinear.jl:62
DefaultLinSolve at linear_nonlinear.jl:53 [inlined]
perform_step!(::OrdinaryDiffEq.ODEIntegrator{Rodas4{0,true,DefaultLinSolve,DataType},true,Array{Float64,1},Float64,DiffEqBase.NullParameters,Float64,Float64,Float64,Array{Array{Float64,1},1},ODESolution{Float64,2,Array{Array{Float64,1},1},Nothing,Nothing,Array{Float64,1},Array{Array{Array{Float64,1},1},1},ODEProblem{Array{Float64,1},Tuple{Float64,Float64},true,DiffEqBase.NullParameters,ODEFunction{true,typeof(tf!),Array{Float64,2},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Base.Iterators.Pairs{Union{},Union{},Tuple{},NamedTuple{(),Tuple{}}},DiffEqBase.StandardODEProblem},Rodas4{0,true,DefaultLinSolve,DataType},OrdinaryDiffEq.InterpolationData{ODEFunction{true,typeof(tf!),Array{Float64,2},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Array{Float64,1},1},Array{Float64,1},Array{Array{Array{Float64,1},1},1},OrdinaryDiffEq.Rodas4Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.RodasTableau{Float64,Float64},DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(tf!),Array{Float64,2},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},DiffEqBase.NullParameters},DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(tf!),Array{Float64,2},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},DefaultLinSolve,SparseDiffTools.ForwardColorJacCache{Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(tf!),Array{Float64,2},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64}(),Float64,1},1},Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(tf!),Array{Float64,2},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64}(),Float64,1},1},Array{Float64,1},Array{Array{Tuple{Bool},1},1},UnitRange{Int64},Nothing},ForwardDiff.DerivativeConfig{ForwardDiff.Tag{DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(tf!),Array{Float64,2},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},DiffEqBase.NullParameters},Float64},Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(tf!),Array{Float64,2},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},DiffEqBase.NullParameters},Float64},Float64,1},1}}}},DiffEqBase.DEStats},ODEFunction{true,typeof(tf!),Array{Float64,2},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},OrdinaryDiffEq.Rodas4Cache{Array{Float64,1},Array{Float64,1},Array{Float64,1},Array{Float64,2},Array{Float64,2},OrdinaryDiffEq.RodasTableau{Float64,Float64},DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(tf!),Array{Float64,2},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},DiffEqBase.NullParameters},DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(tf!),Array{Float64,2},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},DefaultLinSolve,SparseDiffTools.ForwardColorJacCache{Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(tf!),Array{Float64,2},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64}(),Float64,1},1},Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.UJacobianWrapper{ODEFunction{true,typeof(tf!),Array{Float64,2},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Float64,DiffEqBase.NullParameters},Float64}(),Float64,1},1},Array{Float64,1},Array{Array{Tuple{Bool},1},1},UnitRange{Int64},Nothing},ForwardDiff.DerivativeConfig{ForwardDiff.Tag{DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(tf!),Array{Float64,2},Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing,Nothing},Array{Float64,1},DiffEqBase.NullParameters},Float64},Array{ForwardDiff.Dual{ForwardDiff.Tag{DiffEqDiffTools.TimeGradientWrapper{ODEFunction{true,typeof(tf!),Array{Float6...

[ChrisRackauckas]

What does a mass matrix of zeros mean? Your equation is 0u' = 0 🤷‍♂

[strangeli]

This is a small part of a larger equation. The idea was that another function would eventually modify the right hand side to provide a constraint so we obtain 0u' = f(u). I had expected that the solver picks u' = 0 if no f(u) is provided, but I see now that this is truly underspecified, and thus the error is appropriate. I'll close the issue.

Our solution indeed was to modify it to 0u' = u and require the controller to provide f(u) - u.

[ChrisRackauckas]

There might be an issue to work through, but this example was too simplified to work out, because it truly has infinitely many solutions. You'd need to have something like, 2 variables, 1 ODE, and 1 constraint term, and then it's something that's satisfiable. That could possibly be the issue, that there are too many constraints so the problem is ill-posed.