Skip to content

Type promotion error when differentiating SteadyStateProblem #3430

@SebastianM-C

Description

@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.

Metadata

Metadata

Labels

bugSomething isn't working

Type

No type

Projects

No projects

Milestone

No milestone

Relationships

None yet

Development

No branches or pull requests

Issue actions