Power flow convergence and timing from 9 to 9241 buses

Newton-Raphson against the fast-decoupled method on MATPOWER cases from 9 to 2383 buses: the mismatch left after each iteration, the iterations each needs, and the time per solve measured against an independent SciPy power flow, up to the 9241-bus PEGASE model of the European grid.

The mismatch, iteration by iteration

max_it=k stops the solver after k iterations, and mismatch is the largest power imbalance left at any bus, in per unit of the 100 MVA base. Newton's method roughly squares the error each step; the fast-decoupled method, whose matrices stay fixed, gains a constant factor per step. Each line ends at the iteration where the mismatch first falls below the tolerance, 1e-8.

import json

import numpy as np
import vipyr as vp
from vipyr.figure import figure

cases = ["case30", "case118", "case300", "case2383wp"]
hist = {}
for name in cases:
    for method, top in (("nr", 6), ("fdxb", 18)):
        ms = []
        for k in range(top + 1):
            r = vp.power.runpf(name, method=method, max_it=k)
            ms.append(r.mismatch)
            if r.success:
                break
        hist[name, method] = ms
for name in cases:
    print("%-10s nr  " % name, " ".join("%.0e" % m for m in hist[name, "nr"]))
print("case2383wp fdxb", " ".join("%.0e" % m for m in hist["case2383wp", "fdxb"]))

fig = figure("Largest mismatch after k iterations", "iteration k", "mismatch (per unit)")
fig.yscale("log")
for name in ("case118", "case2383wp"):
    for method in ("nr", "fdxb"):
        m = np.array(hist[name, method])
        fig.plot(np.arange(len(m)), m, "%s %s" % (name, method))
fig.hline("tolerance", 1e-8)

case30     nr   4e-01 2e-02 7e-05 1e-09
case118    nr   1e+00 2e-02 1e-05 1e-12
case300    nr   9e+00 3e+00 4e-01 8e-03 5e-06 1e-12
case2383wp nr   1e+03 1e+02 6e+00 6e-01 8e-03 1e-06 6e-12
case2383wp fdxb 1e+03 2e+02 6e-01 1e-01 2e-02 5e-03 1e-03 3e-04 1e-04 4e-05 2e-05 6e-06 2e-06 8e-07 3e-07 1e-07 4e-08 1e-08 5e-09

Iterations and time per solve

The table comes from lib/power/compare.py, which runs every case by both libraries and keeps the best of several runs; it was measured on an Apple M2 Pro under load from other jobs, so the times are indicative. Viv is the compiled solver called through vipyr; the SciPy reference builds the same sparse matrices with scipy.sparse and solves with SuperLU. case9241pegase is not shipped with vipyr (it is 1.5 MB); the comparison downloads it from the MATPOWER repository.

bench = json.load(open("lib/power/compare.json"))
rows = [r for r in bench["cases"]]
print("%-15s %6s  %-4s %5s  %10s  %10s" % ("case", "buses", "", "iter", "Viv (ms)", "SciPy (ms)"))
for r in rows:
    print("%-15s %6d  %-4s %5d  %10.2f  %10s" % (r["case"], r["nb"], r["method"], r["iterations"], 1e3 * r["t_viv"],
                                                 "%.2f" % (1e3 * r["t_ref"]) if r["method"] != "fdxb" else ""))

fig = figure("Time per power flow", "buses", "seconds")
fig.xscale("log")
fig.yscale("log")
for method, label in (("nr", "Viv NR"), ("fdxb", "Viv FDXB"), ("dc", "Viv DC")):
    sel = [r for r in rows if r["method"] == method]
    fig.plot([r["nb"] for r in sel], [r["t_viv"] for r in sel], label)
sel = [r for r in rows if r["method"] == "nr"]
fig.plot([r["nb"] for r in sel], [r["t_ref"] for r in sel], "SciPy NR")

case             buses        iter    Viv (ms)  SciPy (ms)
case9                9  nr       4        0.15        5.29
case9                9  fdxb     6        0.04            
case9                9  dc       1        0.03        0.46
case14              14  nr       2        0.05        3.18
case14              14  fdxb     6        0.05            
case14              14  dc       1        0.03        0.46
case30              30  nr       3        0.11        4.43
case30              30  fdxb    11        0.08            
case30              30  dc       1        0.04        0.48
case118            118  nr       3        0.31        4.99
case118            118  fdxb     8        0.21            
case118            118  dc       1        0.10        0.61
case300            300  nr       5        1.80       10.95
case300            300  fdxb     9        0.57            
case300            300  dc       1        0.24        1.04
case2383wp        2383  nr       6       22.71       57.02
case2383wp        2383  fdxb    18        6.77            
case2383wp        2383  dc       1        2.45        4.10
case9241pegase    9241  nr       6      134.82      202.22
case9241pegase    9241  fdxb    14       44.98            
case9241pegase    9241  dc       1       15.73       17.10

Credits

The methods follow MATPOWER's (BSD-3-Clause). The test cases are MATPOWER's case files, each with its source named in its header.

Source

gallery/power/plot_pf_convergence.py, run as written when these pages were built. The output and the figures above are from that run.

The Viv program it calls is lib/power/api.viv: its atlas, the defs and the calls between them, or open it in the workbench.