Add new algorithm, the semi-implicit Continuous Newton method

[rzyu45]

What kind of problems is it mostly used for? Please describe.

Nonlinear algebraic equations, including but not limited to the AC power flow models.

Describe the algorithm you’d like

As depicted in [1], for equations

$$0=g(y),$$

we can convert it to

$$\begin{aligned} \dot{y}&=z\\\ 0&=J(y)z+g(y) \end{aligned},$$

where $J$ is the jaobian, and then solve the DAEs with the new 4-stage 3rd-order L-stable RODAS3D constructed in [1]. The values of $z$ can be initialized by $-J^{-1}(y_0)g(y_0)$ for consistency.

It is proved in [1] that if the initial values are within the region of attraction of the converted dynamic systems, then the continuous Newton method is bounded to converge. So that the semi-implicit Continuous Newton method is very stable.

Other implementations to know about

Currently there is no other implementation to the best of my knowledge. Maybe I can do it here. It would be nice to have some developer notes.

References
[1] Semi-implicit Continuous Newton Method for Power Flow Analysis
[2] Implicit Continuous Newton Method for Power Flow Analysis
[3] Continuous Newton's Method for Power Flow Analysis

[avik-pal]

Yeah these would be interesting to have. Though we might have to implement this in something like SteadyStateDiffEq, since we can't take a dep on OrdinaryDiffEq (circular dep).

Let us know if you need help getting started. I am guessing the DAE solver goes in OrdinaryDiffEq (@ChrisRackauckas)?

[ChrisRackauckas]

Yeah the DAE solver could go into OrdinaryDiffEq and we can have this in SteadyStateDiffEq.jl, I was thinking the same thing. It is a special DAE solver so it's not one that we have right now, but I think it's the right choice to focus it entirely on stability and not accuracy given accurately solving the ODE isn't the goal but making it stable for big dts is.

[rzyu45]

@avik-pal @ChrisRackauckas
I have understood the dep relationships. Could you please offer me some dev tutorials? The way to use the existing interfaces to generate Jac and Hessian-Vector Products will be of significance.

[avik-pal]

For the Jacobian See NewtonRaphson implementation https://github.com/SciML/NonlinearSolve.jl/blob/aa282b71019eb7a80b91464242f7d80599a5b0da/src/raphson.jl#L76-L122C4. Essentially you have jacobian_caches create the cache for you and you call cache.J = jacobian!!(cache.J, cache) to construct the jacobian. This handles all possible cases of Krylov Methods, Sparse Jacobian, etc...

None of the algorithms currently use Hessian Vector Products. But we have them implemented here