Recipes#

Each recipe answers one task with the minimum code and the value that confirms it. Every recipe assumes the imports and the demo_case() builder below.

import json
import pathlib
import tempfile

import gpf
from gpf import GpfCase, Bus, Generator, Load, Branch
from gpf.scenario import Study

The demo case#

A small two-area system serves every recipe. In practice a case is parsed from a file (GpfCase.from_raw, GpfCase.from_rawx, or GpfCase.from_matpower); here it is built in memory so each recipe runs on its own. Areas 1 and 2 are joined by tie lines, a radial spur reaches bus 8, and the generator on bus 2 regulates remote bus 3.

def demo_case() -> GpfCase:
    # Two areas (1 and 2) joined by tie lines, with a radial spur to bus 8.
    # The bus-2 machine regulates remote bus 3 to its setpoint (1.02 pu).
    return GpfCase(
        sbase=100.0,
        buses=[
            Bus(ibus=1, ide=3, vm=1.02, area=1, baskv=230.0),  # slack
            Bus(ibus=2, ide=2, vm=1.02, area=1, baskv=230.0),  # PV, regulates bus 3
            Bus(ibus=3, ide=1, vm=1.0, area=1, baskv=230.0),
            Bus(ibus=4, ide=1, vm=1.0, area=1, baskv=230.0),
            Bus(ibus=5, ide=2, vm=1.01, area=2, baskv=230.0),  # PV
            Bus(ibus=6, ide=1, vm=1.0, area=2, baskv=230.0),
            Bus(ibus=7, ide=1, vm=1.0, area=2, baskv=230.0),
            Bus(ibus=8, ide=1, vm=1.0, area=2, baskv=230.0),  # radial spur
        ],
        generators=[
            Generator(ibus=1, machid="1", pg=0.0, vs=1.02, qb=-300, qt=300),
            Generator(ibus=2, machid="1", pg=60.0, vs=1.02, ireg=3, qb=-200, qt=200),
            Generator(ibus=5, machid="1", pg=60.0, vs=1.01, qb=-200, qt=200),
        ],
        loads=[
            Load(ibus=3, pl=50, ql=15, area=1), Load(ibus=4, pl=40, ql=12, area=1),
            Load(ibus=6, pl=45, ql=14, area=2), Load(ibus=7, pl=35, ql=10, area=2),
            Load(ibus=8, pl=25, ql=8, area=2),
        ],
        branches=[
            Branch(ibus=1, jbus=2, r=0.01, x=0.08),
            Branch(ibus=2, jbus=3, r=0.01, x=0.08),
            Branch(ibus=3, jbus=4, r=0.01, x=0.08),
            Branch(ibus=4, jbus=1, r=0.01, x=0.08),
            Branch(ibus=5, jbus=6, r=0.01, x=0.08),
            Branch(ibus=6, jbus=7, r=0.01, x=0.08),
            Branch(ibus=7, jbus=5, r=0.01, x=0.08),
            Branch(ibus=1, jbus=5, r=0.01, x=0.10),  # tie
            Branch(ibus=3, jbus=6, r=0.01, x=0.10),  # tie
            Branch(ibus=4, jbus=8, r=0.01, x=0.08),  # radial to bus 8
        ],
    )

Modify and re-solve#

Scale load in one area, then re-solve#

case = demo_case()
base = gpf.solve_power_flow(case, backend="cpu")
gpf.scale_loads(case, areas=[2], multiplier=1.3)   # +30% load in area 2, in place
after = gpf.solve_power_flow(case, backend="cpu")
# the extra demand sags voltage in the scaled area
assert after.get_bus_voltage(6)[0] < base.get_bus_voltage(6)[0]

Edit one device by identity, then re-solve#

case = demo_case()
gen = next(g for g in case.generators if g.ibus == 5 and g.machid == "1")
gen.vs = 1.04                                       # raise the bus-5 setpoint
sol = gpf.solve_power_flow(case, backend="cpu")
assert round(sol.get_bus_voltage(5)[0], 2) == 1.04  # bus 5 holds the new setpoint

Trip a branch, then re-solve#

case = demo_case()
base = gpf.solve_power_flow(case, backend="cpu")
tie = next(br for br in case.branches if (br.ibus, br.jbus) == (1, 5))
tie.stat = 0                                        # open the area-tie line
after = gpf.solve_power_flow(case, backend="cpu")
assert after.converged                              # tie loss is survivable here
# area-2 load buses now feed through the one remaining tie, so voltages shift
assert after.get_bus_voltage(7)[0] != base.get_bus_voltage(7)[0]

Screen and compare#

Rank single-branch contingencies#

result = Study(demo_case(), name="N-1").generate_n1().run(backend="cpu")
ranked = result.top_voltage_drops(3)        # worst first: [(scenario, dV pu, bus), ...]
failed = [s.label for s in result.failed()]  # contingencies that could not solve
# tripping the radial spur islands bus 8, so its contingency lands in `failed`
assert ranked and failed

Compare two solutions bus by bus#

base = gpf.solve_power_flow(demo_case(), backend="cpu")
stressed = demo_case()
gpf.scale_loads(stressed, areas=[2], multiplier=1.3)
alt = gpf.solve_power_flow(stressed, backend="cpu")
dv = {b: round(alt.get_bus_voltage(b)[0] - base.get_bus_voltage(b)[0], 4)
      for b in range(1, 9)}
worst_bus = min(dv, key=dv.get)                    # bus with the biggest drop
assert worst_bus in (6, 7, 8)                      # one of the stressed area-2 buses

Snapshot a solution and guard against drift#

sol = gpf.solve_power_flow(demo_case(), backend="cpu")
snapshot = {b.ibus: round(sol.get_bus_voltage(b.ibus)[0], 6) for b in demo_case().buses}
path = pathlib.Path(tempfile.mkdtemp()) / "voltages.json"
path.write_text(json.dumps(snapshot))

# later, re-solve and confirm nothing has drifted
golden = json.loads(path.read_text())
again = gpf.solve_power_flow(demo_case(), backend="cpu")
for ibus, v in golden.items():
    assert abs(again.get_bus_voltage(int(ibus))[0] - v) < 1e-6

Controls and limits#

Hold a remote bus and check it held#

case = demo_case()                                 # bus-2 machine regulates bus 3
sol = gpf.solve_power_flow(case, backend="cpu")
assert round(sol.get_bus_voltage(3)[0], 3) == 1.02  # the remote target is held

Find a generator pegged at a reactive limit#

# A PV machine with a tight Q range cannot hold its voltage under heavy
# reactive demand, so it pegs at its limit.
case = GpfCase(
    sbase=100.0,
    buses=[
        Bus(ibus=1, ide=3, vm=1.0, area=1, baskv=230.0),
        Bus(ibus=2, ide=2, vm=1.05, area=1, baskv=230.0),
        Bus(ibus=3, ide=1, vm=1.0, area=1, baskv=230.0),
    ],
    generators=[
        Generator(ibus=1, machid="1", pg=0.0, vs=1.0, qb=-500, qt=500),
        Generator(ibus=2, machid="1", pg=50.0, vs=1.05, qb=-10, qt=10),  # tight Q
    ],
    loads=[Load(ibus=3, pl=150, ql=120)],          # heavy reactive demand
    branches=[
        Branch(ibus=1, jbus=2, r=0.01, x=0.08),
        Branch(ibus=2, jbus=3, r=0.01, x=0.08),
        Branch(ibus=1, jbus=3, r=0.01, x=0.10),
    ],
)
sol = gpf.solve_power_flow(case, backend="cpu")
assert sol.any_limit_binding                       # the bus-2 machine cannot hold 1.05

Inspect a case#

Inventory what a case contains#

counts = gpf.analyze_case(demo_case()).counts
assert (counts.n_buses, counts.n_generators) == (8, 3)
assert counts.n_generators_remote == 1             # one remote regulator
assert counts.n_isolated_buses == 0

Validate, and catch a setup error before solving#

assert gpf.validate_case(demo_case()).can_solve    # ready to solve
broken = demo_case()
broken.buses[3].ide = 3                             # mark bus 4 as a swing bus...
# ...but bus 4 has no generator to act as the slack source
assert 4 in gpf.analyze_case(broken).orphan_slack_buses

Read per-bus and per-generator quantities#

case = demo_case()
sol = gpf.solve_power_flow(case, backend="cpu")
voltages = {b.ibus: round(sol.get_bus_voltage(b.ibus)[0], 4) for b in case.buses}
gen_q_mvar = {g.ibus: round(sol.get_generator_q(g.ibus) * case.sbase, 1)
              for g in case.generators}
assert voltages[1] == 1.02                         # slack held at its setpoint
assert gen_q_mvar[2] != 0.0                        # the bus-2 machine carries reactive power

Read the convergence trace#

sol = gpf.solve_power_flow(demo_case(), backend="cpu")
trace = [(r.outer_iter, r.inner_iter, r.max_mismatch) for r in sol.iteration_history]
assert sol.converged                               # the solve reached tolerance
assert trace[0][2] > trace[-1][2]                  # mismatch shrank across the trace

Files and studies#

Convert a case between formats#

case = demo_case()
base = gpf.solve_power_flow(case, backend="cpu")
path = pathlib.Path(tempfile.mkdtemp()) / "case.rawx"
case.to_rawx(path)                                  # write a different format
roundtrip = gpf.solve_power_flow(GpfCase.from_rawx(path), backend="cpu")
assert round(roundtrip.get_bus_voltage(1)[0], 6) == round(base.get_bus_voltage(1)[0], 6)

Import a .con / .sub / .mon contingency study#

d = pathlib.Path(tempfile.mkdtemp())
(d / "study.con").write_text(
    "CONTINGENCY 'tie-1-5'\n"
    "DISCONNECT BRANCH FROM BUS 1 TO BUS 5 CIRCUIT '1'\n"
    "END\n"
    "END\n"
)
(d / "study.sub").write_text("SUBSYSTEM 'AREA1'\nAREA 1\nEND\n")
(d / "study.mon").write_text("MONITOR BRANCHES IN SUBSYSTEM 'AREA1'\nEND\n")
study = Study.from_psse(d / "study.con", sub=d / "study.sub", mon=d / "study.mon")
result = study.run(case=demo_case(), backend="cpu")  # .con carries no base case
status = result.status_counts()                      # how the contingencies resolved
assert status                                        # the imported study ran