Comparing variants of NR

[luke-kiernan]

I did a comparison of the 3 variants of Newton Raphson: vanilla Newton Raphson (NR), Trust Region (TR), and Levenberg-Marquardt (LM). On the matpower_ACTIVSg2000_sys system, NR (top) and TR (bottom) are comparable in runtime. Both take 4 iterations.

On the larger matpower_ACTIVSg10k_sys system, TR (bottom) converges while NR (top) doesn't. TR is a little over 2x faster. Here, TR converges in 14 iterations, while NR does 30 iterations (the default maxIterations), so per-iteration runtime is fairly comparable. I'm not sure why NR has much higher memory usage: that might be worth looking into.

Here's LM on the 2k bus system, where it takes 10 iterations to converge. It's significantly slower, even on a per-iteration basis. It fails to converge for the larger 10k bus system. (I didn't include the screenshot for that one.) I'm convinced something's amiss with my implementation or choice of parameters: at this time, my code is still configured for debugging, not for performance.

I'm planning to do a side-by-side comparison of LM and TR this week: if the two choose wildly different search directions, then hopefully that'll point me in the right direction as to what's amiss (weighting, choice of damping parameter, etc.).

Here's the script I'm using for these profiling comparisons: I ran NR and TR on the main branch, and LM on the levenberg-marquardt branch.

using PowerSystemCaseBuilder
using PowerFlows
const PSB = PowerSystemCaseBuilder
const PF = PowerFlows
using BenchmarkTools
systems = Dict(
    "ACTIVSg 2k buses" => PSB.build_system(MatpowerTestSystems, "matpower_ACTIVSg2000_sys"),
    "ACTIVSg 10k buses" => PSB.build_system(MatpowerTestSystems, "matpower_ACTIVSg10k_sys"),
)
pf_methods = Dict(
    "Newton-Raphon" => NewtonRaphsonACPowerFlow,
    "Trust Region" => PF.TrustRegionACPowerFlow,
    "L-M" => PF.LevenbergMaquardtACPowerFlow,
)
for (pf_name, pf_type) in pf_methods
    for (sys_name, sys) in systems
        pf_inst = ACPowerFlow{pf_type}()
        data = PowerFlowData(pf_inst, sys)
        display(@benchmark solve_powerflow!(data_copy, pf = $pf_inst) setup =
            (data_copy = deepcopy($data)))
        println("$pf_name on $sys_name")
    end
end

[luke-kiernan]

In a similar vein, I made flame graphs for each method. To open them, run using ProfileView; ProfileView.view(nothing), then click the folder icon and navigate to the file. GitHub won't let me upload them with the .jlprof file extension, so I zipped them.

newton-raphson-flame-graph.jlprof.zip
trust-region-flame-graph.jlprof.zip

Keep in mind that each takes a different number of iterations: 14 for TR, 30 for NR. Also, I'm using the same x0 for all trials.

Newton-Raphson: largest chunks are _update_jacobian_values! at about 50% and klu_refactor at 25-30%. _update_residual_values! is more like 10%.

Trust region: the two largest chunks are _update_jacobian_values! and klu_refactor, both at around 30%-40% of runtime. Things with impact around more like 5%: _update_residual_values!, solving for Δx, klu_factor, and klu_analyze.

Script I used to generate the flame graphs:

using Profile, ProfileView
using PowerSystemCaseBuilder, PowerFlows
const PF = PowerFlows
const PSB = PowerSystemCaseBuilder
sys = PSB.build_system(MatpowerTestSystems, "matpower_ACTIVSg10k_sys")
const PF_TYPE = NewtonRaphsonACPowerFlow # change as desired.
pf = PF.ACPowerFlow{PF_TYPE}()
data = PF.PowerFlowData(pf, sys)
p_inj = deepcopy(data.bus_activepower_injection)
p_with = deepcopy(data.bus_activepower_withdrawals)
q_inj = deepcopy(data.bus_reactivepower_injection)
q_with = deepcopy(data.bus_reactivepower_withdrawals)
Va = deepcopy(data.bus_angles)
Vm = deepcopy(data.bus_magnitude)
# single in-place solve is ~0.2 seconds
function inplace_solve(data::PowerFlowData,
    powerflow::PF.ACPowerFlow,
    N::Int,
)
    for _ in 1:N
        # reset from last in-place powerflow before solving again
        copyto!(data.bus_activepower_injection, p_inj)
        copyto!(data.bus_activepower_withdrawals, p_with)
        copyto!(data.bus_reactivepower_injection, q_inj)
        copyto!(data.bus_reactivepower_withdrawals, q_with)
        copyto!(data.bus_angles, Va)
        copyto!(data.bus_magnitude, Vm)
        data.converged .= false
        PF.solve_powerflow!(data; pf = powerflow)
    end
end
@profview inplace_solve(data, pf, 1) # compile: ignore this one
@profview inplace_solve(data, pf, 20)