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

[ChrisRackauckas-Claude]

I reproduced this on current master (453c45c08e, ModelingToolkitBase v1.61.0, NonlinearSolve v4.26.0, Julia 1.12.6). The results split three ways: one MWE doesn't reproduce, one does but for a different reason than the report gives, and there's a third genuine NaN bug in generate_nonlinear_bounds that the report's title describes accurately but whose stated mechanism is not the actual one.

1. The pendulum MWE does not reproduce — the init problem carries no bounds

Running the pendulum script verbatim on master:

init unknowns : [θp(t), var"b#0"(t)]
initprob.lb   : nothing
initprob.ub   : nothing
init retcode  : Success  u = [0.08867835208848147, 0.42403103949074045]

initprob.lb === nothing, so nothing is forwarded to the initialization NonlinearProblem and initialization succeeds. Note b appears in the init system as the renamed var"b#0"(t), which does not carry the bounds metadata. So the "auto-generated bounds forwarded to init NonlinearProblem" framing in the title doesn't hold for this system on current master.

2. The second MWE does reproduce a real failure, but -Inf is not the cause

The solve failure is real:

prob.lb  : [-Inf, 0.0]        # u = [θ, a]
prob.ub  : [Inf, 10.0]
retcode  : Stalled            u = [-0.955, 1.75e-39]   |resid| = 0.316

and stripping the bounds fixes it (retcode = Success). But the diagnosis in the issue — that -Inf gets substituted into cos — isn't what happens. NonlinearSolveBase/src/bounds_transform.jl branches per component on isfinite:

_from_unbounded(t, lb, ub) = has_lb && has_ub ? lb + (ub - lb) * logistic(t) :
                             has_lb ? lb + exp(t) : has_ub ? ub - exp(t) : t

An unbounded component maps through the identity, so the solver never evaluates the residual at -Inf. The prob.f(res, prob.lb, prob.p) call in the report evaluates the residual at -Inf by hand; no solver step does that, so [NaN, NaN] there isn't evidence of the failure mechanism.

Swapping the infinities for finite numbers shows what's actually responsible:

unbounded                  retcode=Success   u=[-6.5278, -0.7345]  |resid|=2.3e-15
a ∈ [0,10]                 retcode=Stalled   u=[-0.9550,  1.7e-39] |resid|=0.316
a ∈ [0,10], θ ∈ [-100,100] retcode=Stalled   u=[-0.9287,  0.0111]  |resid|=0.310   # all finite, still stalls
θ ∈ [-100,100] only        retcode=Success   u=[-6.5278, -0.7345]  |resid|=1.5e-14   # a unbounded, fine

With all-finite bounds it still stalls; with only the unbounded variable left free it succeeds. The culprit is the finite bound on a and the log/logit reparameterization, not ±Inf. And a root does exist strictly inside the box (a = 5.5487, θ = -0.2446, residual ~4e-16), so the solver is failing on a problem it should solve.

This reproduces with no ModelingToolkit in the picture at all:

using NonlinearSolve
f!(r, u, p) = (r[1] = cos(u[2]) + sin(u[1]) - 0.5; r[2] = sin(u[2]) + cos(u[1]) - 0.3; r)
u0 = [0.3, 1.0]

solve(NonlinearProblem(NonlinearFunction(f!), u0, nothing))                              # Success
solve(NonlinearProblem(NonlinearFunction(f!), u0, nothing;
                       lb = [-100.0, 0.0], ub = [100.0, 10.0]))                          # Stalled

So that half is a NonlinearSolve.jl bounds-transform robustness issue, not an MTK one. I've filed it separately: SciML/NonlinearSolve.jl#1147.

The MTK-side question it raises is a design one rather than a bug: auto-forwarding bounds metadata turns a solve that previously succeeded into one that stalls. #4308 asked for bounds "to help the solver process", and here they hurt. Whether NonlinearProblem should opt users into a reparameterized constrained solve purely on the presence of bounds metadata — which is often there for optimization/control purposes rather than as a constraint on a root-find — is worth a separate decision. I haven't changed that behavior.

3. A real NaN bug in generate_nonlinear_bounds (symbolic-bounds path)

Separately, while checking the above I found a genuine defect in generate_nonlinear_bounds that produces exactly the NaN in this issue's title, via a different route. It fires when any bound is symbolic, so the numeric-bounds MWEs above never hit it.

using ModelingToolkit, NonlinearSolve

@parameters lo = -1.0
@variables c [bounds = (lo, 1.0)] d          # `lo` is symbolic; `d` is unbounded
@mtkcompile psys = System([0 ~ c^2 + d^2 - 1, 0 ~ c^2 - d^2 - 0.5])
prob = NonlinearProblem(psys, [c => 0.8, d => 0.4])

prob.lb   # [-1.0, 0.4]   <- 0.4 is d's *initial guess*, not a bound
prob.ub   # [ 1.0, 0.4]
solve(prob)  # retcode = Unstable, u = [0.8, NaN]

The cause is at lib/ModelingToolkitBase/src/problems/nonlinearproblem.jl:95-109. The map is populated only for finite bounds:

for (var, b) in zip(dvs, lb)
    b === -Inf && continue                      # unbounded unknowns get no entry
    write_possibly_indexed_array!(lbmap, var, Symbolics.SConst(b), COMMON_NOTHING)
end
left_merge!(lbmap, op)                          # ...so `op` supplies them — as u0 values
lb = varmap_to_vars(lbmap, dvs; tofloat = false)

Merging the operating point is what lets a symbolic bound like lo resolve, but it also hands every unbounded unknown its operating-point value. That pins the unknown to lb == ub == u0, and the bounds transform then computes logit((u - lb) / (ub - lb)) = logit(0/0) = NaN. The all-numeric branch is unaffected and correctly yields [-1.0, -Inf] / [1.0, Inf].

PR: #4902 — restores the ±Inf sentinel after resolution, with a test that fails before the fix (retcode = Unstable, u = [0.8, NaN]) and passes after.

Summary

• Pendulum/init MWE: not reproducible on master; init problem gets lb = ub = nothing.
• Coupled NonlinearProblem MWE: real failure, but caused by NonlinearSolve's handling of the finite bound, not by -Inf. Filed upstream.
• The suggested fix ("only build bounds when all entries are finite") would not address either, and would regress the deliberate partial-bounds case asserted at lib/ModelingToolkitBase/test/nonlinearsystem.jl:564.
• A separate, genuine NaN bug does exist in generate_nonlinear_bounds on the symbolic-bounds path; PR linked above.