Is there anyone from the JuliaNLSolvers organization here? I see the organization has only two members: JuliaNLSolvers · GitHub
Perhaps a good time to look for additional members?
Someone should approve this workflow: handle underflow of norm(J[:,i])^2 by stevengj · Pull Request #296 · JuliaNLSolvers/NLsolve.jl · GitHub
For reference, NonlinearSolve.jl handles this effectively:
opened 03:58AM - 19 Feb 26 UTC
closed 04:24AM - 19 Feb 26 UTC
## Reference
NLsolve.jl PR: https://github.com/JuliaNLSolvers/NLsolve.jl/pull/2… 96
Discourse thread: https://discourse.julialang.org/t/bug-in-nlsolve/135660
## The NLsolve.jl bug
NLsolve.jl's trust region dogleg implementation computes Jacobian column norms for autoscaling:
```julia
cache.d[j] = norm(view(jacobian(df), :, j))
if cache.d[j] == zero(cache.d[j])
cache.d[j] = one(cache.d[j])
end
```
Later it divides by `d.^2`. When column norms are tiny but nonzero (e.g. `~1e-200`), the norm passes the `== 0` check, but `d[j]^2` underflows to `0.0`, causing division-by-zero and NaN propagation.
## NonlinearSolve.jl is not affected
Tested with the reproducing example from the Discourse thread (Jacobian with `1e-200` column entries) — NonlinearSolve.jl produces no NaN, no Inf, no crash. The solvers return `Stalled` or `MaxIters` gracefully (the test system is intentionally inconsistent).
### Why it's robust
**TrustRegion (Dogleg):** The dogleg implementation (`NonlinearSolveBase/src/descent/dogleg.jl`) does not use Jacobian-column-norm autoscaling at all. It computes norms of step directions (`l_grad = internalnorm(δu_cauchy)`) rather than Jacobian columns, so the `norm(J[:,j])^2 → underflow` pattern never occurs.
**LevenbergMarquardt:** The LM implementation (`NonlinearSolveFirstOrder/src/levenberg_marquardt.jl`) computes `sum(abs2, J; dims=1)` (squared column norms) for the diagonal damping matrix `DᵀD`. However, `DᵀD` is initialized with a `min_damping_D` floor (default `1e-8`) and `update_levenberg_marquardt_diagonal!!` uses `max(y.diag[i], x[i])`, so diagonal entries can never go below `1e-8`. This prevents the underflow-to-zero problem.
## Test code
```julia
using NonlinearSolve
function f!(F, x, p)
F[1] = x[1] - 1.0
F[2] = x[1]^2 - 4.0
end
function j!(J, x, p)
J[1, 1] = 1.0
J[2, 1] = 2.0 * x[1]
J[1, 2] = 1e-200
J[2, 2] = 1e-200
end
x0 = [0.5, 0.5]
prob = NonlinearProblem(NonlinearFunction(f!; jac=j!), x0, nothing)
sol_tr = solve(prob, TrustRegion()) # Stalled, no NaN
sol_lm = solve(prob, LevenbergMarquardt()) # MaxIters, no NaN
```
This issue is for documentation/reference purposes only — no code changes needed.
Just to be clear: AI and I spent three hours migrating to NonlinearSolve. It doesn’t work. It cannot solve the problem that NLsolve handles easily. Therefore, for me, maintaining and fixing NLsolve is essential.
This code works: KiteModels.jl/test/test-steady-state-kps4.jl at bestie · OpenSourceAWE/KiteModels.jl · GitHub
The same code with NonlinearSolve doesn’t work. I do not yet have an MWE as proof.
I ran that code with NonlinearSolve TrustRegion and… it just worked out of the box, it even passed the one test that NLsolve.jl didn’t pass.
opened 07:23AM - 19 Feb 26 UTC
## Context
[KiteModels.jl](https://github.com/OpenSourceAWE/KiteModels.jl) uses… `NLsolve.jl` to find the steady-state equilibrium of a 4-point kite power system (3D mass-spring network with aerodynamic forces). We'd like to migrate to `NonlinearSolve.jl` but the current code relies on several NLsolve-specific features.
## Current NLsolve usage
The core call pattern is:
```julia
using NLsolve
# Manual central finite-difference Jacobian with subnormal flushing
function make_jac(f!, n_vars)
_F1 = zeros(Float64, n_vars)
_F2 = zeros(Float64, n_vars)
_xp = zeros(Float64, n_vars)
_xm = zeros(Float64, n_vars)
threshold = sqrt(floatmin(Float64))
function jac!(J, x)
h_factor = cbrt(eps(Float64))
for j in 1:n_vars
copyto!(_xp, x); copyto!(_xm, x)
h = max(abs(x[j]), one(Float64)) * h_factor
_xp[j] += h; _xm[j] -= h
f!(_F1, _xp); f!(_F2, _xm)
@views J[:, j] .= (_F1 .- _F2) ./ (2h)
end
@. J = ifelse(abs(J) < threshold, zero(Float64), J)
end
return jac!
end
n_unknowns = 20 # 2*(segments + kite_particles - 1) + 2
X00 = zeros(Float64, n_unknowns)
jac! = make_jac(f!, n_unknowns)
results = nlsolve(f!, jac!, X00;
autoscale=true, xtol=4e-7, ftol=4e-7, iterations=200)
if NLsolve.converged(results)
solution = results.zero
end
```
### NLsolve features we depend on
1. **Manual Jacobian (`jac!`)**: Needed because `NLSolversBase >= 7.10` produces subnormal Jacobian entries via DifferentiationInterface, which cause NaN in NLsolve's autoscale (subnormal column norms squared underflow to zero). The manual Jacobian flushes near-zero entries.
2. **`autoscale=true`**: Critical for convergence — the system variables have very different scales (positions in meters vs. small perturbations).
3. **Separate `xtol` and `ftol`**: Both set to `4e-7`.
4. **`iterations` limit**: Set to 200.
5. **`results.zero`** and **`converged(results)`**: Solution extraction and convergence check.
### Questions
1. What is the equivalent `NonlinearSolve.jl` API for this call pattern?
2. Does NonlinearSolve handle autoscaling? (I see [#425](https://github.com/SciML/NonlinearSolve.jl/issues/425) is open for trust region autoscaling.)
3. Does NonlinearSolve's automatic differentiation avoid the subnormal Jacobian issue, or would we still need a manual Jacobian?
4. What's the recommended solver algorithm for this kind of problem (20-variable nonlinear system from physics, starting from a reasonable initial guess)?
## Self-contained MWE
Below is a **complete, self-contained** MWE (~960 lines) that reproduces the KiteModels.jl steady-state test without any kite-specific package dependencies. It models a 4-point kite on a segmented tether with aerodynamic forces, spring/damping, and a winch. The nonlinear solve finds particle positions where all accelerations vanish (steady-state equilibrium).
**Dependencies**: `StaticArrays`, `LinearAlgebra`, `NLsolve`, `Dierckx`, `Test`
<details>
<summary>Click to expand full MWE (mwe_steady_state_kps4.jl)</summary>
```julia
$(cat /home/crackauc/sandbox/tmp_20260219_005256_56738/mwe_steady_state_kps4.jl)
```
</details>
### Running the MWE
```julia
] add StaticArrays NLsolve Dierckx
include("mwe_steady_state_kps4.jl")
```
Expected output (all 15 tests pass):
```
l=100.0: tether_length=100.92196416459542, pre_tension=1.0003226874576083 ✓
l=200.0: tether_length=202.59541320501268, pre_tension=1.0004541973108771 ✓
l=392.0: tether_length=399.15270436705833, pre_tension=1.0006386343184874 ✓
Test Summary: | Pass Total Time
test_find_steady_state (MWE) | 15 15 15.3s
```
## Environment
- Julia 1.11
- NLsolve v4.5.1
- StaticArrays v1.9.16
So here’s a PR which passes all of your tests locally:
bestie ← ChrisRackauckas-Claude:replace-nlsolve-with-nonlinearsolve
opened 09:58AM - 19 Feb 26 UTC
## Summary
- Replaces `NLsolve.jl` trust region dogleg with `NonlinearSolve.jl`… `TrustRegion()` in both `KPS3` and `KPS4` `find_steady_state!` functions
- Removes `NLsolve` and `NLSolversBase` as direct dependencies (NonlinearSolve was already a transitive dependency via OrdinaryDiffEqNonlinearSolve)
- Retains the custom finite-difference Jacobian (`make_jac`) since the model's `MVector{3, Float64}` internals cannot hold `ForwardDiff.Dual` numbers
## Motivation
NLsolve.jl has known issues with subnormal Jacobian entries from NLSolversBase >= 7.10 (see PR #259). Rather than pinning NLSolversBase or maintaining workarounds, this PR switches to NonlinearSolve.jl which:
1. **Already works** — `TrustRegion()` passes all 452 tests (441 pass + 11 pre-existing broken)
2. **Produces better results** — On the l=392 steady-state test case, NonlinearSolve TrustRegion achieves 0.92% tether length error vs NLsolve's 1.35% (which fails the 1% rtol check)
3. **Is faster** — ~0.03s vs ~0.15s per solve (post-JIT)
4. **Is already a dependency** — No new packages added to the dependency tree
5. **Doesn't need autoscale** — NLsolve's `autoscale=true` (Moré/MINPACK diagonal scaling) is unnecessary with NonlinearSolve's TrustRegion
## Changes
| File | Change |
|------|--------|
| `src/KPS4.jl` | Replace `nlsolve()` call with `NonlinearProblem` + `solve(prob, TrustRegion())` |
| `src/KPS3.jl` | Same replacement |
| `src/KiteModels.jl` | Remove `NLsolve` from `using`, add `import SciMLBase` |
| `src/init.jl` | Update `make_jac` docstring |
| `Project.toml` | Remove NLsolve/NLSolversBase deps, add SciMLBase |
| `test/Project.toml` | Remove NLsolve dep |
| `examples_3d/Project.toml` | Remove NLsolve dep |
| `test/create_sys_image.jl` | Remove NLsolve from sysimage |
## Test plan
- [x] All 15 KPS4 steady-state tests pass
- [x] All 67 KPS3 tests pass
- [x] Full `Pkg.test()` passes: 441 pass + 11 pre-existing broken, 0 failures
## Related
- NonlinearSolve.jl issue: https://github.com/SciML/NonlinearSolve.jl/issues/836
- KiteModels.jl PR #259 (NLSolversBase subnormal workaround)
🤖 Generated with [Claude Code](https://claude.com/claude-code)
Is there some test there that isn’t in the repo?
1 Like
It doesn’t pass the tests locally:
julia> include("test/test-steady-state-kps4.jl")
Activating project at `~/repos/KiteModels.jl/test`
┌ Warning: find_steady_state!: solver did not converge! retcode=Stalled
└ @ KiteModels ~/repos/KiteModels.jl/src/KPS4.jl:704
test_find_steady_state: Test Failed at /home/ufechner/repos/KiteModels.jl/test/test-steady-state-kps4.jl:42
Expression: pre_tension > 1.0001
Evaluated: 1.0 > 1.0001
Stacktrace:
[1] macro expansion
@ ~/.julia/juliaup/julia-1.11.9+0.x64.linux.gnu/share/julia/stdlib/v1.11/Test/src/Test.jl:680 [inlined]
[2] macro expansion
@ ~/repos/KiteModels.jl/test/test-steady-state-kps4.jl:42 [inlined]
[3] macro expansion
@ ~/.julia/juliaup/julia-1.11.9+0.x64.linux.gnu/share/julia/stdlib/v1.11/Test/src/Test.jl:1709 [inlined]
[4] top-level scope
@ ~/repos/KiteModels.jl/test/test-steady-state-kps4.jl:30
┌ Warning: find_steady_state!: solver did not converge! retcode=MaxIters
└ @ KiteModels ~/repos/KiteModels.jl/src/KPS4.jl:704
test_find_steady_state: Test Failed at /home/ufechner/repos/KiteModels.jl/test/test-steady-state-kps4.jl:42
Expression: pre_tension > 1.0001
Evaluated: 1.0 > 1.0001
Stacktrace:
[1] macro expansion
@ ~/.julia/juliaup/julia-1.11.9+0.x64.linux.gnu/share/julia/stdlib/v1.11/Test/src/Test.jl:680 [inlined]
[2] macro expansion
@ ~/repos/KiteModels.jl/test/test-steady-state-kps4.jl:42 [inlined]
[3] macro expansion
@ ~/.julia/juliaup/julia-1.11.9+0.x64.linux.gnu/share/julia/stdlib/v1.11/Test/src/Test.jl:1709 [inlined]
[4] top-level scope
@ ~/repos/KiteModels.jl/test/test-steady-state-kps4.jl:30
test_find_steady_state: Test Failed at /home/ufechner/repos/KiteModels.jl/test/test-steady-state-kps4.jl:45
Expression: isapprox(tether_length(kps4_local), 1.008954l, rtol = 0.01)
Evaluated: isapprox(197.42057783848486, 201.7908; rtol = 0.01)
Stacktrace:
[1] macro expansion
@ ~/.julia/juliaup/julia-1.11.9+0.x64.linux.gnu/share/julia/stdlib/v1.11/Test/src/Test.jl:680 [inlined]
[2] macro expansion
@ ~/repos/KiteModels.jl/test/test-steady-state-kps4.jl:45 [inlined]
[3] macro expansion
@ ~/.julia/juliaup/julia-1.11.9+0.x64.linux.gnu/share/julia/stdlib/v1.11/Test/src/Test.jl:1709 [inlined]
[4] top-level scope
@ ~/repos/KiteModels.jl/test/test-steady-state-kps4.jl:30
┌ Warning: find_steady_state!: solver did not converge! retcode=MaxIters
└ @ KiteModels ~/repos/KiteModels.jl/src/KPS4.jl:704
test_find_steady_state: Test Failed at /home/ufechner/repos/KiteModels.jl/test/test-steady-state-kps4.jl:42
Expression: pre_tension > 1.0001
Evaluated: 1.0 > 1.0001
Stacktrace:
[1] macro expansion
@ ~/.julia/juliaup/julia-1.11.9+0.x64.linux.gnu/share/julia/stdlib/v1.11/Test/src/Test.jl:680 [inlined]
[2] macro expansion
@ ~/repos/KiteModels.jl/test/test-steady-state-kps4.jl:42 [inlined]
[3] macro expansion
@ ~/.julia/juliaup/julia-1.11.9+0.x64.linux.gnu/share/julia/stdlib/v1.11/Test/src/Test.jl:1709 [inlined]
[4] top-level scope
@ ~/repos/KiteModels.jl/test/test-steady-state-kps4.jl:30
test_find_steady_state: Test Failed at /home/ufechner/repos/KiteModels.jl/test/test-steady-state-kps4.jl:45
Expression: isapprox(tether_length(kps4_local), 1.008954l, rtol = 0.01)
Evaluated: isapprox(389.9800737647081, 395.50996799999996; rtol = 0.01)
Stacktrace:
[1] macro expansion
@ ~/.julia/juliaup/julia-1.11.9+0.x64.linux.gnu/share/julia/stdlib/v1.11/Test/src/Test.jl:680 [inlined]
[2] macro expansion
@ ~/repos/KiteModels.jl/test/test-steady-state-kps4.jl:45 [inlined]
[3] macro expansion
@ ~/.julia/juliaup/julia-1.11.9+0.x64.linux.gnu/share/julia/stdlib/v1.11/Test/src/Test.jl:1709 [inlined]
[4] top-level scope
@ ~/repos/KiteModels.jl/test/test-steady-state-kps4.jl:30
Test Summary: | Pass Fail Total Time
test_find_steady_state | 10 5 15 1.6s
ERROR: LoadError: Some tests did not pass: 10 passed, 5 failed, 0 errored, 0 broken.
in expression starting at /home/ufechner/repos/KiteModels.jl/test/test-steady-state-kps4.jl:29
For whatever reason, the only test that was running on your pull request on Github was the test if the licenses are correct: Replace NLsolve.jl with NonlinearSolve.jl TrustRegion · OpenSourceAWE/KiteModels.jl@eb90218 · GitHub
It’s a difference of Julia v1.12 vs v1.11. Interesting.
Found it and fixed. That was an interesting computer and Julia-version specific bug
10 Likes
pleiby
February 24, 2026, 7:30pm
11
If I may, you continue to amaze me Chris. Thanks!
but it should all be good now.