Notebooks · lib/power/pf

AC and DC power flow by Newton-Raphson, fast-decoupled and DC

A power flow finds the bus voltages V = Vm e^(j Va) at which the power each bus injects into the network, S = V conj(Y V), equals what its generators produce minus what its loads draw. Each bus fixes two of its four quantities (P, Q, Vm, Va) by its type:

PQ (a load bus): P and Q are given, Vm and Va are unknown;
PV (a generator bus): P and Vm are given, Va is unknown;
REF (the slack bus): Vm and Va are given; its P and Q make up the losses and are known only after the solve.

runpf(c, method, enforce_q, tol, max_it) solves a case by one of three methods, the three of MATPOWER's runpf (R. D. Zimmerman, C. E. Murillo-Sanchez and R. J. Thomas, "MATPOWER: steady-state operations, planning and analysis tools for power systems research and education", IEEE Trans. Power Systems 26 (2011) 12-19; the MATPOWER manual, chapter 4):

NR: Newton's method on the real and imaginary mismatches in polar coordinates, with the Jacobian dS/dVa, dS/dVm in closed form, factored by the sparse LU of power.splu each iteration;
FDXB: the fast-decoupled method of Stott and Alsac (IEEE Trans. PAS 93 (1974) 859-869) in the XB form of van Amerongen (IEEE Trans. Power Systems 4 (1989) 760-770): two constant real matrices B' and B'', factored once, half-steps in Va and in Vm;
DC: the linear model with Vm = 1, no resistance and sin(dVa) = dVa, one sparse solve for the angles.

The result reports the largest mismatch left, per unit (for the fast-decoupled method, divided by Vm as that method measures it). A flow that does not converge in max_it iterations reports its last iterate with success false.

With enforce_q a generator whose reactive output leaves [Qmin, Qmax] is held at the limit, its bus becomes a PQ bus, and the flow is solved again until no limit is violated (MATPOWER's ENFORCE_Q_LIMS = 1).

The result has MATPOWER's columns filled: bus Vm and Va (degrees), generator Pg and Qg (MW, MVAr) and the four branch flow columns (power entering each end, MW and MVAr).

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

This program declares a record type (Net), and records run only in compiled C for now, so the Run buttons here are off. Open it in the workbench and use its C tab to compile it.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 45–870

from power.case import (
    BR_STATUS, BR_X, BRANCH_COLS, BUS_I, BUS_TYPE, Case, F_BUS, GEN_BUS, GEN_STATUS, GS, NONE, PD, PF, PG, PQ, PT,
    PV, QD, QF, QG, QMAX, QMIN, QT, REF, SHIFT, T_BUS, TAP, VA, VG, VM,
)
from power.ybus import DEG, Ybus, YbusOpts, full_network, make_ybus, ybus_mul
from power.splu import SpLU, lu_solve, mindeg_order, splu

#: methods
NR: Int = 1
FDXB: Int = 2
DC: Int = 3

#: the diagonal is kept as pivot when within this factor of the column's largest entry
PIVOT_TOL: Num = 0.001
#: MATLAB's eps, in MATPOWER's splitting of Q between generators at one bus
MACHINE_EPS: Num = 2.220446049250313e-16
#: a generator's Q may exceed its limit by this much, in MVAr (MATPOWER's opf.violation)
QLIM_TOL: Num = 5e-6

record Net {
    nb: Int
    base: Num
    bus: Arr
    gen: Arr
    branch: Arr
    gbus: [Int]
    fbus: [Int]
    tbus: [Int]
}

record PFResult {
    bus: Arr
    gen: Arr
    branch: Arr
    success: Bool
    iterations: Int
    method: Int
    mismatch: Num
}

record Newton {
    v: CArr
    converged: Bool
    iterations: Int
    mismatch: Num
}

#| The case with every bus, generator end and branch end numbered 0..nb-1
#| in bus order, and the branch matrix widened to hold the flow columns.
def internal(c: Case) -> Net {
    nb: Int = int(shape(c.bus)[0])
    ng: Int = int(shape(c.gen)[0])
    nl: Int = int(shape(c.branch)[0])
    idx = map()
    for i to nb {
        idx[c.bus[i, BUS_I]] = i
    }
    gbus: [Int] = [int(idx[c.gen[g, GEN_BUS]]) for g to ng]
    fbus = int_zeros(nl)
    tbus = int_zeros(nl)

from power.case import (
    BR_STATUS, BR_X, BRANCH_COLS, BUS_I, BUS_TYPE, Case, F_BUS, GEN_BUS, GEN_STATUS, GS, NONE, PD, PF, PG, PQ, PT,
    PV, QD, QF, QG, QMAX, QMIN, QT, REF, SHIFT, T_BUS, TAP, VA, VG, VM,
)
from power.ybus import DEG, Ybus, YbusOpts, full_network, make_ybus, ybus_mul
from power.splu import SpLU, lu_solve, mindeg_order, splu

#: methods
NR: Int = 1
FDXB: Int = 2
DC: Int = 3

#: the diagonal is kept as pivot when within this factor of the column's largest entry
PIVOT_TOL: Num = 0.001
#: MATLAB's eps, in MATPOWER's splitting of Q between generators at one bus
MACHINE_EPS: Num = 2.220446049250313e-16
#: a generator's Q may exceed its limit by this much, in MVAr (MATPOWER's opf.violation)
QLIM_TOL: Num = 5e-6

record Net {
    nb: Int
    base: Num
    bus: Arr
    gen: Arr
    branch: Arr
    gbus: [Int]
    fbus: [Int]
    tbus: [Int]
}

record PFResult {
    bus: Arr
    gen: Arr
    branch: Arr
    success: Bool
    iterations: Int
    method: Int
    mismatch: Num
}

record Newton {
    v: CArr
    converged: Bool
    iterations: Int
    mismatch: Num
}

#| The case with every bus, generator end and branch end numbered 0..nb-1
#| in bus order, and the branch matrix widened to hold the flow columns.
def internal(c: Case) -> Net {
    nb: Int = int(shape(c.bus)[0])
    ng: Int = int(shape(c.gen)[0])
    nl: Int = int(shape(c.branch)[0])
    idx = map()
    for i to nb {
        idx[c.bus[i, BUS_I]] = i
    }
    gbus: [Int] = [int(idx[c.gen[g, GEN_BUS]]) for g to ng]
    fbus = int_zeros(nl)
    tbus = int_zeros(nl)