Auto-generated bounds forwarded to init NonlinearProblem inject -Inf into cos -> NaN (regression from #4308)

[baggepinnen]

Description

Since SciML/ModelingToolkit.jl#4308 ("Forward unknown bounds metadata to NonlinearProblem/NonlinearLeastSquaresProblem"), NonlinearProblem construction attaches lb/ub vectors built from the unknowns' bounds metadata. The problem: generate_nonlinear_bounds emits a full lb/ub vector as soon as any unknown has a finite bound, filling in -Inf/Inf for every unbounded unknown.

# lib/ModelingToolkitBase/src/problems/nonlinearproblem.jl
function generate_nonlinear_bounds(sys::AbstractSystem, op)
    dvs = unknowns(sys)
    isempty(dvs) && return nothing, nothing
    lb = first.(getbounds.(dvs))   # unbounded unknown -> -Inf
    ub = last.(getbounds.(dvs))    # unbounded unknown ->  Inf
    ...
    if all(==(-Inf), lb) && all(==(Inf), ub)
        return nothing, nothing
    end
    return lb, ub                  # one finite bound anywhere => keep the -Inf/Inf entries
end

These spurious ±Inf bounds are passed straight into the (initialization) NonlinearProblem. They cause two problems:

1. NaN from the residual/Jacobian. Any solver step that evaluates at a bound substitutes -Inf into the model. For an unbounded angle appearing inside cos/sin, cos(-Inf) == NaN, so the residual becomes NaN and initialization fails.
2. The ±Inf bounds also confuse the nonlinear solve itself — several solvers stall, hit MaxIters, or clamp the unbounded component onto the spurious bound instead of converging.

The user never requested a bounded solve here — the bounds are auto-generated, and an unbounded unknown should not contribute a -Inf lower bound to a bounds-aware solver.

Minimal reproduction (pendulum)

A simple pendulum (x = L*sin(θ), y = -L*cos(θ)) with an initialization that specifies the tip location. θ is unbounded and appears inside cos; an auxiliary state b carries a finite physical bound. Two mutually-coupled initialization equations keep θ and b in the same SCC, so the init becomes a single NonlinearProblem holding both unknowns.

using ModelingToolkit
using ModelingToolkit: t_nounits as t, D_nounits as D
using NonlinearSolve

@parameters L = 1.0 g = 9.81
# θ: pendulum angle — unbounded (default bounds (-Inf, Inf)), appears inside cos/sin.
# b: an auxiliary state carrying a finite (physical) bound.
@variables θ(t) ω(t) x(t) y(t) b(t) [bounds = (0.0, 10.0)]

eqs = [
    D(θ) ~ ω
    D(ω) ~ -(g / L) * sin(θ)
    D(b) ~ 0
    x ~ L * sin(θ)
    y ~ -L * cos(θ)
]

# Two mutually-coupled initialization equations so that θ and b end up in the
# same strongly-connected component => one NonlinearProblem with both unknowns.
init_eqs = [
    cos(b) + sin(θ) ~ 1.0
    sin(b) + x ~ 0.5          # x = L*sin(θ): the tip location enters here
]

@named pend = System(eqs, t; initialization_eqs = init_eqs,
    guesses = [θ => 0.3, b => 0.5, ω => 0.0])
pend = mtkcompile(pend)

prob = ODEProblem(pend, [ω => 0.0], (0.0, 1.0))
initprob = prob.f.initialization_data.initializeprob

println("init unknowns : ", unknowns(initprob.f.sys))
println("initprob.lb   : ", initprob.lb)   # => [-Inf, 0.0]  (θ is unbounded)
println("initprob.ub   : ", initprob.ub)   # => [Inf, 10.0]

# The unbounded unknown's lower bound is -Inf. Any solver step that evaluates the
# residual/Jacobian at the bound substitutes -Inf into cos -> NaN.
res = similar(initprob.u0)
initprob.f(res, initprob.lb, initprob.p)
println("residual at lb: ", res)                          # => [NaN, NaN]
println("has NaN/Inf   : ", any(v -> isnan(v) || isinf(v), res))

Output:

init unknowns : SymbolicUtils...[θ(t), b(t)]
initprob.lb   : [-Inf, 0.0]
initprob.ub   : [Inf, 10.0]
residual at lb: [NaN, NaN]
has NaN/Inf   : true

Even smaller, solver-level reproduction

Stripping away the ODE, the same defect shows up in a plain coupled NonlinearProblem, where it also breaks the solve:

using ModelingToolkit, NonlinearSolve

@variables a [bounds = (0.0, 10.0)] θ      # a bounded, θ unbounded and inside cos
eqs = [
    0 ~ cos(a) + sin(θ) - 0.5
    0 ~ sin(a) + cos(θ) - 0.3
]
@named ns = System(eqs); ns = mtkcompile(ns)

prob = NonlinearProblem(ns, [a => 1.0, θ => 0.3])
@show prob.lb           # [-Inf, 0.0]
@show prob.ub           # [Inf, 10.0]

res = similar(prob.u0)
prob.f(res, prob.lb, prob.p)
@show res               # [NaN, NaN]   (cos(-Inf))

solve(prob)             # retcode = Stalled (clamps θ onto the spurious bound)

Suggested fix

generate_nonlinear_bounds should not hand -Inf/Inf entries to the solver for unbounded unknowns. Options:

• Only build/pass bounds when all relevant entries are finite, or otherwise leave unbounded entries as the solver's "no bound" sentinel rather than literal ±Inf.
• Or skip forwarding bounds entirely unless the chosen solver actually requires/consumes them.

Version info

• ModelingToolkit v11.28.3 (also present on current master; introduced by Add bounds generation to NonlinearProblems #4308)
• NonlinearSolve v4.20.0
• OrdinaryDiffEqCore v4.4.0
• Julia v1.12.6