Steady state solution in ModelingToolkit fails with SSRootFind due to current_time

While playing with a small example to test different functionalities of MTK it seems that SSRootFind fails to solve for the steady state, contrarily to DynamicSS, on a current_time call.

Here is a MWE to reproduce

using ModelingToolkit, DifferentialEquations, ModelingToolkitStandardLibrary
import ModelingToolkit: connect, t_nounits as t, D_nounits as D
import ModelingToolkitStandardLibrary.Mechanical.TranslationalPosition: Mass, Damper, Spring, Fixed
import ModelingToolkitStandardLibrary.Blocks: Sine
using SymbolicIndexingInterface: parameter_values, state_values
using SciMLStructures: Tunable, canonicalize, replace, replace!
using PreallocationTools

# Build model
@named ds = Sine(amplitude=-1, frequency=10)
@named mass = Mass(m=1)
@named spring = Spring(k=4, l=0.)
@named damper = Damper(d=1)
@named fix = Fixed()

eqs = Equation[
    connect(fix.flange, spring.flange_a, damper.flange_a)
    connect(mass.flange, spring.flange_b, damper.flange_b)
]

@named model = ODESystem(eqs, t; systems=[mass, spring, damper, fix])
sys = mtkcompile(model)


# Define transient problem
initial_positions = Dict{Num,Any}([
    mass.s => 1.
    mass.v => 0.
]
)
prob = ODEProblem(sys, initial_positions, (0.0, 100.0))
sol_mass = solve(prob, Rodas5())

# Define steady state problem
prob_ss = SteadyStateProblem(prob)
sol_ss = solve(prob_ss, SSRootfind(), abstol=1e-8, reltol=1e-8)

The error reads

julia> sol_ss = solve(prob_ss, SSRootfind(), abstol=1e-8, reltol=1e-8)
ERROR: MethodError: no method matching current_time(::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithmCache{…})
The function `current_time` exists, but no method is defined for this combination of argument types.

Closest candidates are:
  current_time(::Any, ::Colon)
   @ SymbolicIndexingInterface D:\bamboo\.julia\packages\SymbolicIndexingInterface\hhVqD\src\value_provider_interface.jl:143
  current_time(::Any, ::Any)
   @ SymbolicIndexingInterface D:\bamboo\.julia\packages\SymbolicIndexingInterface\hhVqD\src\value_provider_interface.jl:142
  current_time(::SciMLBase.AbstractJumpProblem)
   @ SciMLBase D:\bamboo\.julia\packages\SciMLBase\hHrRi\src\problems\problem_interface.jl:25
  ...

Stacktrace:
  [1] (::ModelingToolkit.ObservedWrapper{…})(prob::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithmCache{…})
    @ ModelingToolkit D:\bamboo\.julia\packages\ModelingToolkit\Ay2JZ\src\systems\problem_utils.jl:695
  [2] call_composed(fs::Tuple{…}, x::Tuple{…}, kw::@Kwargs{})
    @ Base .\operators.jl:1054
  [3] call_composed
    @ .\operators.jl:1053 [inlined]
  [4] (::ComposedFunction{…})(x::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithmCache{…}; kw::@Kwargs{})
    @ Base .\operators.jl:1050
  [5] (::ModelingToolkit.var"#_getter#952"{…})(valp::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithmCache{…}, initprob::NonlinearProblem{…})
    @ ModelingToolkit D:\bamboo\.julia\packages\ModelingToolkit\Ay2JZ\src\systems\problem_utils.jl:846
  [6] update_initializeprob!(initprob::NonlinearProblem{…}, prob::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithmCache{…})
    @ ModelingToolkit D:\bamboo\.julia\packages\ModelingToolkit\Ay2JZ\src\systems\problem_utils.jl:981
  [7] get_initial_values(prob::NonlinearProblem{…}, valp::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithmCache{…}, f::NonlinearFunction{…}, alg::OverrideInit{…}, iip::Val{…}; nlsolve_alg::NonlinearSolvePolyAlgorithm{…}, abstol::Float64, reltol::Float64, kwargs::@Kwargs{})
    @ SciMLBase D:\bamboo\.julia\packages\SciMLBase\hHrRi\src\initialization.jl:259
  [8] _run_initialization!(cache::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithmCache{…}, initalg::OverrideInit{…}, prob::NonlinearProblem{…}, isinplace::Bool)
    @ NonlinearSolveBase D:\bamboo\.julia\packages\NonlinearSolveBase\ejkqA\src\initialization.jl:26
  [9] _run_initialization!(cache::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithmCache{…}, ::NonlinearSolveBase.NonlinearSolveDefaultInit, prob::NonlinearProblem{…}, isinplace::Bool)
    @ NonlinearSolveBase D:\bamboo\.julia\packages\NonlinearSolveBase\ejkqA\src\initialization.jl:10
 [10] run_initialization!(cache::NonlinearSolveFirstOrder.GeneralizedFirstOrderAlgorithmCache{…}, initializealg::NonlinearSolveBase.NonlinearSolveDefaultInit, prob::NonlinearProblem{…})
    @ NonlinearSolveBase D:\bamboo\.julia\packages\NonlinearSolveBase\ejkqA\src\initialization.jl:4
 [11] __init(::NonlinearProblem{…}, ::GeneralizedFirstOrderAlgorithm{…}; stats::SciMLBase.NLStats, alias_u0::Bool, maxiters::Int64, abstol::Float64, reltol::Float64, maxtime::Nothing, termination_condition::Nothing, internalnorm::typeof(NonlinearSolveBase.L2_NORM), verbose::NonlinearVerbosity{…}, linsolve_kwargs::@NamedTuple{}, initializealg::NonlinearSolveBase.NonlinearSolveDefaultInit, kwargs::@Kwargs{})
    @ NonlinearSolveFirstOrder D:\bamboo\.julia\packages\NonlinearSolveFirstOrder\inFag\src\solve.jl:238
 [12] __solve(::NonlinearProblem{…}, ::GeneralizedFirstOrderAlgorithm{…}; kwargs::@Kwargs{…})
    @ NonlinearSolveBase D:\bamboo\.julia\packages\NonlinearSolveBase\ejkqA\src\solve.jl:269
 [13] macro expansion
    @ D:\bamboo\.julia\packages\NonlinearSolveBase\ejkqA\src\solve.jl:475 [inlined]
 [14] __generated_polysolve(::NonlinearProblem{…}, ::NonlinearSolvePolyAlgorithm{…}; stats::SciMLBase.NLStats, alias::SciMLBase.NonlinearAliasSpecifier, verbose::NonlinearVerbosity{…}, initializealg::NonlinearSolveBase.NonlinearSolveDefaultInit, kwargs::@Kwargs{…})
    @ NonlinearSolveBase D:\bamboo\.julia\packages\NonlinearSolveBase\ejkqA\src\solve.jl:420
 [15] __solve(::NonlinearProblem{…}, ::NonlinearSolvePolyAlgorithm{…}; kwargs::@Kwargs{…})
    @ NonlinearSolveBase D:\bamboo\.julia\packages\NonlinearSolveBase\ejkqA\src\solve.jl:393
 [16] __solve(::NonlinearProblem{…}, ::Nothing; kwargs::@Kwargs{…})
    @ NonlinearSolve D:\bamboo\.julia\packages\NonlinearSolve\JTN4J\src\default.jl:27
 [17] solve_call(_prob::NonlinearProblem{…}, args::Nothing; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ NonlinearSolveBase D:\bamboo\.julia\packages\NonlinearSolveBase\ejkqA\src\solve.jl:168
 [18] solve_up(prob::NonlinearProblem{…}, sensealg::Nothing, u0::Vector{…}, p::MTKParameters{…}, args::Nothing; originator::SciMLBase.ChainRulesOriginator, kwargs::@Kwargs{…})
    @ NonlinearSolveBase D:\bamboo\.julia\packages\NonlinearSolveBase\ejkqA\src\solve.jl:116
 [19] solve(prob::NonlinearProblem{…}, args::Nothing; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, verbose::NonlinearVerbosity{…}, kwargs::@Kwargs{…})
    @ NonlinearSolveBase D:\bamboo\.julia\packages\NonlinearSolveBase\ejkqA\src\solve.jl:87
 [20] __solve(::SteadyStateProblem{…}, ::SSRootfind{…}; kwargs::@Kwargs{…})
    @ SteadyStateDiffEq D:\bamboo\.julia\packages\SteadyStateDiffEq\FbOma\src\solve.jl:4
 [21] solve_call(_prob::SteadyStateProblem{…}, args::SSRootfind{…}; merge_callbacks::Bool, kwargshandle::Nothing, kwargs::@Kwargs{…})
    @ NonlinearSolveBase D:\bamboo\.julia\packages\NonlinearSolveBase\ejkqA\src\solve.jl:168
 [22] #solve_up#147
    @ D:\bamboo\.julia\packages\NonlinearSolveBase\ejkqA\src\solve.jl:116 [inlined]
 [23] solve_up
    @ D:\bamboo\.julia\packages\NonlinearSolveBase\ejkqA\src\solve.jl:109 [inlined]
 [24] solve(prob::SteadyStateProblem{…}, args::SSRootfind{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, verbose::NonlinearVerbosity{…}, kwargs::@Kwargs{…})
    @ NonlinearSolveBase D:\bamboo\.julia\packages\NonlinearSolveBase\ejkqA\src\solve.jl:87
Some type information was truncated. Use `show(err)` to see complete types.

This is with

[0c46a032] DifferentialEquations v7.17.0
[961ee093] ModelingToolkit v10.31.2
[16a59e39] ModelingToolkitStandardLibrary v2.25.0
[8913a72c] NonlinearSolve v4.12.0
[9672c7b4] SteadyStateDiffEq v2.9.0

Open an issue for this. This is just a bug in the codegen.

1 Like

Will do. Bug issue here Steady state solution in ModelingToolkit fails with SSRootFind due to current_time · Issue #4172 · SciML/ModelingToolkit.jl · GitHub