ModelingToolkit broken for SteadyStateProblems

[jd-lara]

There is no way of solving SteadyState problems with ModelingToolKit. Here are two ways I tried and they return different errors

MWE

using ModelingToolkit
using DiffEqBase
using SteadyStateDiffEq
using LinearAlgebra

@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t

eqs = [D(D(x)) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

sys = ODESystem(eqs)
sys = ode_order_lowering(sys)

u0 = [D(x) => 2.0,
      x => 1.0,
      y => 0.0,
      z => 0.0]

p  = [σ => 28.0,
      ρ => 10.0,
      β => 8/3]

First try

prob =  SteadyStateProblem(sys, u0, p, mass_matrix=LinearAlgebra.I)
solve(prob, SSRootfind())

julia> solve(prob, SSRootfind())
ERROR: MethodError: no method matching Pair{Operation,Float64}(::Float64)
Closest candidates are:
  Pair{Operation,Float64}(::Any, ::Any) where {A, B} at pair.jl:7
Stacktrace:
 [1] alloc_DF(::Array{Pair{Operation,Float64},1}, ::Array{Pair{Operation,Float64},1}) at /Users/jdlara/.julia/packages/NLSolversBase/mGaJg/src/objective_types/abstract.jl:19
 [2] nlsolve(::Function, ::Array{Pair{Operation,Float64},1}; method::Symbol, autodiff::Symbol, inplace::Bool, kwargs::Base.Iterators.Pairs{Symbol,Float64,Tuple{Symbol},NamedTuple{(:ftol,),Tuple{Float64}}}) at /Users/jdlara/.julia/packages/NLsolve/ZBTu4/src/nlsolve/nlsolve.jl:49
 [3] (::SteadyStateDiffEq.var"#2#4")(::Function, ::Array{Pair{Operation,Float64},1}, ::Float64) at /Users/jdlara/.julia/packages/SteadyStateDiffEq/q0UMK/src/algorithms.jl:6
 [4] __solve(::SteadyStateProblem{Array{Pair{Operation,Float64},1},true,Array{Pair{Operation,Float64},1},ODEFunction{true,ModelingToolkit.var"#f#84"{GeneralizedGenerated.NGG.RuntimeFn{TypeEncoding(list(##MTKArg#255, ##MTKArg#256, ##MTKArg#257)),TypeEncoding(nil(GeneralizedGenerated.NGG.Argument)),TypeEncoding(begin

Second try

prob = ODEProblem(sys,u0,(0.0,10.0),p)
ss_prob = SteadyStateProblem(prob)
solve(prob, SSRootfind())

ERROR: Default algorithm choices require DifferentialEquations.jl. Please specify an algorithm or import DifferentialEquations directly.

[ovile27]

In the second try you solve for prob, but it should be

solve(ss_prob, SSRootfind())

This works for me:

using ModelingToolkit
using DifferentialEquations

@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t

eqs = [D(D(x)) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]

sys = ODESystem(eqs)
sys = ode_order_lowering(sys)

u0 = [D(x) => 2.0,
      x => 1.0,
      y => 0.0,
      z => 0.0]

p  = [σ => 28.0,
      ρ => 10.0,
      β => 8/3]

tspan = [0.0, 10.0]

prob = ODEProblem(sys,u0,tspan,p;jac=true,sparse=true)
ssprob = SteadyStateProblem(prob)
sol = solve(ssprob,SSRootfind())

Zero solution

julia> sol = solve(ssprob,SSRootfind())
u: 4-element Array{Float64,1}:
  2.985167668612121e-12
  8.18789480661053e-12
  7.93898280448957e-12
 -6.8833827526759706e-15

u0 = [D(x) => 2.0,
      x => 7,
      y => 7,
      z => 20]

prob = ODEProblem(sys,u0,tspan,p)
ssprob = SteadyStateProblem(prob)
sol = solve(ssprob,SSRootfind())

Non-Zero solution

julia> sol = solve(ssprob,SSRootfind())
u: 4-element Array{Float64,1}:
 0.0
 4.89897948556666
 9.000000000000574
 4.89897948556666

However the order of x y z is changed since it should be [x,y,z] = [4.899, 4.899, 9].

Julia Version 1.3.1
DifferentialEquations v6.13.0
ModelingToolkit v3.1.1

[ChrisRackauckas]

Fixed