Replacement of `NLsolve.jl`

[langestefan]

Hi!

PowerModels.jl currently depends on NLsolve.jl to compute the solution of the non-linear system with compute_ac_pf. NLsolve hasn't had a release since December 2020 and is currently not actively maintained. The current release v4.5.1 is limited to NLsolversBase v7, with the change to v8 requiring a breaking release of NLsolve.jl: JuliaNLSolvers/NLsolve.jl#295 (comment)

I spoke to Patrick on slack and he adviced that it is best to move off from NLsolve since nobody is really looking at it anymore, not for years.

Then the question is, what to replace it with? I did some experiments with NonlinearSolve.jl, that seems to be the most natural replacement. Alternatively we could roll our own implementation of Newton Rhapson and just drop the dependency entirely.

Please let me know your opinion and I can make a PR.

[odow]

I did some experiments with NonlinearSolve.jl, that seems to be the most natural replacement

I'm a little nervous of NonlinearSolve.jl because of the many dependencies. But it would be the most standard option that's probably the easiest.

Another downside is that it is a breaking change because we currently pass kwargs to NLsolve.nlsolve, and I assume these don't transfer to NonlinearSolve.jl.

Alternatively we could roll our own implementation of Newton Rhapson and just drop the dependency entirely

This is also an option. We could just vendor NLsolve.jl, removing the stuff we don't use? It seems to be working okay.

[langestefan]

I'm a little nervous of NonlinearSolve.jl because of the many dependencies. But it would be the most standard option that's probably the easiest.

That's fair, it is a huge dependency. It would be nice to support as an opt in, through an extension to SciMLBase. That would give access to the entire suit of solvers in SciML. The cache reuse functionality already gives quite nice speedups (~2x) for repeated solves on 1000 bus networks compared to NLsolve.jl (see also: https://sciml.ai/news/2024/01/23/nlsolve/).

This is also an option. We could just vendor NLsolve.jl, removing the stuff we don't use? It seems to be working okay.

The implementation in NLsolve.jl seems to rely on a number of other dependencies: https://github.com/JuliaNLSolvers/NLsolve.jl/blob/84a3e60d00068458310cb62f24f1ef9b55b8e5ef/src/solvers/newton.jl#L2

it seems to be using functions from OnceDifferentiable.jl, LineSearches.jl and Distances.jl. And we don't need a generic implementation, we have an analytical jacobian already. So writing our own version will be faster and cleaner I think.

[langestefan]

The following could be a nice compromise:

• We remove NLsolve as a dependency. This will be a breaking change, because we also stop passing kwargs to NLsolve.nlsolve. I don't see a way to implement this without making it a breaking change.
• We update the solve interface to compute_ac_pf(data; solver = NativeNewton(), ...).
• We ship an implementation of NativeNewton that is at least as good as the one in NLsolve, but can probably be better. Does not have to be more than 100 lines because it will be PowerModels.jl specific.
• We add two new structs, PowerFlowSystem and PowerFlowSolution to standardize input/output data.
• PowerFlowSystem is an object that holds the system construction stuff that is now inside _compute_ac_pf. Such as f!, j!, x0, p0, and the sparsity pattern which some solvers like to have. It allows us to easily cache things such as the factorization. PowerFlowSolution holds a solver agnostic result that any solver backend can define its own post processing for. THe post processing that is now part of compute_ac_pf can consume this object.
• The actual dispatch to a solver is handled by _solve_nl(sys::PowerFlowSystem, alg). By default this would dispatch to our implementation NativeNewton.
• An extension can now easily add new solvers to PowerModels by dispatching on _solve_nl.
• For SciML, we would add an extension PowerModelsSciMLExt that would dispatch on _solve_nl(sys::PowerFlowSystem, alg::SciMLBase.AbstractNonlinearAlgorithm)

Effectively this removes a dependency, adds no new dependencies, and makes the solver interface extensible.

The effort can be split across two PRs:

1. Deprecation of NLsolve and implementation of NativeNewton
2. Adding the extension for PowerModelsSciMLExt

[jd-lara]

@langestefan we already have serveral implementations of NR in this package https://github.com/Sienna-Platform/PowerFlows.jl/ native to PowerSystems infrastructure

[langestefan]

Great, they can be made available through the same extension path

[odow]

It looks like there was a new release of NLsolve.jl, so the need to do something here seems less urgent.

[langestefan]

Yes, I'm happy Patrick responded. However, he did indicate we should move away from NLsolve. As most of the ecosystem has already moved on NLsolve will not be maintained in the future.

I did a benchmark with NonlinearSolve.jl to show that moving away from NLsolve to an extension architecture has the potential to bring significant speedups: langestefan#1 (comment), anywhere from 4x in single solves to 10x in repeated solves with caching.

Maybe we can consider making one big release that just rips off the bandaid:

• Drop finite_differencing (see: CompatHelper: bump compat for NLsolve to 5, (keep existing compat) #1015 (comment))
• Consider this PR
• Adopt reduced Jacobian formulation (see: Change compute_ac_pf iteration to use more standard (smaller) Jacobian #1011 (comment))