Type promotion error when differentiating SteadyStateProblem

[SebastianM-C]

Describe the bug 🐞

A clear and concise description of what the bug is.

Expected behavior
It looks like the propagation of duals through a SteadyStateProblem hits some type promotion issue in initialization.

Minimal Reproducible Example 👇

Without MRE, we would only be able to help you to a limited extent, and attention to the issue would be limited. to know more about MRE refer to wikipedia and stackoverflow.

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using DifferentiationInterface
using SciMLBase
using SymbolicIndexingInterface
using NonlinearSolve

function reactionsystem()
    @variables s1(t) = 2.0 s1s2(t) = 2.0 s2(t) = 2.0
    @parameters k1 = 1.0 c1 = 2.0
    eqs = [D(s1) ~ -0.25 * c1 * k1 * s1 * s2
        D(s1s2) ~ 0.25 * c1 * k1 * s1 * s2
        D(s2) ~ -0.25 * c1 * k1 * s1 * s2]

    return structural_simplify(ODESystem(eqs, t; name=:sys))
end

sys = reactionsystem()
prob = SteadyStateProblem{true}(sys, [k1 => 1.5])
sol = solve(prob, FastShortcutNonlinearPolyalg())

data = sol.u
oop_update = setsym_oop(prob, [sys.k1])

function loss(x, opt_ps)
    prob, oop_update, data = opt_ps

    u0, p = oop_update(prob, x)
    new_prob = remake(prob; u0, p)

    new_sol = solve(new_prob, FastShortcutNonlinearPolyalg())

    !SciMLBase.successful_retcode(new_sol) && return Inf
    sum(abs.(new_sol .- data))
end

opt_ps = (prob, oop_update, data);

of = OptimizationFunction(loss, AutoForwardDiff())
op = OptimizationProblem(of, [1.1], opt_ps)
solve(op, Optimization.LBFGS())

Error & Stacktrace ⚠️

ERROR: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{DifferentiationInterface.FixTail{…}, Float64}, Float64, 1})
The type `Float64` exists, but no method is defined for this combination of argument types when trying to construct it.

Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat
   @ Base rounding.jl:265
  (::Type{T})(::T) where T<:Number
   @ Core boot.jl:900
  Float64(::IrrationalConstants.Inv4π)
   @ IrrationalConstants ~/.julia/packages/IrrationalConstants/lWTip/src/macro.jl:131
  ...

Stacktrace:
  [1] convert(::Type{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{DifferentiationInterface.FixTail{…}, Float64}, Float64, 1})
    @ Base ./number.jl:7
  [2] setindex!(A::Vector{Float64}, x::ForwardDiff.Dual{ForwardDiff.Tag{…}, Float64, 1}, i::Int64)
    @ Base ./array.jl:987
  [3] set_parameter!(p::MTKParameters{…}, val::ForwardDiff.Dual{…}, pidx::ModelingToolkit.ParameterIndex{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/8S2W1/src/systems/parameter_buffer.jl:322
  [4] set_parameter!(sys::NonlinearProblem{…}, val::ForwardDiff.Dual{…}, idx::ModelingToolkit.ParameterIndex{…})
    @ SymbolicIndexingInterface ~/.julia/packages/SymbolicIndexingInterface/gQU1g/src/value_provider_interface.jl:67
  [5] (::SymbolicIndexingInterface.SetParameterIndex{…})(prob::NonlinearProblem{…}, val::ForwardDiff.Dual{…})
    @ SymbolicIndexingInterface ~/.julia/packages/SymbolicIndexingInterface/gQU1g/src/parameter_indexing.jl:678
  [6] (::SymbolicIndexingInterface.ParameterHookWrapper{…})(prob::NonlinearProblem{…}, args::ForwardDiff.Dual{…})
    @ SymbolicIndexingInterface ~/.julia/packages/SymbolicIndexingInterface/gQU1g/src/parameter_indexing.jl:646
  [7] (::SymbolicIndexingInterface.var"#44#45"{…})(s!::SymbolicIndexingInterface.ParameterHookWrapper{…}, v::ForwardDiff.Dual{…})
    @ SymbolicIndexingInterface ~/.julia/packages/SymbolicIndexingInterface/gQU1g/src/parameter_indexing.jl:695
  [8] (::Base.var"#4#5"{…})(a::Tuple{…})
    @ Base ./generator.jl:37
  [9] iterate
    @ ./generator.jl:48 [inlined]
 [10] collect
    @ ./array.jl:791 [inlined]
 [11] map
    @ ./abstractarray.jl:3495 [inlined]
 [12] MultipleSetters
    @ ~/.julia/packages/SymbolicIndexingInterface/gQU1g/src/parameter_indexing.jl:695 [inlined]
 [13] (::ModelingToolkit.UpdateInitializeprob{…})(initializeprob::NonlinearProblem{…}, prob::NonlinearProblem{…})
    @ ModelingToolkit ~/.julia/packages/ModelingToolkit/8S2W1/src/systems/problem_utils.jl:490
 [14] get_initial_values(prob::NonlinearProblem{…}, valp::NonlinearProblem{…}, f::Function, alg::SciMLBase.OverrideInit{…}, iip::Val{…}; nlsolve_alg::NonlinearSolvePolyAlgorithm{…}, abstol::ForwardDiff.Dual{…}, reltol::ForwardDiff.Dual{…}, kwargs::@Kwargs{})
    @ SciMLBase ~/.julia/packages/SciMLBase/sYmAV/src/initialization.jl:235
 [15] get_initial_values
    @ ~/.julia/packages/SciMLBase/sYmAV/src/initialization.jl:222 [inlined]
 [16] _run_initialization!
    @ ~/.julia/packages/NonlinearSolveBase/jA9TW/src/initialization.jl:27 [inlined]
 [17] _run_initialization!
    @ ~/.julia/packages/NonlinearSolveBase/jA9TW/src/initialization.jl:11 [inlined]
 [18] run_initialization!
    @ ~/.julia/packages/NonlinearSolveBase/jA9TW/src/initialization.jl:4 [inlined]
 [19] macro expansion
    @ ~/.julia/packages/NonlinearSolveBase/jA9TW/src/solve.jl:149 [inlined]
 [20] __generated_polysolve(::NonlinearProblem{…}, ::NonlinearSolvePolyAlgorithm{…}; stats::SciMLBase.NLStats, alias_u0::Bool, verbose::Bool, initializealg::NonlinearSolveBase.NonlinearSolveDefaultInit, kwargs::@Kwargs{})
    @ NonlinearSolveBase ~/.julia/packages/NonlinearSolveBase/jA9TW/src/solve.jl:130
 [21] __generated_polysolve
    @ ~/.julia/packages/NonlinearSolveBase/jA9TW/src/solve.jl:130 [inlined]
 [22] #__solve#154
    @ ~/.julia/packages/NonlinearSolveBase/jA9TW/src/solve.jl:127 [inlined]
 [23] __solve
    @ ~/.julia/packages/NonlinearSolveBase/jA9TW/src/solve.jl:124 [inlined]
 [24] #solve_call#35
    @ ~/.julia/packages/DiffEqBase/HGITF/src/solve.jl:635 [inlined]
 [25] solve_call
    @ ~/.julia/packages/DiffEqBase/HGITF/src/solve.jl:592 [inlined]
 [26] #solve_call#45
    @ ~/.julia/packages/DiffEqBase/HGITF/src/solve.jl:1136 [inlined]
 [27] solve_call
    @ ~/.julia/packages/DiffEqBase/HGITF/src/solve.jl:1133 [inlined]
 [28] #solve_up#44
    @ ~/.julia/packages/DiffEqBase/HGITF/src/solve.jl:1112 [inlined]
 [29] solve_up
    @ ~/.julia/packages/DiffEqBase/HGITF/src/solve.jl:1106 [inlined]
 [30] solve(prob::SteadyStateProblem{…}, args::NonlinearSolvePolyAlgorithm{…}; sensealg::Nothing, u0::Nothing, p::Nothing, wrap::Val{…}, kwargs::@Kwargs{})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HGITF/src/solve.jl:1043
 [31] solve(prob::SteadyStateProblem{…}, args::NonlinearSolvePolyAlgorithm{…})
    @ DiffEqBase ~/.julia/packages/DiffEqBase/HGITF/src/solve.jl:1033
 [32] loss(x::Vector{…}, opt_ps::Tuple{…})
    @ Main ~/dev/fww_ss.jl:32
 [33] FixTail
    @ ~/.julia/packages/DifferentiationInterface/F5K7v/src/utils/context.jl:7 [inlined]
 [34] vector_mode_dual_eval!
    @ ~/.julia/packages/ForwardDiff/UBbGT/src/apiutils.jl:24 [inlined]
 [35] vector_mode_gradient!(result::Vector{…}, f::DifferentiationInterface.FixTail{…}, x::Vector{…}, cfg::ForwardDiff.GradientConfig{…})
    @ ForwardDiff ~/.julia/packages/ForwardDiff/UBbGT/src/gradient.jl:98
 [36] gradient!
    @ ~/.julia/packages/ForwardDiff/UBbGT/src/gradient.jl:39 [inlined]
 [37] gradient!(f::typeof(loss), grad::Vector{…}, prep::DifferentiationInterfaceForwardDiffExt.ForwardDiffGradientPrep{…}, backend::AutoForwardDiff{…}, x::Vector{…}, contexts::Constant{…})
    @ DifferentiationInterfaceForwardDiffExt ~/.julia/packages/DifferentiationInterface/F5K7v/ext/DifferentiationInterfaceForwardDiffExt/onearg.jl:362
 [38] (::OptimizationBase.var"#grad#16"{…})(res::Vector{…}, θ::Vector{…})
    @ OptimizationBase ~/.julia/packages/OptimizationBase/gvXsf/src/OptimizationDIExt.jl:28
 [39] (::LBFGSB.L_BFGS_B)(func::Optimization.var"#16#28"{…}, grad!::OptimizationBase.var"#grad#16"{…}, x0::Vector{…}, bounds::Matrix{…}; m::Int64, factr::Float64, pgtol::Float64, iprint::Int64, maxfun::Int64, maxiter::Int64)
    @ LBFGSB ~/.julia/packages/LBFGSB/UZibA/src/wrapper.jl:62
 [40] __solve(cache::OptimizationCache{…})
    @ Optimization ~/.julia/packages/Optimization/qX4vR/src/lbfgsb.jl:240
 [41] solve!(cache::OptimizationCache{…})
    @ SciMLBase ~/.julia/packages/SciMLBase/sYmAV/src/solve.jl:187
 [42] solve(::OptimizationProblem{…}, ::Optimization.LBFGS; kwargs::@Kwargs{})
    @ SciMLBase ~/.julia/packages/SciMLBase/sYmAV/src/solve.jl:95
 [43] solve(::OptimizationProblem{…}, ::Optimization.LBFGS)
    @ SciMLBase ~/.julia/packages/SciMLBase/sYmAV/src/solve.jl:92
 [44] top-level scope
    @ ~/dev/fww_ss.jl:42
Some type information was truncated. Use `show(err)` to see complete types.

Environment (please complete the following information):

• Output of using Pkg; Pkg.status()

  [961ee093] ModelingToolkit v9.64.3
  [8913a72c] NonlinearSolve v4.4.0
  [0bca4576] SciMLBase v2.75.1
  [2efcf032] SymbolicIndexingInterface v0.3.38

• Output of using Pkg; Pkg.status(; mode = PKGMODE_MANIFEST)

• Output of versioninfo()

Julia Version 1.11.3
Commit d63adeda50d (2025-01-21 19:42 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: Linux (x86_64-linux-gnu)
  CPU: 32 × Intel(R) Core(TM) i9-14900K
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, alderlake)
Threads: 32 default, 0 interactive, 16 GC (on 32 virtual cores)
Environment:
  JULIA_EDITOR = code

Additional context

Add any other context about the problem here.

[AayushSabharwal]

SteadyStateProblem doesn't have a custom remake in SciMLBase, so it doesn't promote the initialization problem.