Examples¶
Thirteen runnable scripts live under examples/ in the repository. Each one is self-contained,
reads only files under fixtures/, prints a short deterministic summary and exits 0 in about a
second. They are executed on every push — by tests/unit/test_examples_run.py inside the test
matrix and by the dedicated examples CI job — and this page embeds them with
pymdownx.snippets, so the code you read here and the code CI ran are the same bytes.
Run any of them from a clone, from the repository root:
| Script | Shows | Manual page |
|---|---|---|
01_load_and_validate.py |
load_with_report, a hand-built network, all-issues validation |
Network model, File formats |
02_ac_power_flow.py |
solve_ac with and without Q-limits, voltages, loading table, warm start |
Power flow |
03_dc_power_flow.py |
solve_dc on case300, DC vs AC flows |
Power flow |
04_jobs_api.py |
SolveRequest / run / run_json, structured failures, captured warnings |
Jobs API |
05_roles_and_islands.py |
effective bus roles, NoSlackGeneratorError, island repair vs strict model |
Power flow, File formats |
06_network_matrices.py |
NetworkArrays, Ybus sparsity, PTDF, LODF with the bridge NaN column, a 3-bus case in full |
Numerics |
07_results_and_export.py |
JSON round trip, to_arrays(), CSV export |
Results |
08_opf_and_n1.py |
solve_dc_opf dispatch/duals/LMP, ac_check, congestion, contingency.n1 screen-then-confirm |
DC-OPF, N-1 screening |
09_nodal_market.py |
market.solve_nodal on a Scenario: elastic-demand dispatch, LMPs split by congestion, settlement |
Nodal market |
10_multiperiod_market.py |
market.solve_multiperiod over a 24-period Scenario: ramp coupling, storage SoC with efficiency, the cyclic horizon, per-period LMPs and settlement |
Multiperiod market |
11_zonal_redispatch.py |
market.solve_zonal: zonal clearing, min-cost redispatch, the nodal reference, corridor duals, the three gap figures and the settlement identity |
Zonal market |
12_agent_market.py |
market.solve_agents: generators offering through a Strategy, price-takers reproducing solve_nodal bit-exactly, a pivotal markup stopping at demand's own bid, the duopoly, termination_reason, StrategyConfig crossing jobs as JSON |
Agent-based bidding |
13_interop.py |
pandapower JSON export solved by pp.rundcpp, pandapower's case14 imported losslessly, PyPSA export optimised by PyPSA against solve_dc_opf, PSS/E RAW import with RAW_NO_COSTS, a bit-exact CSV bundle round trip, and one deliberately lossy conversion with its ExportReport |
File formats |
1. Load and validate¶
Load case14 with load_with_report and read the typed ImportIssue entries; build a 3-bus
network by hand; make the model reject a broken document with every issue listed at once.
"""Load a MATPOWER case, read the import report, and see validation fail with every issue.
What this shows:
* ``io.matpower.load_with_report`` returns the ``Network`` **and** an ``ImportReport`` whose
typed ``ImportIssue`` entries name every repair the importer made (here: case14 stores
``BASE_KV = 0`` on every bus, which becomes 1.0 kV with one issue per bus).
* A network can be built by hand from the entity classes; the model validates on construction.
* Validation reports **all** issues at once: one ``NetworkValidationError`` carrying a
``ValidationIssue`` (code, path, message) per problem, never just the first one.
* ``validate_network`` re-checks a mutated network and returns the issues instead of raising.
Run from the repository root: ``uv run python examples/01_load_and_validate.py``.
"""
from __future__ import annotations
from mambo_power.io import matpower
from mambo_power.model import (
Branch,
Bus,
Generator,
Load,
Network,
NetworkValidationError,
validate_network,
)
# --- 1. Load a case with its import report -------------------------------------------------
net, report = matpower.load_with_report("fixtures/matpower/case14.m")
print("case14:", len(net.buses), "buses,", len(net.branches), "branches,", end=" ")
print(len(net.generators), "generators,", len(net.loads), "loads,", len(net.shunts), "shunts")
print("base_mva:", net.base_mva, "| slack:", [b.id for b in net.buses if b.type == "slack"])
print("import report:", len(report.warnings), "issue(s), codes", sorted(report.codes))
first = report.warnings[0]
print("first issue:", first.code, "| buses", first.bus_ids, "|", first.message)
print("legacy string form:", report.as_strings()[0])
# --- 2. Build a tiny network by hand ------------------------------------------------------
mini = Network(
base_mva=100,
buses=[
Bus(id="b1", base_kv=110, type="slack"),
Bus(id="b2", base_kv=110, type="pv"),
Bus(id="b3", base_kv=110, type="pq"),
],
branches=[
Branch(id="l12", from_bus="b1", to_bus="b2", r=0.01, x=0.10, b=0.02),
Branch(id="l13", from_bus="b1", to_bus="b3", r=0.02, x=0.20, b=0.02),
Branch(id="l23", from_bus="b2", to_bus="b3", r=0.01, x=0.10, b=0.02),
],
generators=[
Generator(
id="g1",
bus="b1",
p_mw=0,
q_mvar=0,
p_min_mw=0,
p_max_mw=300,
q_min_mvar=-100,
q_max_mvar=100,
v_set_pu=1.02,
),
Generator(
id="g2",
bus="b2",
p_mw=60,
q_mvar=0,
p_min_mw=0,
p_max_mw=100,
q_min_mvar=-40,
q_max_mvar=40,
v_set_pu=1.01,
),
],
loads=[Load(id="d3", bus="b3", p_mw=120, q_mvar=40)],
)
print("hand-built network is valid:", validate_network(mini) == [])
# --- 3. Trigger a validation error and print every issue ----------------------------------
broken = {
"base_mva": 100,
"buses": [
{"id": "a", "base_kv": 110, "type": "pq"},
{"id": "a", "base_kv": 0, "type": "pq"},
],
"branches": [{"id": "x", "from_bus": "a", "to_bus": "zz", "r": 0.0, "x": 0.0, "b": 0.0}],
"generators": [
{
"id": "g",
"bus": "a",
"p_mw": 10,
"q_mvar": 0,
"p_min_mw": 50,
"p_max_mw": 20,
"q_min_mvar": 0,
"q_max_mvar": 0,
"v_set_pu": 1.0,
}
],
}
try:
Network.model_validate(broken)
except NetworkValidationError as err:
print(f"NetworkValidationError with {len(err.issues)} issue(s); codes {sorted(err.codes)}")
for issue in err.issues:
print(f" {issue.code:16s} {issue.path:22s} {issue.message}")
# --- 4. Mutation does not re-validate; validate_network does -----------------------------
mini.buses[2].base_kv = -1.0
issues = validate_network(mini)
print("after mutation:", [(i.code, i.path) for i in issues])
2. AC power flow¶
Newton-Raphson on case14 and case118 with Q-limits on and off: iterations, rounds, mismatch, the pinned generators, the first bus voltages, a branch-loading table and a warm start that needs zero iterations.
"""AC Newton-Raphson power flow on case14 and case118, with and without Q-limit enforcement.
What this shows:
* ``pf.solve_ac(net, options=AcOptions(...))`` — flat start, tolerance 1e-8 pu on the mismatch
infinity norm, sparse Jacobian factorised with ``scipy.sparse.linalg.splu``.
* ``AcPowerFlowResult`` diagnostics: ``iterations`` (summed over Q-limit rounds),
``q_limit_rounds``, ``max_mismatch_mva``, and ``GenResult.q_limited`` (``"min"``/``"max"``)
for the generators pinned at a reactive limit. On case118 with limits on, six buses are
pinned (the same set pandapower pins); with limits off nothing is pinned and bus 103 sits
at its 1.01 setpoint.
* Bus voltages and a branch-loading table from the typed result.
* Warm start: copy the solved state into ``Bus.vm_pu`` / ``Bus.va_deg`` and solve again with
``init="auto"`` — a start already inside tolerance reports 0 iterations.
Run from the repository root: ``uv run python examples/02_ac_power_flow.py``.
"""
from __future__ import annotations
from mambo_power import pf
from mambo_power.io import matpower
from mambo_power.results import AcPowerFlowResult
def summarise(label: str, result: AcPowerFlowResult) -> None:
pinned = [(g.id, g.bus, g.q_limited) for g in result.generators if g.q_limited != "none"]
print(f"--- {label}")
print(
f"converged={result.converged} iterations={result.iterations} "
f"q_limit_rounds={result.q_limit_rounds} max_mismatch={result.max_mismatch_mva:.2e} MVA"
)
print("pinned generators:", pinned if pinned else "none")
for name in ("case14", "case118"):
net = matpower.load(f"fixtures/matpower/{name}.m")
for q_limits in (True, False):
options = pf.AcOptions(init="flat", q_limits=q_limits)
summarise(f"{name}, q_limits={q_limits}", pf.solve_ac(net, options=options))
# --- A closer look at case118 with limits enforced ------------------------------------------
net = matpower.load("fixtures/matpower/case118.m")
result = pf.solve_ac(net, options=pf.AcOptions(init="flat"))
print("\nfirst 5 bus voltages (case118, q_limits on):")
for bus in result.buses[:5]:
print(f" {bus.id:8s} {bus.vm_pu:7.4f} pu {bus.va_deg:8.3f} deg role={bus.role_effective}")
print("\nbus 103 (limited gen) with and without Q-limits:")
off = pf.solve_ac(net, options=pf.AcOptions(init="flat", q_limits=False))
vm_on = next(b.vm_pu for b in result.buses if b.id == "bus-103")
vm_off = next(b.vm_pu for b in off.buses if b.id == "bus-103")
print(f" vm on={vm_on:.5f} pu vm off={vm_off:.5f} pu (= its 1.01 setpoint)")
# case118 ships no thermal ratings (RATE_A = 0 -> rating_mva None -> loading_pct None), so
# stamp a uniform 250 MVA rating to show the loading column; ratings do not affect the solve.
for branch in net.branches:
branch.rating_mva = 250.0
result = pf.solve_ac(net, options=pf.AcOptions(init="flat"))
print("\nfive most loaded branches (case118, q_limits on, 250 MVA on every branch):")
rated = [b for b in result.branches if b.loading_pct is not None]
for br in sorted(rated, key=lambda b: b.loading_pct or 0.0, reverse=True)[:5]:
flow = complex(br.p_from_mw, br.q_from_mvar)
print(
f" {br.id:10s} {br.from_bus:8s}->{br.to_bus:8s} "
f"P={br.p_from_mw:8.2f} MW Q={br.q_from_mvar:8.2f} MVAr "
f"|S|={abs(flow):7.2f} MVA loading={br.loading_pct:5.2f} %"
)
losses = sum(b.p_from_mw + b.p_to_mw for b in result.branches)
print(f"total active losses: {losses:.3f} MW")
# --- Warm start ------------------------------------------------------------------------------
# Copy the q_limits=False solution into the buses; "auto" then starts from it and the mismatch
# is already inside tolerance: 0 iterations. (With limits on, a pinned bus is PV again at the
# start, its magnitude snaps back to the setpoint, and one re-pin round is needed.)
state = {b.id: (b.vm_pu, b.va_deg) for b in off.buses}
for bus in net.buses:
bus.vm_pu, bus.va_deg = state[bus.id]
warm = pf.solve_ac(net, options=pf.AcOptions(init="auto", q_limits=False))
print(
f"\nwarm start from the solved state: iterations={warm.iterations} "
f"rounds={warm.q_limit_rounds} converged={warm.converged}"
)
print("provenance:", warm.provenance.kind, warm.provenance.solver, warm.provenance.options)
3. DC power flow¶
solve_dc on case300: angle range, the largest flows, the slack balance, and how the lossless
linear flows compare with the AC solution on the same case.
"""DC power flow on case300 and how far it sits from the AC solution.
What this shows:
* ``pf.solve_dc(net)`` — the lossless linear model ``B'θ = P − P_shift`` with the slack
angle fixed at 0, branch flows via ``B_f``; always converges on a connected network.
* The ``DcPowerFlowResult`` summary: angle range, the largest flows, the slack-generator
balance (the first in-service slack-bus generator absorbs it).
* A comparison with the AC Newton-Raphson flows on the same case: the DC approximation is
close on most branches and off by tens of MW on a few, and it knows nothing about the
AC losses.
Run from the repository root: ``uv run python examples/03_dc_power_flow.py``.
"""
from __future__ import annotations
import numpy as np
from mambo_power import pf
from mambo_power.io import matpower
net = matpower.load("fixtures/matpower/case300.m")
dc = pf.solve_dc(net)
print("case300 DC:", dc.provenance.kind, dc.provenance.solver, "converged =", dc.converged)
angles = np.array([b.va_deg for b in dc.buses])
print(f"angles: min {angles.min():.2f} deg, max {angles.max():.2f} deg (slack at 0)")
print("largest DC flows:")
for br in sorted(dc.branches, key=lambda b: abs(b.p_from_mw), reverse=True)[:5]:
print(f" {br.id:10s} {br.from_bus:9s}->{br.to_bus:9s} {br.p_from_mw:9.2f} MW")
slack_bus = next(b for b in dc.buses if b.role_effective == "slack").id
slack_gens = [g for g in dc.generators if g.bus == slack_bus]
total_load = sum(ld.p_mw for ld in net.loads if ld.in_service)
total_shunt = sum(sh.g_mw for sh in net.shunts if sh.in_service)
total_gen = sum(g.p_mw for g in dc.generators)
print(f"slack bus {slack_bus}: generators {[(g.id, round(g.p_mw, 2)) for g in slack_gens]}")
print(f"generation {total_gen:.2f} MW = load {total_load:.2f} + shunt G {total_shunt:.2f} MW")
# --- Compare with the AC solution ------------------------------------------------------------
ac = pf.solve_ac(net, options=pf.AcOptions(init="flat", q_limits=False))
p_dc = dc.to_arrays().p_from_mw
p_ac = ac.to_arrays().p_from_mw
diff = p_ac - p_dc
worst = int(np.argmax(np.abs(diff)))
print(f"\nAC (flat start, no Q-limits): {ac.iterations} iterations, converged = {ac.converged}")
print(f"AC losses: {sum(b.p_from_mw + b.p_to_mw for b in ac.branches):.2f} MW")
print(f"|P_ac - P_dc| on from-side flows: median {np.median(np.abs(diff)):.2f} MW,", end=" ")
print(f"95th pct {np.percentile(np.abs(diff), 95):.2f} MW, max {np.abs(diff).max():.2f} MW")
print(f"largest gap on {dc.branches[worst].id}: AC {p_ac[worst]:.2f} MW vs DC {p_dc[worst]:.2f} MW")
print("(that branch feeds the slack bus, whose generator picks up every MW of AC losses)")
4. Jobs API¶
pf.ac and pf.dc through jobs.run, the JSON-in / JSON-out path, three failures that come
back as structured results instead of exceptions, and a solver warning captured on the result.
"""The jobs API: one stateless, JSON-serialisable call for every analysis kind.
What this shows:
* ``jobs.SolveRequest(kind, network, options, job_id)`` → ``jobs.run`` → ``jobs.SolveResult``
for ``pf.ac`` and ``pf.dc``; the validated options come back in the provenance.
* ``jobs.run_json`` — text in, text out, exactly what an HTTP handler or a queue worker does —
and the round trip back to typed results.
* Failures are **data**: an unknown kind, bad options and an invalid network each give
``status="failed"`` with a ``StructuredError`` (stable ``code``, message, and the full issue
list for validation), never an exception across the boundary.
* Warnings raised during the solve (here a ``SetpointConflictWarning`` on the case14 variant
with two generators at different setpoints) are captured on the result, not printed.
Run from the repository root: ``uv run python examples/04_jobs_api.py``.
"""
from __future__ import annotations
import json
from mambo_power import jobs
from mambo_power.io import matpower
net = matpower.load("fixtures/matpower/case14.m")
# --- 1. Run pf.ac and pf.dc through the same entry point ----------------------------------
print("registered kinds:", jobs.kinds())
ac = jobs.run(jobs.SolveRequest(kind="pf.ac", network=net, options={"init": "flat"}, job_id="a1"))
dc = jobs.run(jobs.SolveRequest(kind="pf.dc", network=net, job_id="d1"))
for outcome in (ac, dc):
assert outcome.result is not None and outcome.provenance is not None
print(
f"{outcome.kind:5s} job_id={outcome.job_id} status={outcome.status} "
f"result={type(outcome.result).__name__} converged={outcome.result.converged} "
f"slack P={outcome.result.generators[0].p_mw:.3f} MW"
)
print("pf.ac options as run:", ac.provenance.options)
# --- 2. JSON in, JSON out --------------------------------------------------------------------
request_text = jobs.SolveRequest(kind="pf.dc", network=net, job_id="json-1").model_dump_json()
reply_text = jobs.run_json(request_text)
payload = json.loads(reply_text)
print("reply keys:", sorted(payload), "| status:", payload["status"])
typed = jobs.SolveResult.model_validate_json(reply_text) # back to typed models
assert typed.result is not None and dc.result is not None
print("round trip gives", type(typed.result).__name__, end="; ")
print("equal to the direct run:", typed.result.buses == dc.result.buses)
# --- 3. Failures are structured results ----------------------------------------------------
# `market.zonal` stood here until it was registered, at which point this example stopped
# demonstrating an unknown kind and the assertion below failed -- the same day it broke
# `docs/manual/jobs.md` and `tests/unit/test_jobs.py`. `pf.telepathy` is deliberately fictional
# so that cannot happen again: an unknown-kind demo must name a kind that can never become real.
unknown = jobs.run(jobs.SolveRequest(kind="pf.telepathy", network=net))
assert unknown.error is not None
print("\nunknown kind ->", unknown.status, unknown.error.code, "|", unknown.error.message)
bad_options = jobs.run(jobs.SolveRequest(kind="pf.ac", network=net, options={"tol": -1}))
assert bad_options.error is not None and bad_options.error.details is not None
print("bad options ->", bad_options.status, bad_options.error.code, "|", end=" ")
print([(d["loc"], d["type"]) for d in bad_options.error.details])
request = jobs.SolveRequest(kind="pf.ac", network=net)
request.network.branches[0].to_bus = "nowhere" # mutate after construction: no re-validation
invalid = jobs.run(request)
assert invalid.error is not None and invalid.error.issues is not None
print("invalid network ->", invalid.status, invalid.error.code, "|", end=" ")
print([(i.code, i.path) for i in invalid.error.issues])
# --- 4. Warnings travel with the result ----------------------------------------------------
roles = matpower.load("fixtures/matpower/derived/case14_roles.m")
outcome = jobs.run(jobs.SolveRequest(kind="pf.ac", network=roles, options={"init": "flat"}))
assert outcome.result is not None
print("\ncase14_roles ->", outcome.status, "converged", outcome.result.converged)
for line in outcome.warnings:
print(" warning:", line)
5. Roles and islands¶
The derived case14 fixtures: a PV bus solved as PQ, a setpoint conflict resolved by the last generator, a slack without a generator, and an island the importer repairs but the model rejects.
"""Effective bus roles and the island policy on the derived case14 fixtures.
What this shows:
* ``numerics.effective_roles`` derives the roles a solver must use from the declared ones:
a PV bus whose only generator is out of service solves as PQ; a bus with two in-service
generators takes the **last** one's ``v_set_pu`` (MATPOWER's rule) and a
``SetpointConflictWarning`` names the bus when the setpoints differ.
* A slack bus without an in-service generator is a named error, ``NoSlackGeneratorError``.
* Islands: ``load_with_report`` deactivates the buses the slack cannot reach (and their
elements) and reports an ``ISLAND_DEACTIVATED`` issue listing them; the solve then runs on
the main island. The model itself stays strict: constructing the same ``Network`` with the
island switched back on raises ``DISCONNECTED_BUS``.
Run from the repository root: ``uv run python examples/05_roles_and_islands.py``.
"""
from __future__ import annotations
import warnings
from mambo_power import pf
from mambo_power.io import matpower
from mambo_power.model import Network, NetworkValidationError
from mambo_power.numerics import (
NetworkArrays,
NoSlackGeneratorError,
SetpointConflictWarning,
effective_roles,
)
ROLE = {1: "pq", 2: "pv", 3: "slack"}
# --- 1. Effective roles on case14_roles ----------------------------------------------------
net = matpower.load("fixtures/matpower/derived/case14_roles.m")
arr = NetworkArrays.from_network(net)
with warnings.catch_warnings(record=True) as caught:
warnings.simplefilter("always")
roles = effective_roles(arr)
print("case14_roles: declared vs effective role where they differ")
for i, bus_id in enumerate(arr.bus_ids):
if arr.bus_type[i] != roles.bus_type[i]:
print(f" {bus_id}: declared {ROLE[int(arr.bus_type[i])]}, effective", end=" ")
print(f"{ROLE[int(roles.bus_type[i])]} (no in-service generator)")
print("demoted PV buses:", [arr.bus_ids[i] for i in roles.demoted_pv])
for bus_id, gen_ids, setpoints in roles.setpoint_conflicts:
print(f"setpoint conflict at {bus_id}: {gen_ids} -> {setpoints}; last wins:", end=" ")
print(roles.v_set[arr.bus_index[bus_id]])
print("warnings raised:", [type(w.message).__name__ for w in caught])
assert any(issubclass(w.category, SetpointConflictWarning) for w in caught)
with warnings.catch_warnings():
warnings.simplefilter("ignore", SetpointConflictWarning)
result = pf.solve_ac(net, options=pf.AcOptions(init="flat"))
bus6 = next(b for b in result.buses if b.id == "bus-6")
bus2 = next(b for b in result.buses if b.id == "bus-2")
print(f"solved: bus-6 role_effective={bus6.role_effective} vm={bus6.vm_pu:.4f} pu;", end=" ")
print(f"bus-2 vm={bus2.vm_pu:.4f} pu (the last generator's setpoint)")
# --- 2. Slack without a generator --------------------------------------------------------
noslack = matpower.load("fixtures/matpower/derived/case14_noslackgen.m")
try:
pf.solve_ac(noslack)
except NoSlackGeneratorError as err:
print(f"\ncase14_noslackgen: NoSlackGeneratorError for {err.bus_id}: {err}")
# --- 3. Island repair by the importer ------------------------------------------------------
island, report = matpower.load_with_report("fixtures/matpower/derived/case14_island.m")
issue = next(w for w in report.warnings if w.code == "ISLAND_DEACTIVATED")
print("\ncase14_island:", issue.code, "| buses", issue.bus_ids, "| elements", issue.element_ids)
print(" message:", issue.message)
live_buses = sum(b.in_service for b in island.buses)
live_gens = sum(g.in_service for g in island.generators)
print(f" in service: {live_buses} of {len(island.buses)} buses,", end=" ")
print(f"{live_gens} of {len(island.generators)} generators")
main = pf.solve_ac(island, options=pf.AcOptions(init="flat"))
print(f" solve on the main island: converged={main.converged}, {len(main.buses)} bus rows")
# --- 4. The model stays strict --------------------------------------------------------------
raw = island.model_dump()
for bus in raw["buses"]:
if bus["id"] in issue.bus_ids:
bus["in_service"] = True
try:
Network.model_validate(raw)
except NetworkValidationError as err:
print("direct Network with the island re-enabled ->", sorted(err.codes), end=": ")
print([i.path for i in err.issues])
6. Network matrices¶
NetworkArrays, the sparse Ybus, a 100 MW transfer through the PTDF, the LODF's NaN
column on case14's only bridge, and a 3-bus network whose matrices fit on the screen.
"""Network matrices: NetworkArrays, Ybus, PTDF and LODF on case14, then a 3-bus case in full.
What this shows:
* ``NetworkArrays.from_network`` — the positional, per-unit, in-service view that every
matrix builder consumes (the only place positions live).
* ``ybus`` is sparse (CSC): on case14 about 25 % of the entries are non-zero, on case300
about 1 %.
* ``ptdf`` (dense ``n_branch × n_bus``, zero slack column) and a 100 MW transfer through it.
* ``lodf`` with ``bridges``: the single bridge of case14 (branch 7-8) has a ``NaN`` column,
because outaging a bridge islands a bus and no redistribution exists.
* A 3-bus hand example small enough to print every matrix in full.
Run from the repository root: ``uv run python examples/06_network_matrices.py``.
"""
from __future__ import annotations
import numpy as np
from mambo_power.io import matpower
from mambo_power.model import Branch, Bus, Generator, Load, Network
from mambo_power.numerics import NetworkArrays, bbus, bridges, lodf, ptdf, ybus
np.set_printoptions(precision=4, suppress=True, linewidth=100)
# --- 1. case14 -----------------------------------------------------------------------------
net = matpower.load("fixtures/matpower/case14.m")
arr = NetworkArrays.from_network(net)
print(f"case14 arrays: {arr.n_bus} buses, {arr.n_branch} branches, slack position {arr.slack}")
print(f" p_load_pu[:4] = {arr.p_load_pu[:4]} (MW / base_mva = {arr.base_mva})")
y = ybus(arr)
density = y.nnz / (arr.n_bus * arr.n_bus)
print(f"Ybus: {y.shape}, format {y.format}, nnz {y.nnz}, density {density:.1%}")
print(f" diagonal of bus-1: {y[0, 0]:.4f}")
h = ptdf(arr)
print(f"PTDF: shape {h.shape}, slack column all zero: {np.allclose(h[:, arr.slack], 0.0)}")
src, dst = arr.bus_index["bus-5"], arr.bus_index["bus-14"]
flows = 100.0 * (h[:, src] - h[:, dst]) # 100 MW from bus-5 to bus-14
top = np.argsort(-np.abs(flows))[:4]
print(" 100 MW bus-5 -> bus-14, largest branch flows (MW):")
for k in top:
ends = f"{arr.bus_ids[arr.f[k]]:7s}->{arr.bus_ids[arr.t[k]]:7s}"
print(f" {arr.branch_ids[k]:10s} {ends} {flows[k]:8.2f}")
bridge_positions = bridges(arr)
lodf_matrix = lodf(arr, h)
print("bridges:", [(k, arr.branch_ids[k]) for k in bridge_positions], end=" ")
print([f"{arr.bus_ids[arr.f[k]]}-{arr.bus_ids[arr.t[k]]}" for k in bridge_positions])
nan_columns = [k for k in range(arr.n_branch) if np.isnan(lodf_matrix[:, k]).all()]
print(f"LODF: shape {lodf_matrix.shape}; NaN columns {nan_columns} == bridges {bridge_positions}")
print(" (outaging a bridge disconnects a bus: there is no valid redistribution to report)")
k = arr.branch_index["branch-1"]
biggest = np.argsort(-np.nan_to_num(np.abs(lodf_matrix[:, k])))[:3]
print(f" outage of {arr.branch_ids[k]}: largest LODF entries", end=" ")
print([(arr.branch_ids[j], round(float(lodf_matrix[j, k]), 4)) for j in biggest if j != k])
# --- 2. A 3-bus network, every matrix in full ----------------------------------------------
mini = Network(
base_mva=100,
buses=[
Bus(id="b1", base_kv=110, type="slack"),
Bus(id="b2", base_kv=110, type="pv"),
Bus(id="b3", base_kv=110, type="pq"),
],
branches=[
Branch(id="l12", from_bus="b1", to_bus="b2", r=0.01, x=0.10, b=0.02),
Branch(id="l13", from_bus="b1", to_bus="b3", r=0.02, x=0.20, b=0.02),
Branch(id="l23", from_bus="b2", to_bus="b3", r=0.01, x=0.10, b=0.02),
],
generators=[
Generator(
id="g1",
bus="b1",
p_mw=0,
q_mvar=0,
p_min_mw=0,
p_max_mw=300,
q_min_mvar=-100,
q_max_mvar=100,
v_set_pu=1.02,
),
Generator(
id="g2",
bus="b2",
p_mw=60,
q_mvar=0,
p_min_mw=0,
p_max_mw=100,
q_min_mvar=-40,
q_max_mvar=40,
v_set_pu=1.01,
),
],
loads=[Load(id="d3", bus="b3", p_mw=120, q_mvar=40)],
)
small = NetworkArrays.from_network(mini)
print("\n3-bus network, bus order", small.bus_ids, "branch order", small.branch_ids)
print("Ybus (dense):")
print(ybus(small).toarray())
print("B' (DC susceptance):")
print(bbus(small).toarray())
print("PTDF (rows = branches, columns = buses; slack column is 0):")
print(ptdf(small))
print("LODF (no bridges in a triangle, so no NaN column; diagonal -1):")
print(lodf(small))
7. Results and export¶
A result's exact JSON round trip, the positional to_arrays() view, and a CSV export of the
bus and branch tables with the standard library.
"""Results: JSON round trip, the positional ``to_arrays()`` view, and a CSV export.
What this shows:
* A result is a pydantic model: ``model_dump_json`` / ``model_validate_json`` round-trip it
exactly (provenance included), so a result can be stored, queued or returned from a
service as-is.
* ``to_arrays()`` gives one numpy array per column in table order, for numeric consumers
that want positions rather than ids — e.g. finding the most loaded branch or the voltage
envelope in one line.
* Exporting the bus and branch tables to CSV needs nothing beyond the standard library:
every row is a pydantic model, so ``model_dump()`` feeds ``csv.DictWriter`` directly.
Run from the repository root: ``uv run python examples/07_results_and_export.py``.
"""
from __future__ import annotations
import csv
import json
import tempfile
from pathlib import Path
import numpy as np
from mambo_power import pf
from mambo_power.io import matpower
from mambo_power.results import AcPowerFlowResult
net = matpower.load("fixtures/matpower/case14.m")
result = pf.solve_ac(net, options=pf.AcOptions(init="flat"))
# --- 1. JSON round trip --------------------------------------------------------------------
text = result.model_dump_json()
back = AcPowerFlowResult.model_validate_json(text)
print(f"JSON: {len(text)} bytes; round trip equal: {back == result}")
document = json.loads(text)
print("top-level keys:", sorted(document))
print("provenance keys:", sorted(document["provenance"]))
print("a bus row:", document["buses"][0])
# --- 2. The positional view ------------------------------------------------------------------
arrays = result.to_arrays()
print(f"\nto_arrays: {len(arrays.bus_ids)} buses, {len(arrays.branch_ids)} branches,", end=" ")
print(f"{len(arrays.gen_ids)} generators")
lo, hi = int(np.argmin(arrays.vm_pu)), int(np.argmax(arrays.vm_pu))
print(f"voltage envelope: {arrays.bus_ids[lo]} {arrays.vm_pu[lo]:.4f} pu", end=" ... ")
print(f"{arrays.bus_ids[hi]} {arrays.vm_pu[hi]:.4f} pu")
apparent = np.hypot(arrays.p_from_mw, arrays.q_from_mvar)
k = int(np.argmax(apparent))
print(f"largest from-side flow: {arrays.branch_ids[k]} {apparent[k]:.2f} MVA")
print(f"total losses from arrays: {(arrays.p_from_mw + arrays.p_to_mw).sum():.3f} MW")
print("loading_pct is NaN where the branch is unrated:", int(np.isnan(arrays.loading_pct).sum()))
# --- 3. CSV export with the standard library ------------------------------------------------
with tempfile.TemporaryDirectory() as tmp:
out = Path(tmp)
for name, rows in (("buses", result.buses), ("branches", result.branches)):
records = [row.model_dump() for row in rows]
with (out / f"{name}.csv").open("w", newline="", encoding="utf-8") as handle:
writer = csv.DictWriter(handle, fieldnames=list(records[0]))
writer.writeheader()
writer.writerows(records)
for path in sorted(out.iterdir()):
lines = path.read_text(encoding="utf-8").splitlines()
print(f"\n{path.name}: {len(lines) - 1} rows")
print(" " + lines[0])
print(" " + lines[1])
8. OPF and N-1¶
Cost-minimising DC-OPF dispatch and duals on case14, the ac_check AC-feasibility re-solve
finding real voltage violations on an otherwise-clean dispatch, a tightened branch rating
splitting the LMP into energy and congestion, and contingency.n1's LODF screen next to its
confirming DC re-solve on a flagged outage.
"""DC-OPF dispatch and duals, then N-1 branch-contingency screening, on case14.
What this shows:
* ``opf.solve_dc_opf(net, options=...)`` — the cost-minimising LP/QP dispatch, its shadow
prices (the balance dual and each generator's bound reduced cost), and ``options.ac_check``:
a DC-OPF-optimal dispatch is not automatically AC-feasible, and this fixture's own clean
base case already has two buses outside their declared voltage band once AC-solved.
* Locational marginal prices via ``lmp_decomposition``: on case14 as shipped no branch is
rated (no bundled MATPOWER fixture carries a real ``RATE_A``), so every bus's LMP is pure
energy; tightening one branch's rating until it binds splits the price into energy +
congestion, and buses on the constrained side of the network pay more.
* ``contingency.n1`` — the LODF fast screen followed by a confirming DC re-solve, on a copy of
case14 with synthetic ratings (case14 as shipped has none). One outage's screened estimate
and DC-re-solve-confirmed flow, side by side, on a branch the screen correctly flagged.
Run from the repository root: ``uv run python examples/08_opf_and_n1.py``.
"""
from __future__ import annotations
from mambo_power import contingency, opf, pf
from mambo_power.io import matpower
net = matpower.load("fixtures/matpower/case14.m")
# --- 1. DC-OPF dispatch, duals, and the AC-feasibility check --------------------------------
result = opf.solve_dc_opf(net, options=opf.OpfDcOptions(ac_check=True))
print(f"status: {result.status} cost: {result.objective_cost:.2f} $/h", end=" ")
print(f"balance dual (energy price): {result.balance_dual:.4f} $/MWh")
print("dispatch:")
for g in result.generators:
pinned = " (pinned)" if g.bound_dual != 0.0 else ""
print(f" {g.id:8s} {g.bus:8s} {g.p_mw:8.3f} MW bound dual {g.bound_dual:7.4f}{pinned}")
assert result.ac_check is not None # options.ac_check=True guarantees it, once status == Optimal
print(f"\nac_check: converged = {result.ac_check.converged}", end=" ")
print(f"thermal violations: {len(result.ac_check.thermal_violations)}", end=" ")
print(f"voltage violations: {len(result.ac_check.voltage_violations)}")
for v in result.ac_check.voltage_violations:
print(f" {v.bus_id}: {v.vm_pu:.4f} pu vs limit {v.limit_pu:.4f} pu")
print("(the DC-OPF dispatch minimises cost with no voltage constraint at all: an AC re-solve")
print(" of the identical injections can still land outside the declared voltage band)")
# --- 2. Congestion: tighten one branch's rating until the OPF's own dispatch is forced off it -
base = pf.solve_dc(net)
busiest = max(base.branches, key=lambda b: abs(b.p_from_mw))
congested = net.model_copy(deep=True)
tight_rating_mva = abs(busiest.p_from_mw) * 0.5
for br in congested.branches:
if br.id == busiest.id:
br.rating_mva = tight_rating_mva
congested_result = opf.solve_dc_opf(congested)
binding = next(b for b in congested_result.branches if b.id == busiest.id)
print(
f"\n{busiest.id} rated down to {tight_rating_mva:.2f} MVA (from {busiest.p_from_mw:.2f} MW",
"base-case flow)",
)
print(f"congested flow: {binding.p_from_mw:.2f} MW == rating, dual {binding.flow_limit_dual:.4f}")
congested_lmp = [b for b in congested_result.buses if abs(b.congestion) > 1e-9]
print(f"buses with nonzero congestion price: {len(congested_lmp)} of {len(congested_result.buses)}")
for b in sorted(congested_lmp, key=lambda b: b.congestion, reverse=True)[:3]:
print(f" {b.id:8s} lmp {b.lmp:7.3f} = energy {b.energy:7.3f} + congestion {b.congestion:6.3f}")
# --- 3. N-1: LODF screen, then a confirming DC re-solve --------------------------------------
# case14 ships no real branch ratings either (research: RATE_A == 0 everywhere) — derive
# synthetic ones from the base-case flow, the same "test-time transformation of an
# already-owned fixture" this wave's own test suite uses, at 20% headroom above the base flow.
rated = net.model_copy(deep=True)
base_flow_mw = {b.id: abs(b.p_from_mw) for b in base.branches}
for br in rated.branches:
if br.id in base_flow_mw:
br.rating_mva = max(1.2 * base_flow_mw[br.id], 1.0)
n1_result = contingency.n1(rated)
print(
f"\nN-1: {len(n1_result.outages)} outages flagged by the LODF screen (of "
f"{len(rated.branches) - len(n1_result.bridge_branch_ids)} screenable branches)"
)
outage = n1_result.outages[0]
print(f"outage {outage.outage_branch_id}: confirmed violating = {outage.confirmed_violating}")
for flag in outage.flagged_branches[:3]:
print(f" {flag.branch_id:8s} rating {flag.rating_mva:7.2f} MVA ", end="")
print(
f"screened {flag.estimated_flow_mw:7.2f} MW confirmed {flag.confirmed_flow_mw:7.2f} MW",
end=" ",
)
print(f"violating = {flag.confirmed_violating}")
print("(the screen's LODF estimate and the confirming DC re-solve agree to five decimal places:")
print(" AC-6 proves that holds for every outage on every bundled fixture, not just this one)")
9. Nodal market¶
market.solve_nodal on a hand-built 2-bus network wrapped in a Scenario: a bid load whose
elastic response is capped by a binding branch rating, a fixed load with no bid, the LMP split
into energy and congestion, and the settlement identity between load payment, generator
receipts, and congestion rent.
"""Nodal-market clearing: welfare-maximizing DC-OPF with elastic demand, LMPs, and settlement.
What this shows:
* ``market.solve_nodal(scenario)`` on a hand-built 2-bus network -- the same fixture-free style
as the wave's own AC-1 hand-KKT test: a cheap slack generator, an expensive generator behind a
rated branch, one load that bids a 2-segment piecewise-linear demand curve and one that stays
fixed (unbid). ``MarketNodalResult`` reports every load, bid or fixed -- the unbid load still
gets a dispatch/LMP row, just no reduced cost.
* The binding branch rating splits the LMP into energy + congestion, the identical decomposition
``opf.lmp_decomposition`` gives ``opf.solve_dc_opf`` (reused verbatim -- ADR-006, now exercised
by its intended second consumer).
* The settlement identity: total load payment minus total generator receipts equals the
congestion rent, which equals ``-sum(mu_k * flow_k)`` over the binding branches.
Run from the repository root: ``uv run python examples/09_nodal_market.py``.
"""
from __future__ import annotations
from mambo_power import market
from mambo_power.model import (
Branch,
Bus,
Generator,
Load,
Network,
PiecewiseBid,
PolynomialCost,
Scenario,
)
net = Network(
base_mva=100.0,
buses=[
Bus(id="b1", base_kv=138.0, type="slack"),
Bus(id="b2", base_kv=138.0, type="pq"),
],
branches=[
Branch(id="br12", from_bus="b1", to_bus="b2", r=0.0, x=0.1, b=0.0, rating_mva=20.0),
],
generators=[
Generator(
id="g1",
bus="b1",
p_mw=0,
q_mvar=0,
p_min_mw=0,
p_max_mw=100,
q_min_mvar=0,
q_max_mvar=0,
v_set_pu=1.0,
cost=PolynomialCost(coefficients=[10.0, 0.0]), # 10 $/MWh, linear
),
Generator(
id="g2",
bus="b2",
p_mw=0,
q_mvar=0,
p_min_mw=0,
p_max_mw=100,
q_min_mvar=0,
q_max_mvar=0,
v_set_pu=1.0,
cost=PolynomialCost(coefficients=[50.0, 0.0]), # 50 $/MWh, linear
),
],
loads=[
Load(id="d0", bus="b1", p_mw=10.0, q_mvar=0.0), # fixed -- no bid
Load(
id="d1",
bus="b2",
p_mw=100.0,
q_mvar=0.0,
# 2-segment concave PWL bid: marginal value 45 $/MWh on [0, 50], 20 $/MWh on
# [50, 100] (the wave's own hand-KKT example, see AC-1 in the wave spec).
bid=PiecewiseBid(points=[(0.0, 0.0), (50.0, 2250.0), (100.0, 3250.0)]),
),
],
)
result = market.solve_nodal(Scenario(network=net))
print(f"status: {result.status}")
print("dispatch:")
for g in result.generators:
print(f" gen {g.id:4s} bus {g.bus:3s} {g.p_mw:7.3f} MW bound dual {g.bound_dual:7.3f}")
for d in result.loads:
tag = " (fixed, no bid)" if d.id == "d0" else ""
print(f" load {d.id:4s} bus {d.bus:3s} {d.p_mw:7.3f} MW bound dual {d.bound_dual:7.3f}{tag}")
print("LMPs:")
for b in result.buses:
print(f" {b.id}: lmp {b.lmp:7.3f} energy {b.energy:7.3f} congestion {b.congestion:7.3f}")
print(
f"settlement: load payment {result.total_load_payment:.2f} "
f"generator receipts {result.total_generator_receipts:.2f} "
f"congestion rent {result.congestion_rent:.2f}"
)
identity_holds = (
abs((result.total_load_payment - result.total_generator_receipts) - result.congestion_rent)
< 1e-6
)
print(f"settlement identity (payment - receipts == congestion rent) holds: {identity_holds}")
10. Multiperiod market¶
market.solve_multiperiod on a 24-hour horizon over case14 with derived ratings, one storage
unit and ramp limits on every generator: the whole day cleared as one coupled LP, the unit
charging through the overnight trough and discharging into the afternoon peak, the cyclic
end-of-horizon SoC, two binding ramp rows with duals of opposite sign, the per-period settlement
with storage as a third participant, and the periods=None degeneracy reproducing
market.solve_nodal bit-exactly.
"""Multiperiod clearing: a 24-hour horizon with ramp coupling, storage SoC and per-period LMPs.
What this shows:
* ``market.solve_multiperiod(scenario)`` on a 24-period ``Scenario`` built over case14: one
``Period`` per hour carrying a per-load ``load_p_mw`` override, a storage unit, and ramp
limits on every generator. The whole horizon is **one** coupled LP -- the periods are not
solved one at a time, which is the entire point: a ramp row ties hour ``t`` to ``t-1`` and the
state-of-charge rows tie all 24 hours into a single energy budget.
* Storage arbitrage the clearing finds by itself: charge through the overnight trough, discharge
into the afternoon peak, with the charge/discharge efficiencies applied in the SoC row and the
cyclic end-of-horizon condition returning the unit to exactly its starting energy.
* Per-period LMPs, split into energy and congestion by the same ``opf.lmp_decomposition`` a
single-period clearing uses, and the per-period settlement -- where storage is a **third
participant**: it pays ``LMP * charge_mw`` and is paid ``LMP * discharge_mw``, and the
identity does not close if a dispatched unit is left out.
* A ramp row binding, and its dual: negative when the ramp-*up* side binds, positive when the
ramp-*down* side does.
* The degenerate case: a ``Scenario`` with ``periods=None`` clears one period and reproduces
``market.solve_nodal`` exactly -- the same dispatch and the same LMPs, not merely close.
Run from the repository root: ``uv run python examples/10_multiperiod_market.py``.
"""
from __future__ import annotations
import math
from mambo_power import market, pf
from mambo_power.io import matpower
from mambo_power.model import Period, Scenario, Storage
# --- 1. Build the 24-period scenario ----------------------------------------------------------
# case14 ships no branch ratings, no storage and no ramp data (every MATPOWER RATE_A and ramp
# column reads 0, the format's "unpopulated" convention), so all three are derived here from the
# fixture's own committed numbers -- the same test-time-transformation discipline the test suite
# uses, and the same 20%-headroom rating rule as `08_opf_and_n1.py`.
net = matpower.load("fixtures/matpower/case14.m")
base_flow_mw = {b.id: abs(b.p_from_mw) for b in pf.solve_dc(net).branches}
for br in net.branches:
if br.id in base_flow_mw:
br.rating_mva = max(1.2 * base_flow_mw[br.id], 1.0)
total_load_mw = sum(ld.p_mw for ld in net.loads)
load_by_bus: dict[str, float] = {}
for ld in net.loads:
load_by_bus[ld.bus] = load_by_bus.get(ld.bus, 0.0) + ld.p_mw
busiest_bus = max(load_by_bus, key=lambda bus: load_by_bus[bus])
net.storage = [
Storage(
id="st-1",
bus=busiest_bus,
p_max_mw=0.15 * total_load_mw, # 4-hour unit at 15% of system load
energy_mwh=0.15 * total_load_mw * 4.0,
soc_initial=0.5, # half-charged: free to move either way from hour 0
# Deliberately *unequal*: the two efficiencies enter the SoC row with different
# coefficients (+eta_c against -1/eta_d), so with an equal pair transposing them is a
# silent no-op and the asymmetry this example exists to show is invisible. The pair
# also has to leave arbitrage worth doing: this day's LMPs swing 33.31 -> 40.88 $/MWh,
# so a round trip below 33.31/40.88 = 0.815 leaves the unit idle for all 24 hours --
# which is why this is 0.9021 and not `tests/_storage.py`'s more pessimistic 0.8096.
efficiency_charge=0.97,
efficiency_discharge=0.93,
)
]
for g in net.generators:
# `None` means unconstrained; a limit must be strictly > 0 (0 would freeze the unit).
g.ramp_up_mw = g.ramp_down_mw = 0.05 * g.p_max_mw
# A raised-cosine day: every load scaled by the same multiplier, 0.7x at 04:00 up to 1.2x at
# 16:00. `Period.load_p_mw` is an id-keyed *override*, not a scale factor -- a load left out of
# the dict keeps its own `Load.p_mw` in that period.
PEAK, TROUGH, TROUGH_HOUR = 1.2, 0.7, 4
def multiplier(hour: int) -> float:
swing = (1.0 - math.cos(2.0 * math.pi * (hour - TROUGH_HOUR) / 24.0)) / 2.0
return TROUGH + (PEAK - TROUGH) * swing
periods = [
Period(load_p_mw={ld.id: ld.p_mw * multiplier(h) for ld in net.loads}) for h in range(24)
]
result = market.solve_multiperiod(Scenario(network=net, periods=periods))
print(f"status: {result.status} periods: {result.n_periods}", end=" ")
print(f"horizon cost: {result.objective_cost:.2f} $")
# --- 2. The horizon, hour by hour -------------------------------------------------------------
unit = net.storage[0]
print(
f"\nstorage {unit.id} at {unit.bus}: {unit.p_max_mw:.2f} MW / {unit.energy_mwh:.2f} MWh",
end=", ",
)
print(f"round trip {unit.efficiency_charge * unit.efficiency_discharge:.4f}")
print(" h load MW LMP@bus energy congestion charge discharge SoC MWh")
for t, period in enumerate(result.periods):
price = next(b for b in period.buses if b.id == unit.bus)
store = period.storage[0]
print(
f" {t:2d} {sum(ld.p_mw for ld in period.loads):8.2f} {price.lmp:9.4f}"
f" {price.energy:9.4f}"
f" {price.congestion:10.4f} {store.charge_mw:7.3f} {store.discharge_mw:9.3f}"
f" {store.soc_mwh:10.3f}"
)
print(f"cyclic end-of-horizon SoC: {result.periods[-1].storage[0].soc_mwh:.3f} MWh", end=" == ")
print(f"soc_initial * energy_mwh = {unit.soc_initial * unit.energy_mwh:.3f} MWh")
congested = [
t for t, p in enumerate(result.periods) if any(abs(b.congestion) > 1e-9 for b in p.buses)
]
print(f"hours with a binding branch rating: {len(congested)} of 24 -- {congested}")
# --- 3. A binding ramp row and its dual -------------------------------------------------------
binding_ramps = [
(t, g.id, g.ramp_dual)
for t, period in enumerate(result.periods)
for g in period.generators
if abs(g.ramp_dual) > 1e-9
]
print(f"\nbinding ramp rows: {len(binding_ramps)}")
for t, gen_id, dual in binding_ramps:
side = "ramp-up" if dual < 0 else "ramp-down"
previous = next(g for g in result.periods[t - 1].generators if g.id == gen_id)
now = next(g for g in result.periods[t].generators if g.id == gen_id)
limit = next(g for g in net.generators if g.id == gen_id)
print(
f" h{t:02d} {gen_id}: {previous.p_mw:8.3f} -> {now.p_mw:8.3f} MW"
f" (delta {now.p_mw - previous.p_mw:+7.3f}, limit +-{limit.ramp_up_mw:.3f})"
f" {side} dual {dual:.6f} $/MWh"
)
# --- 4. Settlement, per period, with storage as the third participant -------------------------
# `congestion_rent` is the operator's merchandising surplus: (load payment + storage charge
# payment) - (generator receipts + storage discharge revenue). It equals congestion rent
# proper, -sum_k(mu_k * flow_k), only where the network has no phase-shifting transformer and
# no bus shunt conductance -- case14 has neither, so it does here.
print("\nsettlement (per period, $/h):")
print(" h load payment receipts st charge st discharge surplus")
for t in (TROUGH_HOUR, 16):
p = result.periods[t]
print(
f" {t:2d} {p.total_load_payment:12.3f} {p.total_generator_receipts:10.3f}"
f" {p.total_storage_charge_payment:11.3f} {p.total_storage_discharge_revenue:12.3f}"
f" {p.congestion_rent:11.3f}"
)
# An hour with no binding rating has one price everywhere, so the surplus is zero -- to the LP
# solver's own precision, which is what the printed exponent below is: a residual, not a real
# imbalance. It reads that way only because storage is settled. Leaving the two storage columns
# out of the sum (M4's nodal form of the identity, which had no storage to settle) reads a large
# number instead: the arbitrage profit the unit is making at the operator's expense on paper.
uncongested = [p for t, p in enumerate(result.periods) if t not in congested]
worst = max(abs(p.congestion_rent) for p in uncongested)
unsettled = max(abs(p.total_load_payment - p.total_generator_receipts) for p in uncongested)
print(f"largest surplus over the {len(uncongested)} uncongested hours: {worst:.3e} $/h", end=" ")
print(f"(storage left unsettled: {unsettled:.3f} $/h)")
storage_profit = result.total_storage_discharge_revenue - result.total_storage_charge_payment
print(f"horizon: surplus {result.congestion_rent:.3f} $", end=" ")
print(f"storage net revenue {storage_profit:.3f} $ (its arbitrage profit)")
print("(the identity's other side, -sum_k(mu_k*flow_k), needs the flow duals, which this result")
print(" type does not carry; tests/unit/test_market_multiperiod.py computes it from a second,")
print(" array-level solve and proves the equality period by period)")
# --- 5. Degeneracy: one period is the nodal clearing, exactly ----------------------------------
plain = matpower.load("fixtures/matpower/case14.m")
nodal = market.solve_nodal(Scenario(network=plain))
single = market.solve_multiperiod(Scenario(network=plain)) # periods=None -> a one-period horizon
print(f"\nperiods=None -> n_periods {single.n_periods}, status {single.status}")
one = single.periods[0]
same_dispatch = [a.p_mw == b.p_mw for a, b in zip(nodal.generators, one.generators, strict=True)]
same_lmp = [a.lmp == b.lmp for a, b in zip(nodal.buses, one.buses, strict=True)]
print(f"dispatch identical to market.solve_nodal: {all(same_dispatch)} ({len(same_dispatch)} gens)")
print(f"LMPs identical to market.solve_nodal: {all(same_lmp)} ({len(same_lmp)} buses)")
print("(bit-exact `==`, not a tolerance: at T=1 the multiperiod builder issues the identical")
print(" calls, in the identical column and row order, that `dc_opf` itself does)")
11. Zonal redispatch¶
market.solve_zonal on a hand-solvable 2-zone/3-bus market and then on case30 with its three
MATPOWER areas promoted to zones: the corridor at its cap and the price split it creates, the
copper plate the lifted cap produces and the islanding that deleting the corridor produces
instead, a corridor binding in the negative direction with a positive capacity price, the zonal
schedule overloading 17 real branches where the redispatched one overloads none, the three
separated gap figures with the unsigned one negative, and both sides of the settlement identity
computed from the result object alone.
"""Zonal clearing, min-cost redispatch, and what the pair costs against the nodal optimum.
What this shows:
* ``market.solve_zonal(scenario, options)`` -- **three** solves chained, not one: a zonal
clearing that ignores the intra-zone grid, a minimum-cost redispatch that puts the resulting
schedule back onto the real network, and ``market.solve_nodal`` as the reference the pair is
measured against.
* A hand-solvable 2-zone/3-bus market first, where every number can be checked by eye: the
corridor at its cap, two zone prices separated by exactly the two generators' cost difference,
and the same market with the cap lifted, where the two prices collapse into one. Also the
trap: *deleting* the corridor is not the copper plate, it islands the zones.
* case30 with its three MATPOWER areas promoted to real zones and corridor capacities derived
from the cut-set branch ratings: which corridors bind, what each zone pays, and how far the
operator has to move the fleet afterwards.
* The point of the whole exercise -- the zonal schedule overloads real branches, the
redispatched one does not, and the redispatch is what that costs.
* The three separated figures, including the one that is **not** sign-constrained:
``generation_cost_gap`` here is *negative*, and reading it as "zonal beat nodal" is exactly
the mistake it is separated out to prevent.
* Two identities computed from the result object alone: the settlement identity's flow-dual
side (this is the first market result type carrying per-branch duals), and the redispatched
point agreeing with ``market.solve_nodal`` -- which is a theorem, not a coincidence.
Run from the repository root: ``uv run python examples/11_zonal_redispatch.py``.
"""
from __future__ import annotations
from mambo_power import market, pf
from mambo_power.io import matpower
from mambo_power.model import (
Branch,
Bus,
Generator,
Load,
Network,
PolynomialCost,
Scenario,
Zone,
)
from mambo_power.numerics import NetworkArrays
from mambo_power.opf import gen_cost_coeffs
from mambo_power.opf.zonal import zonal_dc_opf
RATING_MARGIN = 1.2 # 20% headroom over the base-case flow, as in 08_opf_and_n1.py
RATING_FLOOR_MVA = 1.0 # so a near-zero base-case flow does not become a near-zero rating
def gen(gen_id: str, bus: str, price: float, p_max: float) -> Generator:
"""A generator offering a flat ``price`` \\$/MWh up to ``p_max`` MW."""
return Generator(
id=gen_id,
bus=bus,
p_mw=0.0,
q_mvar=0.0,
p_min_mw=0.0,
p_max_mw=p_max,
q_min_mvar=0.0,
q_max_mvar=0.0,
v_set_pu=1.0,
cost=PolynomialCost(coefficients=[price, 0.0]),
)
# --- 1. A hand-solvable 2-zone/3-bus market ---------------------------------------------------
# Zone A holds two buses joined by an unrated branch (so the zonal LP is right to carry no
# intra-zone flow row at all); zone B is one bus. The A-B corridor is the only thing that can
# stop cheap zone-A power from serving zone B.
small = Network(
base_mva=100.0,
zones=[Zone(id="A"), Zone(id="B")],
buses=[
Bus(id="bus1", base_kv=138.0, type="slack", zone="A"),
Bus(id="bus2", base_kv=138.0, type="pq", zone="A"),
Bus(id="bus3", base_kv=138.0, type="pq", zone="B"),
],
branches=[
Branch(id="br12", from_bus="bus1", to_bus="bus2", r=0.0, x=0.1, b=0.0),
Branch(id="br23", from_bus="bus2", to_bus="bus3", r=0.0, x=0.1, b=0.0, rating_mva=20.0),
],
generators=[gen("genA", "bus1", 10.0, 200.0), gen("genB", "bus3", 50.0, 200.0)],
loads=[
Load(id="loadA", bus="bus2", p_mw=50.0, q_mvar=0.0),
Load(id="loadB", bus="bus3", p_mw=30.0, q_mvar=0.0),
],
)
small_scenario = Scenario(network=small)
print("=== 1. Two zones, three buses, one corridor ===")
print("genA @ zone A: 10 $/MWh genB @ zone B: 50 $/MWh load: 50 MW in A, 30 MW in B")
small_payment: dict[str, float] = {}
for label, caps in (
("corridor capped at 20 MW", [market.CorridorLimit(zone1="A", zone2="B", cap_mw=20.0)]),
# The copper plate: the corridor stays in the LP with no bound, so the two balance rows
# collapse into one and the market clears as if the zones were one. `cap_mw=None` *is*
# unbounded -- a large finite cap would only be unbounded for a network this small.
("cap lifted (cap_mw=None)", [market.CorridorLimit(zone1="A", zone2="B", cap_mw=None)]),
("no corridor at all", []),
):
res = market.solve_zonal(small_scenario, market.MarketZonalOptions(corridors=caps))
prices = {z.id: z.price for z in res.zones}
schedule = {g.id: g.p_mw for g in res.generators}
small_payment[label] = res.redispatch_payment
print(
f" {label:<26} price A {prices['A']:6.2f} price B {prices['B']:6.2f}"
f" genA {schedule['genA']:6.2f} MW genB {schedule['genB']:6.2f} MW"
)
print(" the 40 $/MWh price split is exactly genB's cost minus genA's -- it is the corridor's")
print(" own capacity shadow price, and it vanishes the moment the corridor stops binding.")
print(" deleting the corridor is NOT the copper plate: with no exchange column the two balance")
print(" rows decouple, each zone self-supplies, and the prices separate as far as they can go.")
print(
" redispatch_payment across the three:"
f" capped {small_payment['corridor capped at 20 MW']:+8.2f}"
f" lifted {small_payment['cap lifted (cap_mw=None)']:+8.2f}"
f" deleted {small_payment['no corridor at all']:+8.2f} $/h"
)
print(" the last one is NEGATIVE: the settlement figure is >= 0 only where the zonal LP is a")
print(" relaxation of the nodal one, i.e. where no corridor cap restricts an exchange more than")
print(" the network itself would. Island the zones and the operator collects instead.")
# --- 2. case30: three areas promoted to zones, corridors from the cut-set ratings --------------
# case30's ZONE column is a single group, but its AREA column carries three real ones. Branch
# ratings are derived from the base-case DC flows (case30's RATE_A is not used here, so there is
# one derivation rule and no mixed provenance), and each corridor's capacity is the sum of the
# ratings on the branches that cross it.
net = matpower.load("fixtures/matpower/case30.m")
base_flow_mw = {b.id: abs(b.p_from_mw) for b in pf.solve_dc(net).branches}
for br in net.branches:
if br.id in base_flow_mw:
br.rating_mva = max(RATING_MARGIN * base_flow_mw[br.id], RATING_FLOOR_MVA)
zone_of_bus = {bus.id: str(bus.area) for bus in net.buses}
net.zones = [Zone(id=zone_id) for zone_id in sorted(set(zone_of_bus.values()))]
for bus in net.buses:
bus.zone = zone_of_bus[bus.id]
caps_mw: dict[tuple[str, str], float] = {}
for br in net.branches:
z1, z2 = zone_of_bus[br.from_bus], zone_of_bus[br.to_bus]
if z1 == z2 or br.rating_mva is None:
continue
key = (min(z1, z2), max(z1, z2))
caps_mw[key] = caps_mw.get(key, 0.0) + br.rating_mva
corridors = [
market.CorridorLimit(zone1=z1, zone2=z2, cap_mw=cap)
for (z1, z2), cap in sorted(caps_mw.items())
]
scenario = Scenario(network=net)
result = market.solve_zonal(scenario, market.MarketZonalOptions(corridors=corridors))
buses_per_zone = {z.id: sum(1 for b in zone_of_bus.values() if b == z.id) for z in net.zones}
print("\n=== 2. case30, three zones ===")
print(f"status: {result.status} buses per zone: {buses_per_zone}")
for corridor in corridors:
crossing = sum(
1
for br in net.branches
if {zone_of_bus[br.from_bus], zone_of_bus[br.to_bus]} == {corridor.zone1, corridor.zone2}
)
print(
f" corridor {corridor.zone1}-{corridor.zone2}: cap {corridor.cap_mw:7.3f} MW"
f" ({crossing} crossing branches)"
)
for zone in result.zones:
print(f" zone {zone.id}: price {zone.price:.6f} $/MWh")
spread = max(z.price for z in result.zones) - min(z.price for z in result.zones)
print(f" price spread across the three zones: {spread:.6f} $/MWh")
# A corridor's own flow and capacity shadow price are array-level quantities: MarketZonalResult
# reports zone prices, not corridor rows. Call the zonal builder directly for them.
arr = NetworkArrays.from_network(net)
cost_coeffs, pwl_costs = gen_cost_coeffs(net, arr)
zonal = zonal_dc_opf(
arr,
cost_coeffs,
{bus_id: zone_of_bus[bus_id] for bus_id in arr.bus_ids},
caps_mw,
pwl_costs=pwl_costs or None,
)
for k, key in enumerate(zonal.corridor_ids):
flow = zonal.corridor_flow_mw[k]
price = zonal.duals.corridor_cap[k]
print(f" corridor {key}: flow {flow:+8.4f} MW capacity price {price:.6f} $/MWh")
print(" corridor (2,3) binds NEGATIVE -- zone 3 exports to zone 2, against the sorted key's own")
print(" direction -- and its capacity price is positive all the same: the price is a magnitude.")
print(" zones 1 and 3 are joined by the one slack corridor, so their balance duals are equal;")
print(" zone 2 separates by exactly the two binding corridors' capacity price.")
# --- 3. The zonal schedule is not deliverable; the redispatched one is -------------------------
# Read both dispatches back through pf.dc and compare each branch flow against its own rating.
# The energy balance is checked too, and deliberately: pf.dc pins the slack bus and lets it
# absorb whatever mismatch the declared injections carry, so a rating-respecting flow vector on
# its own is not proof that a dispatch is feasible.
def overloads(dispatch: dict[str, float]) -> tuple[int, float, float]:
"""(branches over rating, worst overload MW, slack absorption MW) for a generator schedule."""
probe = net.model_copy(deep=True)
for generator in probe.generators:
generator.p_mw = dispatch[generator.id]
solved = pf.solve_dc(probe)
rating = {br.id: br.rating_mva for br in probe.branches}
over = [
abs(b.p_from_mw) - rating[b.id]
for b in solved.branches
if rating[b.id] is not None and abs(b.p_from_mw) > rating[b.id] + 1e-6
]
slack_ids = {b.id for b in solved.buses if b.role_effective == "slack"}
absorbed = sum(g.p_mw - dispatch[g.id] for g in solved.generators if g.bus in slack_ids)
return len(over), max(over, default=0.0), absorbed
zonal_dispatch = {g.id: g.p_mw for g in result.generators}
final_dispatch = {g.id: g.p_mw for g in result.generators_final}
n_zonal, worst_zonal, slack_zonal = overloads(zonal_dispatch)
n_final, worst_final, slack_final = overloads(final_dispatch)
print("\n=== 3. Deliverability ===")
print(
f" zonal schedule: {n_zonal:2d} of {len(net.branches)} branches over rating"
f" (worst {worst_zonal:8.4f} MW) slack absorbs {slack_zonal:+.3e} MW"
)
print(
f" redispatched: {n_final:2d} of {len(net.branches)} branches over rating"
f" (worst {worst_final:8.4f} MW) slack absorbs {slack_final:+.3e} MW"
)
moved_up = sum(g.delta_up_mw for g in result.redispatch_generators)
moved_down = sum(g.delta_down_mw for g in result.redispatch_generators)
touched = sum(1 for g in result.redispatch_generators if g.delta_up_mw + g.delta_down_mw > 1e-9)
print(
f" redispatch volume: +{moved_up:.3f} MW up / -{moved_down:.3f} MW down"
f" across {touched} of {len(net.generators)} generators"
)
# --- 4. The three figures, and the one that is not sign-constrained ----------------------------
print("\n=== 4. What the zonal design cost ===")
print(f" redispatch_payment {result.redispatch_payment:+12.6f} $/h settlement figure")
print(f" welfare_gap {result.welfare_gap:+12.3e} $/h exactness row, 0 by construction")
print(f" generation_cost_gap {result.generation_cost_gap:+12.6f} $/h diagnostic, ANY sign")
print(" the third figure is negative here: the zonal clearing burns less fuel than the nodal")
print(" optimum. It is not therefore cheaper -- it is serving the same demand from a dispatch")
print(" the network cannot carry, and the payment above is what un-carrying it costs.")
print(" the first figure is >= 0 here but not in general -- see part 1's deleted corridor.")
# The three figures are two independent quantities plus a check. Under the theorem below,
# cost(final) == cost(nodal), so generation_cost_gap is exactly minus the payment's fuel term
# and the two published figures sum to the curtailment-compensation term alone -- zero on this
# fixture, which has no elastic demand, and the third field's entire independent content.
compensation = result.redispatch_payment + result.generation_cost_gap
print(f" redispatch_payment + generation_cost_gap = {compensation:+.3e} $/h -- the curtailment")
print(" compensation term, and 0 on this fixed-load fixture: with no bids the third figure")
print(" carries nothing the first does not. Put bids on the same case30 and it is +0.94 $/h.")
nodal = market.solve_nodal(scenario)
nodal_dispatch = {g.id: g.p_mw for g in nodal.generators}
worst_gen = max(abs(final_dispatch[g] - nodal_dispatch[g]) for g in nodal_dispatch)
nodal_lmp = {b.id: b.lmp for b in nodal.buses}
lmp_gaps = sorted((abs(b.lmp - nodal_lmp[b.id]) for b in result.buses), reverse=True)
print(f" redispatched point vs market.solve_nodal: dispatch within {worst_gen:.2e} MW")
print(" (a theorem, not a tolerance sweep: the redispatch objective is the true welfare")
print(" function over nodal's own feasible set, so its optimum IS the nodal optimum)")
# The *primal* theorem above is exact. The duals are a different matter on this fixture: more
# branches sit exactly at their rating than carry a price, so the optimum has several valid
# dual solutions and two LPs may legitimately pick different ones. That is a property of the
# nodal problem, not of either builder -- and it is worth seeing rather than averaging away.
rating_by_id = {br.id: br.rating_mva for br in net.branches}
at_rating = [b.id for b in result.branches if abs(abs(b.p_from_mw) - rating_by_id[b.id]) < 1e-6]
priced = [b.id for b in result.branches if abs(b.flow_limit_dual) > 1e-9]
tight = [gap for gap in lmp_gaps if gap < 1e-4]
print(
f" LMPs: {len(tight)} of {len(lmp_gaps)} buses agree within {max(tight):.1e} $/MWh;"
f" the rest differ by up to {lmp_gaps[0]:.3f} $/MWh"
)
print(
f" because the final point is primal-degenerate: {len(at_rating)} branches sit at their"
f" rating, only {len(priced)} carry a nonzero dual"
)
print(" (put elastic bids on this same fixture and the two solves select the same dual solution")
print(" and every LMP agrees to 1e-5 -- the ambiguity is the nodal LP's, not either builder's)")
# --- 5. Both sides of the settlement identity, from the result object alone --------------------
lmp_by_bus = {b.id: b.lmp for b in result.buses}
load_payment = sum(lmp_by_bus[ld.bus] * ld.p_mw for ld in result.loads_final)
gen_receipts = sum(lmp_by_bus[g.bus] * g.p_mw for g in result.generators_final)
flow_dual_side = -sum(br.flow_limit_dual * br.p_from_mw for br in result.branches)
binding = sum(1 for br in result.branches if abs(br.flow_limit_dual) > 1e-9)
print("\n=== 5. Settlement identity, computed from MarketZonalResult alone ===")
print(f" load payment {load_payment:.4f} - generator receipts {gen_receipts:.4f}", end=" ")
print(f"= {load_payment - gen_receipts:.6f} $/h")
print(f" -sum_k(mu_k * flow_k) over {binding} binding branches = {flow_dual_side:.6f} $/h")
print(f" residual: {abs((load_payment - gen_receipts) - flow_dual_side):.3e} $/h")
print(" no second solve and nothing from numerics/ or opf/ -- MarketZonalResult.branches is")
print(" the first market result surface carrying per-branch flows and their shadow prices.")
12. Strategic bidding¶
market.solve_agents on hand-built linear-cost networks (every bundled MATPOWER generator is
quadratic, and a markup agent needs a linear offer to mark up): the overlay proved by a
byte-identical network after every agent marked up, price-takers reproducing
market.solve_nodal with array_equal on dispatch and LMPs and no tolerance anywhere, a pivotal
supplier's markup climbing to the point where demand's own bid stops paying, checked against the
closed-form optimum, the paired control where a rival rather than demand ends the climb, the
two-agent duopoly reporting converged, the same run under an iteration cap reporting
iteration_cap instead, and the StrategyConfig union crossing jobs as JSON data.
"""Strategic bidding: generators offer, the market clears the offers, and the loop reports how it
ended.
What this shows:
* ``market.agents.solve_agents(scenario, options)`` -- the first market mode whose *input* is an
output of a decision. Every other mode reads the supply curve off the network; here each agent
chooses what to offer, and the offered curve is a different object from the true one.
* The overlay, proved rather than asserted: the network is byte-identical after a run in which
every agent marked up, and ``Generator.cost`` still holds the true cost.
* Price-takers reproduce the competitive result **exactly** -- ``array_equal`` on both dispatch
and LMPs against ``market.solve_nodal``, with no tolerance anywhere. This runs the ordinary
loop; there is no price-taker short-circuit for it to take.
* A pivotal supplier's markup walking up to the point where demand's own bid refuses to pay
more, against the closed-form optimum of the same problem -- and the paired control where a
rival, not demand, is what stops the climb.
* The two-agent duopoly, where ``converged`` has to mean something, and the same run under an
iteration cap, where the result says ``iteration_cap`` rather than pretending to have settled.
* Through ``jobs``: the ``StrategyConfig`` union crossing as JSON data, never a callable.
The synthetic networks are built here rather than imported from ``fixtures/``: a markup agent
needs a **linear** cost, and every one of the 147 generators in the six bundled MATPOWER cases
carries a quadratic one. The price-taker section below does use ``case14``, because a price-taker
offers whatever shape its true cost is.
Run from the repository root: ``uv run python examples/12_agent_market.py``.
"""
from __future__ import annotations
import json
import numpy as np
from mambo_power import jobs, market
from mambo_power.io import matpower
from mambo_power.market.agents import MarketAgentsOptions, solve_agents
from mambo_power.model import (
Branch,
Bus,
Generator,
Load,
Network,
PolynomialBid,
PolynomialCost,
Scenario,
)
# The shared demand side of every synthetic market below: marginal value ``100 - 0.1*p``, i.e.
# ``q(price) = 1000 - 10*price``. A ``PolynomialBid``'s coefficients are the *value* curve, whose
# derivative is the marginal value, so v1 = 100.0 and v2 = -0.05. ``p_mw`` is the quantity at
# which marginal value reaches zero -- a smaller cap would truncate the curve before the market
# reached its own optimum.
DEMAND_BID = PolynomialBid(coefficients=[-0.05, 100.0, 0.0])
DEMAND_P_MAX_MW = 1000.0
def star(generators: list[tuple[str, float, float]]) -> Network:
"""A star network: ``b1`` (slack) hosts the first generator, every further one gets its own
bus, and the shared elastic load sits on the last bus.
Each entry is ``(id, p_max_mw, true_marginal_cost)``. Every branch is built with no rating,
so nothing here ever congests -- this example is about bidding, not about flow limits. Every
true cost is linear, ``cost(p) = c1 * p``, which is both what ``MarkupStrategy`` requires and
what makes the profit arithmetic below closed-form.
"""
n = len(generators)
buses = [Bus(id="b1", base_kv=138.0, type="slack")]
buses += [Bus(id=f"b{i}", base_kv=138.0, type="pq") for i in range(2, n + 2)]
branches = [
Branch(id=f"l{i}", from_bus="b1", to_bus=f"b{i}", r=0.0, x=0.05, b=0.0)
for i in range(2, n + 2)
]
gens = [
Generator(
id=gen_id,
bus="b1" if k == 0 else f"b{k + 1}",
p_mw=0.0,
q_mvar=0.0,
p_min_mw=0.0,
p_max_mw=p_max_mw,
q_min_mvar=-9999.0,
q_max_mvar=9999.0,
v_set_pu=1.0,
cost=PolynomialCost(coefficients=[true_cost, 0.0]),
)
for k, (gen_id, p_max_mw, true_cost) in enumerate(generators)
]
load = Load(id="d1", bus=f"b{n + 1}", p_mw=DEMAND_P_MAX_MW, q_mvar=0.0, bid=DEMAND_BID)
return Network(base_mva=100.0, buses=buses, branches=branches, generators=gens, loads=[load])
def price(result) -> float: # noqa: ANN001 - any market result carrying `buses`
"""The clearing price at the load's own bus, which is the last one every `star` builds."""
return result.buses[-1].lmp
# --- 1. Price-takers reproduce the competitive result, exactly --------------------------------
# Every generator on case14 offers its own true cost, unchanged. This is an ordinary run of the
# loop -- the offer map is built, the overlay is handed to the array builder, the clearing comes
# back and is compared -- and it agrees with `market.solve_nodal` bitwise, not to a tolerance.
case14 = matpower.load("fixtures/matpower/case14.m")
scenario14 = Scenario(network=case14)
taker_options = MarketAgentsOptions(
strategies={gen.id: {"kind": "price_taker"} for gen in case14.generators}
)
taker = solve_agents(scenario14, taker_options)
nodal = market.solve_nodal(scenario14)
print("--- 1. price-takers vs market.solve_nodal, on case14 ---")
print(
"dispatch array_equal:",
np.array_equal(
np.array([g.p_mw for g in taker.generators]),
np.array([g.p_mw for g in nodal.generators]),
),
"| LMP array_equal:",
np.array_equal(np.array([b.lmp for b in taker.buses]), np.array([b.lmp for b in nodal.buses])),
)
print(
f"status {taker.status} | converged {taker.converged} | "
f"termination_reason {taker.termination_reason} | iterations {taker.iterations}"
)
print(
"every offer is the true cost object:",
all(offer.offer == offer.true_cost for offer in taker.offers),
"| markups:",
sorted({offer.markup for offer in taker.offers}),
)
# --- 2. A pivotal supplier climbs to demand's own limit ----------------------------------------
# One 900 MW unit at a true $20/MWh, no rival, facing q = 1000 - 10*price. Profit
# (pi - 20)(1000 - 10*pi) peaks at pi = $60.00, q = 400 MW, $16,000/h -- a closed form this
# market has no knowledge of. The agent finds it by climbing on its own observed profit alone.
pivotal = star([("strategic", 900.0, 20.0)])
before_json = pivotal.model_dump_json()
markup_options = MarketAgentsOptions(
strategies={"strategic": {"kind": "markup", "step": 0.5}}, offer_tol=1.5
)
climbed = solve_agents(Scenario(network=pivotal), markup_options)
baseline = solve_agents(
Scenario(network=pivotal),
MarketAgentsOptions(strategies={"strategic": {"kind": "price_taker"}}),
)
peak = climbed.offers[0]
print()
print("--- 2. a pivotal supplier, against a closed-form optimum ---")
print("closed form: offer $60.00/MWh, cleared 400.00 MW, profit $16,000.00/h")
print(
f"the climb: offer ${peak.offer.coefficients[0]:.2f}/MWh, "
f"cleared {peak.cleared_mw:.2f} MW, markup ${peak.markup:,.2f}/h"
)
print(
f"at true cost: price ${price(baseline):.2f}/MWh, "
f"cleared {baseline.offers[0].cleared_mw:.2f} MW, markup ${baseline.offers[0].markup:,.2f}/h"
)
print(f"clearing price ${price(climbed):.2f}/MWh after {climbed.iterations} update rounds")
# The overlay never touched the network. Both halves matter: byte-identity alone would also hold
# for a run in which nothing happened, and the markup above is what rules that out.
print("network byte-identical after the run:", pivotal.model_dump_json() == before_json)
print("Generator.cost still the true curve:", pivotal.generators[0].cost.coefficients)
# --- 3. The paired control: a rival, not demand, stops the climb -------------------------------
# The same unit, now with a 900 MW rival at $22/MWh. The markup is real and nonzero -- market
# power is reduced, not eliminated -- but it is an order of magnitude smaller, and what stops it
# is the rival's cost rather than demand's willingness to pay.
controlled = solve_agents(
Scenario(network=star([("strategic", 900.0, 20.0), ("rival", 900.0, 22.0)])), markup_options
)
rivalled = controlled.offers[0]
print()
print("--- 3. the same agent with a rival at $22/MWh ---")
print(
f"offer ${rivalled.offer.coefficients[0]:.2f}/MWh, cleared {rivalled.cleared_mw:.2f} MW, "
f"markup ${rivalled.markup:,.2f}/h after {controlled.iterations} update rounds"
)
print(f"against the pivotal ${peak.markup:,.2f}/h -- {peak.markup / rivalled.markup:.1f}x smaller")
# --- 4. Two agents, and what `converged` has to mean -------------------------------------------
# Two 300 MW units at $20/MWh: the only shape in this example where best response can fail to
# settle in one round. A fixed-step climber never comes to rest -- it dithers by two steps about
# its optimum, three when the optimum sits halfway between two grid points -- so the loop
# classifies the repetition it finds by its *amplitude*, which is why `offer_tol` must be at
# least `3 * step` and why the options model enforces that.
duopoly = star([("g1", 300.0, 20.0), ("g2", 300.0, 20.0)])
duopoly_strategies = {
"g1": {"kind": "markup", "step": 0.5},
"g2": {"kind": "markup", "step": 0.5},
}
settled = solve_agents(
Scenario(network=duopoly),
MarketAgentsOptions(strategies=duopoly_strategies, offer_tol=1.5),
)
competitive = solve_agents(
Scenario(network=duopoly),
MarketAgentsOptions(strategies={"g1": {"kind": "price_taker"}, "g2": {"kind": "price_taker"}}),
)
print()
print("--- 4. a two-agent duopoly ---")
print(
"offers",
[offer.offer.coefficients[0] for offer in settled.offers],
f"| price ${price(settled):.2f}/MWh | joint markup "
f"${sum(offer.markup for offer in settled.offers):,.2f}/h",
)
print(
f"at true cost: price ${price(competitive):.2f}/MWh, "
f"cleared {[round(offer.cleared_mw, 2) for offer in competitive.offers]}"
)
print(
f"status {settled.status} | converged {settled.converged} | "
f"termination_reason {settled.termination_reason} | iterations {settled.iterations}"
)
# The same run under a cap it cannot meet. `status` is still the LP's and still Optimal; the loop
# reports that it ran out of rounds, and never presents a truncated run as a settled one.
capped = solve_agents(
Scenario(network=duopoly),
MarketAgentsOptions(strategies=duopoly_strategies, offer_tol=1.5, max_iterations=10),
)
print(
f"under max_iterations=10: status {capped.status} | converged {capped.converged} | "
f"termination_reason {capped.termination_reason} | iterations {capped.iterations}"
)
# A tolerance narrower than the settling oscillation would turn every arrival into a false cycle
# report, so it is rejected up front rather than mis-diagnosed later.
try:
MarketAgentsOptions(strategies={"g1": {"kind": "markup", "step": 0.5}}, offer_tol=0.5)
except ValueError as exc:
print("offer_tol below 3 * step is refused:", str(exc).splitlines()[1].strip()[:78])
# --- 5. Through the jobs API -------------------------------------------------------------------
# The strategy configuration crosses as data. A `Strategy` object never does: `solve_agents` has
# an in-process `strategies=` seam for a rule the config union cannot express, and `jobs` cannot
# reach it, so nothing a service sends decides which code runs.
request = jobs.SolveRequest(
kind="market.agents",
network=duopoly,
options={"strategies": duopoly_strategies, "offer_tol": 1.5},
)
reply = json.loads(jobs.run_json(request.model_dump_json()))
print()
print("--- 5. through jobs ---")
print("kinds:", jobs.kinds())
print(
reply["status"],
reply["provenance"]["kind"],
"| converged",
reply["result"]["converged"],
"| termination_reason",
reply["result"]["termination_reason"],
"| iterations",
reply["result"]["iterations"],
)
print("strategies crossed JSON as data:", reply["provenance"]["options"]["strategies"])
bad = jobs.run(
jobs.SolveRequest(
kind="market.agents",
network=duopoly,
options={"strategies": {"nope": {"kind": "markup", "step": 0.5}}, "offer_tol": 1.5},
)
)
print("a strategy naming a generator that does not exist:", bad.status, bad.error.code)
13. Interop¶
Unlike the other twelve, this script takes about a minute rather than a second: it imports pandapower and PyPSA (their cold imports alone are ~20 s) and runs their solvers as the oracles.
One Network (IEEE case14) through every format of wave M8, each conversion returning the
report that says what it could not carry: the pandapower JSON export loaded by pp.from_json
and solved by pandapower's own rundcpp, agreeing with pf.solve_dc to 1e-14 degrees;
pp.networks.case14() imported with an empty report and its neutral-tap transformers kept
as transformers because the source table says so; the PyPSA export optimised by PyPSA, its
objective agreeing with opf.solve_dc_opf to 1e-13 relative (the constant cost term travels
in the marginal_cost_constant column); fixtures/case14_v33.raw imported with
RAW_NO_COSTS rather than invented costs; the CSV bundle round trip load(dump(net)) == net
with no tolerance; and a piecewise-cost generator pushed into PyPSA, which has no piecewise
cost, with the ExportReport naming the generator and what was written instead.
"""Interchange: one ``Network`` in and out of pandapower, PyPSA, PSS/E RAW and a CSV bundle,
with every conversion reporting what it could not carry.
What this shows:
* ``io.pandapower_json.dumps(net)`` -> ``pp.from_json_string`` -> ``pp.rundcpp``: pandapower's
own DC solver on the exported document agrees with ``pf.solve_dc`` on the original, angle by
angle, and the export report names each field pandapower has no column for.
* ``io.pandapower_json.loads_with_report`` on ``pp.networks.case14()``: pandapower's case14 comes
in with an **empty** report -- the conversion was lossless -- and its neutral-tap transformers
keep ``kind="transformer"`` because the source table says so.
* ``io.pypsa.to_network_with_report`` then PyPSA ``optimize()``: the DC-OPF objective agrees
with ``opf.solve_dc_opf`` (the constant cost term travels beside PyPSA's objective in the
``marginal_cost_constant`` column, since ``n.objective`` excludes constants).
* ``io.psse_raw.load_with_report`` on ``fixtures/case14_v33.raw``: the same IEEE case14 spelled
as a RAW file; RAW carries no costs, and the report says so with ``RAW_NO_COSTS`` rather than
the importer inventing any.
* ``io.csv_bundle.dump`` / ``load`` through a temporary directory: ``load(dump(net)) == net``,
bit-exact, with no tolerance.
* One deliberately lossy conversion -- a piecewise-cost generator into PyPSA, which has no
piecewise cost -- and the ``ExportReport`` that names the generator, the field and what was
written instead.
Every conversion is *best effort + report*: an empty report means lossless; anything dropped,
approximated or repaired is an issue naming the element id and the field. Nothing is logged or
printed by the converters themselves.
pandapower and PyPSA are development extras (``uv sync`` installs them); the core package never
imports either.
Run from the repository root: ``uv run python examples/13_interop.py``.
"""
from __future__ import annotations
import logging
import tempfile
import warnings
from pathlib import Path
import numpy as np
import pandapower as pp
import pandapower.networks as pn
from mambo_power import opf, pf
from mambo_power.io import csv_bundle, matpower, pandapower_json, psse_raw, pypsa
from mambo_power.model import PiecewiseCost
# The third-party libraries log freely (numba advice, PyPSA consistency notes, solver banners);
# the converters themselves never do. Keep the output to what this script prints.
for name in ("pandapower", "pypsa", "linopy"):
logging.getLogger(name).setLevel(logging.ERROR)
net = matpower.load("fixtures/matpower/case14.m")
print(f"case14: {len(net.buses)} buses, {len(net.branches)} branches, {len(net.generators)} gens")
# --- 1. pandapower JSON export, solved by pandapower itself ----------------------------------
text, export_report = pandapower_json.dumps_with_report(net)
print(f"\npandapower export: {len(text)} chars, report codes {sorted(export_report.codes)}")
for issue in export_report.warnings[:3]:
print(" ", issue)
print(f" ... {len(export_report.warnings)} issues in all; none touches a carried value")
pp_net = pp.from_json_string(text)
with warnings.catch_warnings(): # pandapower warns about the missing optional numba
warnings.simplefilter("ignore")
pp.rundcpp(pp_net, numba=False, trafo_model="pi")
ours = pf.solve_dc(net)
theirs = {
str(name): float(va) for name, va in zip(pp_net.bus.name, pp_net.res_bus.va_degree, strict=True)
}
worst = max(abs(b.va_deg - theirs[b.id]) for b in ours.buses)
print(f"pp.rundcpp vs pf.solve_dc: worst angle difference {worst:.1e} deg, {len(ours.buses)} buses")
# --- 2. pandapower import: pp.networks.case14() ---------------------------------------------
pp14, import_report = pandapower_json.loads_with_report(pp.to_json(pn.case14()))
trafos = [br.id for br in pp14.branches if br.kind == "transformer"]
print(f"\npandapower import of pp.networks.case14(): report {import_report.as_strings()}")
print(" (an empty report means the conversion was lossless)")
print(f" {len(pp14.buses)} buses, {len(pp14.branches)} branches, transformers {trafos}")
neutral = [
br.id for br in pp14.branches if br.kind == "transformer" and br.tap_ratio in (None, 1.0)
]
print(f" neutral-tap transformers kept as transformers by the source table: {neutral}")
# --- 3. PyPSA export, optimised by PyPSA ------------------------------------------------------
n, pypsa_report = pypsa.to_network_with_report(net)
print(f"\nPyPSA export: {len(n.buses)} buses, {len(n.lines)} lines, {len(n.transformers)} trafos")
print(f" report codes {sorted(pypsa_report.codes)} ({len(pypsa_report.warnings)} issues)")
with warnings.catch_warnings():
warnings.simplefilter("ignore")
status = n.optimize(solver_name="highs", solver_options={"output_flag": False})
c0 = float(n.generators[pypsa.COST_CONSTANT_COLUMN].sum())
pypsa_objective = float(n.objective) + c0
dc_opf = opf.solve_dc_opf(net)
rel = abs(dc_opf.objective_cost - pypsa_objective) / pypsa_objective
print(f" PyPSA optimize {status}: objective {pypsa_objective:.4f} $/h (incl. constant {c0:.1f})")
print(f" opf.solve_dc_opf {dc_opf.status}: objective {dc_opf.objective_cost:.4f} $/h")
print(f" relative difference {rel:.1e}")
# --- 4. PSS/E RAW v33 import ------------------------------------------------------------------
raw_net, raw_report = psse_raw.load_with_report("fixtures/case14_v33.raw")
print(f"\nRAW import: {len(raw_net.buses)} buses, {len(raw_net.branches)} branches")
print(f" report codes {sorted(raw_report.codes)} ({len(raw_report.warnings)} issues)")
print(" ", next(str(w) for w in raw_report.warnings if w.code == "RAW_NO_COSTS"))
raw_dc = pf.solve_dc(raw_net)
raw_angles = np.array([b.va_deg for b in raw_dc.buses])
m_angles = np.array([b.va_deg for b in ours.buses])
raw_worst = float(np.abs(raw_angles - m_angles).max())
print(
f" pf.solve_dc on the RAW network vs the MATPOWER one: worst angle diff {raw_worst:.1e} deg"
)
# --- 5. CSV bundle round trip -----------------------------------------------------------------
with tempfile.TemporaryDirectory() as directory:
csv_bundle.dump(net, directory)
files = sorted(p.name for p in Path(directory).iterdir())
back = csv_bundle.load(directory)
print(f"\nCSV bundle: {files}")
print(f" load(dump(net)) == net: {back == net}")
# --- 6. A deliberately lossy conversion --------------------------------------------------------
lossy = net.model_copy(deep=True)
lossy.generators[1].cost = PiecewiseCost(points=[(0.0, 0.0), (50.0, 1500.0), (140.0, 5000.0)])
_, lossy_report = pypsa.to_network_with_report(lossy)
print("\nPiecewise cost into PyPSA (which has none):")
for issue in lossy_report.warnings:
if issue.code == "PYPSA_PWL_COST_DROPPED":
print(f" {issue.code}: element_ids={issue.element_ids}")
print(f" {issue.message}")
Conventions for examples¶
- Each script is self-contained, runs from the repository root, and reads only files under
fixtures/. It prints what it computes and writes nothing outside a temporary directory. - Numbered
NN_name.py; the module docstring says what the script shows and how to run it. - A script that exits non-zero fails
tests/unit/test_examples_run.pyand theexamplesCI job; a script that is not embedded on this page fails the same test. - Embedding uses the
{ .python }fence form on purpose:ruff formatrewrites Python fences in Markdown and would turn the--8<--marker into an expression.