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, GenCost, Load, Branch
from gpf.scenario import Study
from gpf import parse_raw, parse_rawx, write_raw, write_rawx

The demo case#

A small two-area system serves every recipe. In practice a case is parsed from a file (load_case, which the file-name suffix dispatches); 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. Affine cost curves accompany the generators for the optimal-power-flow recipe; the power-flow recipes ignore them.

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, pb=0, pt=300),
            Generator(ibus=2, machid="1", pg=60.0, vs=1.02, ireg=3, qb=-200, qt=200, pb=0, pt=200),
            Generator(ibus=5, machid="1", pg=60.0, vs=1.01, qb=-200, qt=200, pb=0, pt=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
        ],
        # Affine \$/MWh curves, read only by the OPF layer: bus 1 is cheapest.
        gen_costs=[
            GenCost(ibus=1, machid="1", model=2, cost=[20.0, 0.0]),
            GenCost(ibus=2, machid="1", model=2, cost=[30.0, 0.0]),
            GenCost(ibus=5, machid="1", model=2, cost=[40.0, 0.0]),
        ],
    )

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 inspect the swing setup before solving#

assert gpf.validate_case(demo_case()).can_solve    # ready to solve
odd = demo_case()
odd.buses[3].ide = 3                                # mark bus 4 as a swing bus...
# ...but bus 4 has no machine, so its slack injection answers to none
assert 4 in gpf.analyze_case(odd).orphan_slack_buses

Read per-bus and per-generator quantities#

case = demo_case()
sol = gpf.solve_power_flow(case)
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())
trace = [(r.iteration, r.max_mismatch) for r in sol.step_history]
assert sol.converged                               # the solve reached tolerance
assert trace[0][1] > trace[-1][1]                  # mismatch shrank across the trace

Modify and re-solve#

Editing records by identity is introduced in Working with Cases.

Scale load in one area, then re-solve#

case = demo_case()
base = gpf.solve_power_flow(case)
gpf.scale_loads(case, areas=[2], multiplier=1.3)   # +30% load in area 2, in place
after = gpf.solve_power_flow(case)
# 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)
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)
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)
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").sweep().run()
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

Get PTDF / LODF factors at the solved point#

with gpf.sensitivity(demo_case()) as snap:
    ptdf = snap.ptdf(monitored=[(1, 5, "1")])       # tie flow vs every bus
    lodf = snap.lodf(monitored=[(3, 6, "1")], outage=[(1, 5, "1")])
shift = lodf.values[0, 0]         # share of the tripped tie's flow picked up
assert 0.9 < shift < 1.1          # the only other tie takes (nearly) all of it
assert ptdf.as_dict()["1-5-1"][1] == 0.0            # slack-bus column is zero

Judge N-1 against criteria and read the report#

A branch N-1 monitors every finalized AC branch by default. This example makes the selected slot explicit and refuses missing ratings. The report exposes the verdict and coverage separately from the ranked violations.

case = demo_case()
for branch in case.branches:
    branch.rate1 = 200.0                        # exact selected rating slot
study = (Study(case, name="summer N-1")
         .monitor(voltage_band=(0.95, 1.05),
                  base_rating_set=1, contingency_rating_set=1,
                  missing_rating="error")
         .sweep())
result = study.run()
doc = json.loads(result.to_json())              # one versioned, violations-first report
criteria = doc["meta"]["criteria"]              # the explicit pass/fail criteria
worst = doc["summary"]["worst_contingencies"]   # ranked by severity
assert doc["schema_version"] == "1.3"
assert criteria["voltage_band"] == [0.95, 1.05]
assert criteria["base_rating_set"] == criteria["contingency_rating_set"] == 1
assert doc["summary"]["coverage"]["rated_contingency_branches"] == 10
assert doc["summary"]["verdict"] in {"fail", "incomplete"}
assert worst                                    # the radial-spur trip islands bus 8

Compare two solutions bus by bus#

base = gpf.solve_power_flow(demo_case())
stressed = demo_case()
gpf.scale_loads(stressed, areas=[2], multiplier=1.3)
alt = gpf.solve_power_flow(stressed)
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())
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())
for ibus, v in golden.items():
    assert abs(again.get_bus_voltage(int(ibus))[0] - v) < 1e-6

Dispatch and prices#

Price the least-cost dispatch (DC-OPF)#

sol = gpf.solve_dc_opf(demo_case())           # least-cost DC dispatch, priced
price = {b: round(sol.lmp[b], 1) for b in sol.lmp}          # nodal $/MWh
split = {b: round(sol.lmp_congestion[b], 1) for b in sol.lmp}  # congestion part
# the cheapest unit (bus 1) carries the load; with no binding branch limit the
# price is uniform at that unit's marginal cost and the congestion part is zero
assert round(sol.pg[(1, "1")], 1) == 195.0
assert set(price.values()) == {20.0} and set(split.values()) == {0.0}

For the full arc (security-constrained dispatch, corrective recourse, and multi-period pricing) see the optimal power flow tutorial.

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)
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)
assert sol.any_limit_binding                       # the bus-2 machine cannot hold 1.05

Files and studies#

Save a case to RAW and read it back#

case = demo_case()
base = gpf.solve_power_flow(case)
path = pathlib.Path(tempfile.mkdtemp()) / "case.raw"
write_raw(case, path)                                   # write RAW (v33)
reloaded = parse_raw(str(path))              # read it back
again = gpf.solve_power_flow(reloaded)
assert round(again.get_bus_voltage(1)[0], 6) == round(base.get_bus_voltage(1)[0], 6)

Convert a case between formats#

case = demo_case()
base = gpf.solve_power_flow(case)
path = pathlib.Path(tempfile.mkdtemp()) / "case.rawx"
write_rawx(case, path)                                  # write a different format
roundtrip = gpf.solve_power_flow(parse_rawx(path))
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())  # .con carries no base case
status = result.status_counts()                      # how the contingencies resolved
assert status                                        # the imported study ran

Version history#

Walked end to end in the Case management tutorial.

Track edits to a case file, then roll back#

Tracking lives in a .gpflog sidecar next to the case; the case file stays a standard RAW/RAWX. Opening the ledger reconciles the file on disk, so edits made by hand (or any other tool) become versions too. See also gpf history and gpf checkout in the gpf Reference.

path = pathlib.Path(tempfile.mkdtemp()) / "case.rawx"
write_rawx(demo_case(), path)
gpf.open_ledger(path)                       # start tracking: version 0
case = parse_rawx(path)
case.loads[0].pl += 20.0                    # hand-edit the file and re-save
write_rawx(case, path)
led = gpf.open_ledger(path)                 # the edit is captured as version 1
assert [row["seq"] for row in led.history()] == [0, 1]
original = led.checkout(0)                  # any version, from the log alone
assert original.loads[0].pl == case.loads[0].pl - 20.0

Hand a project to a colleague#

The bundle contains the project and its declared resources. Import verifies its content before restoration and refuses conflicting destination files. The CLI equivalents are gpf export and gpf import.

path = pathlib.Path(tempfile.mkdtemp()) / "case.rawx"
write_rawx(demo_case(), path)
led = gpf.open_ledger(path)                 # tracked: version 0
bundle = path.parent / "case_bundle.zip"
project = gpf.write_project(path.parent / "project.toml", path.name,
                            name="Grid", history=path.name + ".gpflog")
gpf.export_bundle(project, bundle)
inbox = pathlib.Path(tempfile.mkdtemp())    # a colleague's machine
received = gpf.import_bundle(bundle, inbox)  # tamper-verified on import
assert gpf.read_ledger(received.path).head_hash == led.head_hash