Flow-Based Market Coupling¶

Flow-based market coupling (FBMC) is how much of Europe allocates cross-border capacities of electricity trade. This notebook builds it up from scratch on a small three-zone system in PyPSA, in three parts:

The domain - what the flow-based feasible region is and what it represents.
Market clearing - letting PyPSA optimise the market clearing with a flow-based domain.
External borders - ways to handle borders that lie outside the flow-based region.

The core idea¶

A market is split into zones (also called hubs or bidding zones). Each zone has a net position $NP_z$: its total generation minus its load, i.e. its net export (positive) or import (negative), in MW. The grid cannot host every combination of net positions; FBMC describes the ones it can host with a set of linear inequalities:

$$\sum_{z} \mathrm{PTDF}_{c,z}\, NP_z \le \mathrm{RAM}_c \quad\text{for every constraint } c, \qquad\qquad \sum_{z} NP_z = 0.$$

CNEC (Critical Network Element + Contingency) - a monitored grid element, in one flow direction, under one assumed outage. One CNEC = one constraint $c$ = one row.
Zonal PTDF (Power Transfer Distribution Factor) - how strongly zone $z$'s export loads CNEC $c$. A value of $0.1$ means 10% of the export flows over that element; a negative value means the export relieves it.
RAM (Remaining Available Margin) - the MW of headroom left on the CNEC for the market to use.
The second equation is energy conservation: exports and imports must net to zero.

Together these inequalities carve out the flow-based domain: the set of net-position combinations the grid can physically support. Throughout, we take the domain (PTDF + RAM) as a given input. Deriving it from a grid model is a separate task.

In [2]:

Copied!





import linopy
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.spatial import HalfspaceIntersection

import pypsa
import linopy
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from scipy.spatial import HalfspaceIntersection

import pypsa

Part 1 - The flow-based domain¶

Three zones A, B, C sit at the corners of a triangle joined by three lines AB, BC, AC of equal reactance. Each line is monitored in both directions, so the domain has six CNECs. With three zones and $\sum_z NP_z = 0$, only two net positions are free - the domain is a 2-D polygon we can draw.

The PTDFs are the textbook DC-flow values for an equal-reactance triangle with C as the reference zone (its column is zero). A unit export from A, withdrawn at C, splits 2/3 over the direct line AC and 1/3 over the detour AB+BC.

In [3]:

Copied!





zones = ["A", "B", "C"]
lines = {
    "AB": {"A": 1 / 3, "B": -1 / 3, "C": 0.0, "RAM": 1000.0},
    "BC": {"A": 1 / 3, "B": 2 / 3, "C": 0.0, "RAM": 1500.0},
    "AC": {"A": 2 / 3, "B": 1 / 3, "C": 0.0, "RAM": 2000.0},
}
pos = pd.DataFrame(lines).T
neg = pos.copy()
neg[zones] *= -1  # reverse direction: flip PTDF signs
domain = pd.concat(
    [pos.rename(lambda s: f"{s} (+)"), neg.rename(lambda s: f"{s} (-)")]
).sort_index()
domain.index.name = "cnec"
domain
zones = ["A", "B", "C"]
lines = {
    "AB": {"A": 1 / 3, "B": -1 / 3, "C": 0.0, "RAM": 1000.0},
    "BC": {"A": 1 / 3, "B": 2 / 3, "C": 0.0, "RAM": 1500.0},
    "AC": {"A": 2 / 3, "B": 1 / 3, "C": 0.0, "RAM": 2000.0},
}
pos = pd.DataFrame(lines).T
neg = pos.copy()
neg[zones] *= -1  # reverse direction: flip PTDF signs
domain = pd.concat(
    [pos.rename(lambda s: f"{s} (+)"), neg.rename(lambda s: f"{s} (-)")]
).sort_index()
domain.index.name = "cnec"
domain

Out[3]:

	A	B	C	RAM
cnec
AB (+)	0.333333	-0.333333	0.0	1000.0
AB (-)	-0.333333	0.333333	-0.0	1000.0
AC (+)	0.666667	0.333333	0.0	2000.0
AC (-)	-0.666667	-0.333333	-0.0	2000.0
BC (+)	0.333333	0.666667	0.0	1500.0
BC (-)	-0.333333	-0.666667	-0.0	1500.0

Reading one row. AC (+) has $\mathrm{PTDF}_A = 0.67$, $\mathrm{PTDF}_B = 0.33$, $\mathrm{RAM} = 2000$, i.e. $0.67\,NP_A + 0.33\,NP_B \le 2000$. If A exports 3000 MW with B balanced, the loading is $0.67 \times 3000 = 2000$, exactly the limit. A is then pinned by line AC, unable to export more even while other lines still have room. That coupling between zones is what FBMC captures and a single number per border cannot.

Drawing the domain¶

Substituting $NP_C = -(NP_A + NP_B)$ leaves $(NP_A, NP_B)$ as the two free axes; each CNEC becomes a straight edge and the feasible region is their intersection. The reduction is general: the effective coefficient on a free zone $z$ is $\mathrm{PTDF}_{c,z} - \mathrm{PTDF}_{c,\mathrm{ref}}$, where $\mathrm{ref}$ is the chosen reference zone. The picture does not depend on which zone that is.

In [4]:

Copied!





ref = "C"
free = [z for z in zones if z != ref]
a_ub = domain[free].values - domain[[ref]].values  # half-plane coefficients
b_ub = domain["RAM"].values
hs = HalfspaceIntersection(np.column_stack([a_ub, -b_ub]), np.zeros(len(free)))
v = hs.intersections
poly = v[np.argsort(np.arctan2(v[:, 1] - v[:, 1].mean(), v[:, 0] - v[:, 0].mean()))]


def draw_domain(ax: plt.Axes) -> None:
    """Draw the feasible polygon, its CNEC reference lines and edge labels."""
    span = np.abs(poly).max() * 1.4
    t = np.array([-span, span])
    for a1, a2, rhs in zip(a_ub[:, 0], a_ub[:, 1], b_ub, strict=True):
        xy = (t, (rhs - a1 * t) / a2) if abs(a2) > abs(a1) else ((rhs - a2 * t) / a1, t)
        ax.plot(*xy, color="0.75", lw=0.8, zorder=0)
    ax.fill(poly[:, 0], poly[:, 1], color="#440154", alpha=0.12)
    ring = np.vstack([poly, poly[0]])
    ax.plot(ring[:, 0], ring[:, 1], color="#440154", lw=1.5)
    for p, q in zip(poly, ring[1:], strict=True):
        mid = (p + q) / 2
        ax.annotate(
            domain.index[np.argmin(b_ub - a_ub @ mid)],
            mid,
            textcoords="offset points",
            xytext=(4, 4),
            fontsize=8,
            color="#440154",
        )
    ax.axhline(0, linestyle="--", color="k", lw=0.6)
    ax.axvline(0, linestyle="--", color="k", lw=0.6)
    ax.set(
        xlim=(-span, span),
        ylim=(-span, span),
        aspect="equal",
        xlabel=f"NP_{free[0]}  (MW, export +)",
        ylabel=f"NP_{free[1]}  (MW, export +)",
    )


def check(np_free: tuple[float, float]) -> pd.DataFrame:
    """Return loadings, margins and active (binding) CNECs for a net position."""
    netpos = pd.Series(dict(zip(free, np_free, strict=True)))
    netpos[ref] = -netpos.sum()
    load = domain[zones] @ netpos
    out = pd.DataFrame(
        {
            "loading": load.round(1),
            "RAM": domain["RAM"],
            "margin": (domain["RAM"] - load).round(1),
        }
    )
    out["active"] = out["margin"].abs() < 1e-6
    return out
ref = "C"
free = [z for z in zones if z != ref]
a_ub = domain[free].values - domain[[ref]].values  # half-plane coefficients
b_ub = domain["RAM"].values
hs = HalfspaceIntersection(np.column_stack([a_ub, -b_ub]), np.zeros(len(free)))
v = hs.intersections
poly = v[np.argsort(np.arctan2(v[:, 1] - v[:, 1].mean(), v[:, 0] - v[:, 0].mean()))]


def draw_domain(ax: plt.Axes) -> None:
    """Draw the feasible polygon, its CNEC reference lines and edge labels."""
    span = np.abs(poly).max() * 1.4
    t = np.array([-span, span])
    for a1, a2, rhs in zip(a_ub[:, 0], a_ub[:, 1], b_ub, strict=True):
        xy = (t, (rhs - a1 * t) / a2) if abs(a2) > abs(a1) else ((rhs - a2 * t) / a1, t)
        ax.plot(*xy, color="0.75", lw=0.8, zorder=0)
    ax.fill(poly[:, 0], poly[:, 1], color="#440154", alpha=0.12)
    ring = np.vstack([poly, poly[0]])
    ax.plot(ring[:, 0], ring[:, 1], color="#440154", lw=1.5)
    for p, q in zip(poly, ring[1:], strict=True):
        mid = (p + q) / 2
        ax.annotate(
            domain.index[np.argmin(b_ub - a_ub @ mid)],
            mid,
            textcoords="offset points",
            xytext=(4, 4),
            fontsize=8,
            color="#440154",
        )
    ax.axhline(0, linestyle="--", color="k", lw=0.6)
    ax.axvline(0, linestyle="--", color="k", lw=0.6)
    ax.set(
        xlim=(-span, span),
        ylim=(-span, span),
        aspect="equal",
        xlabel=f"NP_{free[0]}  (MW, export +)",
        ylabel=f"NP_{free[1]}  (MW, export +)",
    )


def check(np_free: tuple[float, float]) -> pd.DataFrame:
    """Return loadings, margins and active (binding) CNECs for a net position."""
    netpos = pd.Series(dict(zip(free, np_free, strict=True)))
    netpos[ref] = -netpos.sum()
    load = domain[zones] @ netpos
    out = pd.DataFrame(
        {
            "loading": load.round(1),
            "RAM": domain["RAM"],
            "margin": (domain["RAM"] - load).round(1),
        }
    )
    out["active"] = out["margin"].abs() < 1e-6
    return out

Three example net positions are marked below:

Case 1 NP = (+1000, -500, -500) - well inside the domain; every CNEC keeps spare margin.
Case 2 NP = (+1500, -1500, 0) - a straight transfer from A to B landing on an edge (one binding CNEC).
Case 3 NP = (+3000, 0, -3000) - A exporting to C, landing at a corner (two binding CNECs).

In [5]:

Copied!





case1, case2, case3 = (1000, -500), (1500, -1500), (3000, 0)

fig, ax = plt.subplots(figsize=(5.5, 5.5), layout="constrained")
draw_domain(ax)
for pt, name in [(case1, "Case 1"), (case2, "Case 2"), (case3, "Case 3")]:
    ax.scatter(*pt, color="C3", zorder=3)
    ax.annotate(
        name, pt, textcoords="offset points", xytext=(6, -12), fontsize=9, color="C3"
    )
case1, case2, case3 = (1000, -500), (1500, -1500), (3000, 0)

fig, ax = plt.subplots(figsize=(5.5, 5.5), layout="constrained")
draw_domain(ax)
for pt, name in [(case1, "Case 1"), (case2, "Case 2"), (case3, "Case 3")]:
    ax.scatter(*pt, color="C3", zorder=3)
    ax.annotate(
        name, pt, textcoords="offset points", xytext=(6, -12), fontsize=9, color="C3"
    )

No description has been provided for this image

The shaded hexagon is every net-position combination the grid can host; each grey line is one CNEC and the six touching the hexagon are its edges. The region is not a box: A's room to export depends on what B does.

Case 1 is strictly interior and no CNEC is active, so the market could move in any direction from here:

In [6]:

Copied!

check(case1)
check(case1)

Out[6]:

	loading	RAM	margin	active
cnec
AB (+)	500.0	1000.0	500.0	False
AB (-)	-500.0	1000.0	1500.0	False
AC (+)	500.0	2000.0	1500.0	False
AC (-)	-500.0	2000.0	2500.0	False
BC (+)	-0.0	1500.0	1500.0	False
BC (-)	0.0	1500.0	1500.0	False

Case 2 loads the direct line AB hardest, so AB (+) is the single binding CNEC (zero margin):

In [7]:

Copied!

check(case2)
check(case2)

Out[7]:

	loading	RAM	margin	active
cnec
AB (+)	1000.0	1000.0	0.0	True
AB (-)	-1000.0	1000.0	2000.0	False
AC (+)	500.0	2000.0	1500.0	False
AC (-)	-500.0	2000.0	2500.0	False
BC (+)	-500.0	1500.0	2000.0	False
BC (-)	500.0	1500.0	1000.0	False

Case 3 binds AB (+) and AC (+) at once and the two constraints pin the outcome:

In [8]:

Copied!

check(case3)
check(case3)

Out[8]:

	loading	RAM	margin	active
cnec
AB (+)	1000.0	1000.0	0.0	True
AB (-)	-1000.0	1000.0	2000.0	False
AC (+)	2000.0	2000.0	0.0	True
AC (-)	-2000.0	2000.0	4000.0	False
BC (+)	1000.0	1500.0	500.0	False
BC (-)	-1000.0	1500.0	2500.0	False

Part 2 - Clearing the market¶

So far we probed various net positions by hand. Now we give each zone generators and load and let PyPSA optimize the net position that meets demand at least cost while staying inside the domain. This is market clearing.

Representing FBMC in PyPSA¶

A zone's net position is simply its generation minus its load. We need a device that lets each zone be individually imbalanced while the system balances overall, and that exposes each $NP_z$ as a variable to constrain. We use an auxiliary pool bus:

one bus per zone (A, B, C) carrying that zone's generators and load;
one extra auxiliary bus pool as a modelling device;
one link zone -> pool per zone, where its flow denotes that zone's net position;
the pool's own energy balance then forces $\sum_z NP_z = 0$ automatically;
the flow-based domain is one custom constraint on those link flows.

No physical lines or links between zones must be added: the grid physics lives entirely in the PTDF/RAM domain. For the example, we choose costs to create tension: zone A is cheap (10 EUR/MWh) and wants to export; B is expensive (80) with the largest load; C is in between (50).

We add the custom constraint by taking the linopy model from n.optimize.create_model(), adding the flow-based constraints, and solving with n.optimize.solve_model().

In [9]:

Copied!





loads = pd.Series({"A": 500.0, "B": 1500.0, "C": 1000.0})  # MW
cost = pd.Series({"A": 10.0, "B": 80.0, "C": 50.0})  # EUR/MWh


def add_fb_domain(m: linopy.Model, zones: list[str], ram: pd.Series) -> None:
    """Add PTDF . NP <= RAM to model m over the given domain columns."""
    netpos = m["Link-p"].sel(name=zones)
    m.add_constraints(netpos @ domain[zones] <= ram, name="fb_domain")


n = pypsa.Network()
n.add("Bus", [*zones, "pool"])
n.add("Load", zones, bus=zones, p_set=loads)
n.add("Generator", zones, bus=zones, p_nom=4000, marginal_cost=cost)
n.add("Link", zones, bus0=zones, bus1="pool", p_nom=10000, p_min_pu=-1)  # flow = NP_z
n.sanitize()

m = n.optimize.create_model()
add_fb_domain(m, zones, domain["RAM"])
n.optimize.solve_model(log_to_console=False)
loads = pd.Series({"A": 500.0, "B": 1500.0, "C": 1000.0})  # MW
cost = pd.Series({"A": 10.0, "B": 80.0, "C": 50.0})  # EUR/MWh


def add_fb_domain(m: linopy.Model, zones: list[str], ram: pd.Series) -> None:
    """Add PTDF . NP <= RAM to model m over the given domain columns."""
    netpos = m["Link-p"].sel(name=zones)
    m.add_constraints(netpos @ domain[zones] <= ram, name="fb_domain")


n = pypsa.Network()
n.add("Bus", [*zones, "pool"])
n.add("Load", zones, bus=zones, p_set=loads)
n.add("Generator", zones, bus=zones, p_nom=4000, marginal_cost=cost)
n.add("Link", zones, bus0=zones, bus1="pool", p_nom=10000, p_min_pu=-1)  # flow = NP_z
n.sanitize()

m = n.optimize.create_model()
add_fb_domain(m, zones, domain["RAM"])
n.optimize.solve_model(log_to_console=False)

/home/runner/work/PyPSA/PyPSA/pypsa/network/io.py:2079: FutureWarning: pandas infers the `str` dtype for string data since its version 3.0. PyPSA still converts it back to numpy object dtype on import, but will keep it from PyPSA 2.0 on. Set `pypsa.options.api.legacy_string_dtype` explicitly to suppress this warning.
  new_static = _coerce_string_dtypes(new_static)
INFO:pypsa.consistency:Sanitizing network...

INFO:pypsa.components._types.carriers:Adding 1 missing carriers: ['AC']

INFO:pypsa.components._types.carriers:Assigned colors to 1 carriers using 'tab10' palette.

INFO:pypsa.consistency:Network sanitization complete.

/opt/hostedtoolcache/Python/3.13.16/x64/lib/python3.13/contextlib.py:85: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
  return func(*args, **kwds)

INFO:linopy.model: Solve problem using Gurobi solver

INFO:linopy.model:Solver options:
 - log_to_console: False

INFO:linopy.io: Writing time: 0.03s

Set parameter WLSAccessID

Set parameter WLSSecret

Set parameter LicenseID to value 2537914

Academic license 2537914 - for non-commercial use only - registered to l.___@tu-berlin.de

Read LP format model from file /tmp/linopy-problem-tizteex_.lp

Reading time = 0.00 seconds

obj: 22 rows, 6 columns, 33 nonzeros

Set parameter LogToConsole to value 0

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 6 primals, 22 duals
Objective: 6.50e+04
Solver: gurobi
Runtime: 0.00s
Dual bound: 6.50e+04
Solver model: available
Solver message: 2

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, fb_domain were not assigned to the network.

Out[9]:

('ok', 'optimal')

Net positions and generation¶

The link flows are the net positions; generation shows who supplies whom.

In [10]:

Copied!





pd.DataFrame(
    {
        "net position": n.links_t.p0.iloc[0],
        "generation": n.generators_t.p.iloc[0],
        "load": pd.Series(loads),
    }
).round(1)
pd.DataFrame(
    {
        "net position": n.links_t.p0.iloc[0],
        "generation": n.generators_t.p.iloc[0],
        "load": pd.Series(loads),
    }
).round(1)

Out[10]:

	net position	generation	load
A	2000.0	2500.0	500.0
B	-1000.0	500.0	1500.0
C	-1000.0	0.0	1000.0

A exports 2000 MW; B and C import 1000 MW each. But A does not serve everyone and B still runs its expensive 500 MW. A CNEC capped A's export.

From net positions to physical flows¶

The net position is a commercial quantity. Its physical consequence is the flow it puts on the grid: $F_c = \sum_z \mathrm{PTDF}_{c,z}\, NP_z$, exactly the left-hand side of each domain inequality. In this toy model, every CNEC is a line between two regions, so these loadings are the actual MW flowing between the regions. Collapsing each line's two directions into one signed flow yields:

In [11]:

Copied!





netpos = n.links_t.p0.iloc[0]
flows = pd.DataFrame(
    {
        "physical flow": (pos[zones] @ netpos).round(1),  # signed MW on each line
        "Fmax": pos["RAM"],
    }
)
flows
netpos = n.links_t.p0.iloc[0]
flows = pd.DataFrame(
    {
        "physical flow": (pos[zones] @ netpos).round(1),  # signed MW on each line
        "Fmax": pos["RAM"],
    }
)
flows

Out[11]:

	physical flow	Fmax
AB	1000.0	1000.0
BC	0.0	1500.0
AC	1000.0	2000.0

Line AB carries 1000 MW (its limit) while AC carries 1000 MW of its 2000 MW and BC carries nothing. The single line AB decides the whole outcome; on the domain picture the result lands exactly on the AB (+) edge:

In [12]:

Copied!





fig, ax = plt.subplots(figsize=(5, 5), layout="constrained")
draw_domain(ax)
ax.scatter(netpos["A"], netpos["B"], color="C3", zorder=3, s=40)
ax.annotate(
    "optimum",
    (netpos["A"], netpos["B"]),
    textcoords="offset points",
    xytext=(8, -4),
    fontsize=9,
    color="C3",
);
fig, ax = plt.subplots(figsize=(5, 5), layout="constrained")
draw_domain(ax)
ax.scatter(netpos["A"], netpos["B"], color="C3", zorder=3, s=40)
ax.annotate(
    "optimum",
    (netpos["A"], netpos["B"]),
    textcoords="offset points",
    xytext=(8, -4),
    fontsize=9,
    color="C3",
);

Prices and shadow prices¶

Two sets of dual variables are interesting to review:

zonal prices $\pi_z$ - the dual of each zone's energy balance (EUR/MWh);
CNEC shadow prices $\mu_c$ - the dual of each domain inequality, which describes the welfare gained per extra MW of RAM on that CNEC. A binding CNEC is also called an active constraint. linopy reports the dual of a $\le$ constraint as $\le 0$, so $|\mu_c|$ is that welfare value in EUR/MW.

In [13]:

Copied!

prices = n.buses_t.marginal_price.iloc[0].round(2)
prices.to_frame("price (EUR/MWh)")
prices = n.buses_t.marginal_price.iloc[0].round(2)
prices.to_frame("price (EUR/MWh)")

Out[13]:

	price (EUR/MWh)
name
A	10.0
B	80.0
C	45.0
pool	45.0

In [14]:

Copied!

shadow = n.model.constraints["fb_domain"].dual.sel(snapshot="now").to_pandas().round(2)
shadow[shadow != 0].to_frame("shadow price (EUR/MW)")
shadow = n.model.constraints["fb_domain"].dual.sel(snapshot="now").to_pandas().round(2)
shadow[shadow != 0].to_frame("shadow price (EUR/MW)")

Out[14]:

	shadow price (EUR/MW)
cnec
AB (+)	-105.0

One binding CNEC (AB (+), $|\mu| = 105$) splits the three zones into three prices: cheap A stays at 10, blocked B pays its own 80, C settles at 45.

The base price. The pool bus has no generation or load, so its price $\lambda$ is the multiplier on the balance $\sum_z NP_z = 0$: the copper-plate reference price that every zone would pay if no CNEC bound. Congestion then spreads the zonal prices around $\lambda$. That means zone $z$ sits at $\lambda$ plus its PTDF-weighted shadow prices:

$$\pi_z = \lambda + \sum_c \mu_c\, \mathrm{PTDF}_{c,z}.$$

We can do a plausibility check that must reproduce the solver's zonal prices:

In [15]:

Copied!





lam = prices["pool"]
mu = n.model.constraints["fb_domain"].dual.sel(snapshot="now").to_pandas()
pd.DataFrame(
    {"solver": prices[zones], "identity": (lam + domain[zones].T @ mu).round(2)}
)
lam = prices["pool"]
mu = n.model.constraints["fb_domain"].dual.sel(snapshot="now").to_pandas()
pd.DataFrame(
    {"solver": prices[zones], "identity": (lam + domain[zones].T @ mu).round(2)}
)

Out[15]:

	solver	identity
A	10.0	10.0
B	80.0	80.0
C	45.0	45.0

They match. The spread between A and B is then $\sum_c \mu_c (\mathrm{PTDF}_{c,A} - \mathrm{PTDF}_{c,B}) = (-105)\,(\tfrac13 - (-\tfrac13)) = -70$: B is 70 EUR/MWh higher than A, entirely attributable to the AB (+) constraint. This is how a grid limit becomes a price signal.

Part 3 - External borders¶

Real markets have borders to zones outside the flow-based region via HVDC links and neighbouring AC systems. Their exchanges still load the domain's CNECs, and there are two ways to handle them.

Standard Hybrid Coupling (SHC) - the external exchange is a fixed forecast. The flow it induces on the CNECs is reserved out of RAM before the market optimization runs, through the $F_{uaf}$ term of the RAM formula. The market cannot change it as it is treated as a static capacity reservation.

Advanced Hybrid Coupling (AHC) - the border becomes a virtual hub and receives its own PTDF column in the domain and a net position optimised by the market, against the RAM before that reservation. The reserved capacity moves from $F_{uaf}$ back into usable RAM.

The RAM formula¶

The margin published for each CNEC is built up from the physical line limit as

$$\mathrm{RAM} = F_{max} - \mathrm{FRM} - F_{uaf} - F_{0} - F_{LTN} + \mathrm{AMR} - \mathrm{IVA},$$

where $F_{max}$ is the maximum physical flow; $\mathrm{FRM}$ (flow reliability margin) a safety buffer (e.g. 10%); $F_{uaf}$ the flow from exchanges with non-flow-based borders (the forecast handled here); $F_{0}$ the internal flow at all-zero net positions; $F_{LTN}$ long-term nominations (low relevance); $\mathrm{AMR}$ an added minimum RAM (regulatory virtual capacity, e.g. 70% of $F_{max}$); and $\mathrm{IVA}$ a TSO validation reduction. Only $F_{uaf}$ concerns us; the others we fold into a single physical margin we take as RAM.

We add one external market X: a cheap import (20 EUR/MWh) arriving via an HVDC that lands at node B, so its PTDF column equals B's. To compare the two methods on equal footing we start from that physical margin RAM and, for SHC, subtract $F_{uaf} = \mathrm{PTDF}_{c,X}\, NP_X^{\mathrm{forecast}}$ ourselves. Under SHC the forecast is assumed to be a conservative 300 MW.

In [16]:

Copied!





domain["X"] = domain["B"]  # X lands at B -> same PTDF column
cost_x, forecast = 20.0, 300.0  # X price (EUR/MWh), SHC forecast export (MW)


def solve(mode: str) -> pypsa.Network:
    """Clear the market with the external border treated as SHC or AHC."""
    ext = [*zones, "X"]
    cols = ext if mode == "ahc" else zones  # domain columns
    ram = domain["RAM"] - (domain["X"] * forecast if mode == "shc" else 0.0)

    n = pypsa.Network()
    n.add("Bus", [*ext, "pool"])
    n.add("Load", zones, bus=zones, p_set=loads)
    n.add(
        "Generator",
        ext,
        bus=ext,
        p_nom=6000,
        marginal_cost=[cost.get(z, cost_x) for z in ext],
    )
    n.add("Link", ext, bus0=ext, bus1="pool", p_nom=10000, p_min_pu=-1)
    if mode == "shc":  # freeze X at the forecast
        n.links.loc["X", "p_set"] = forecast
    n.sanitize()

    m = n.optimize.create_model()
    add_fb_domain(m, cols, ram)
    n.optimize.solve_model(log_to_console=False)
    return n
domain["X"] = domain["B"]  # X lands at B -> same PTDF column
cost_x, forecast = 20.0, 300.0  # X price (EUR/MWh), SHC forecast export (MW)


def solve(mode: str) -> pypsa.Network:
    """Clear the market with the external border treated as SHC or AHC."""
    ext = [*zones, "X"]
    cols = ext if mode == "ahc" else zones  # domain columns
    ram = domain["RAM"] - (domain["X"] * forecast if mode == "shc" else 0.0)

    n = pypsa.Network()
    n.add("Bus", [*ext, "pool"])
    n.add("Load", zones, bus=zones, p_set=loads)
    n.add(
        "Generator",
        ext,
        bus=ext,
        p_nom=6000,
        marginal_cost=[cost.get(z, cost_x) for z in ext],
    )
    n.add("Link", ext, bus0=ext, bus1="pool", p_nom=10000, p_min_pu=-1)
    if mode == "shc":  # freeze X at the forecast
        n.links.loc["X", "p_set"] = forecast
    n.sanitize()

    m = n.optimize.create_model()
    add_fb_domain(m, cols, ram)
    n.optimize.solve_model(log_to_console=False)
    return n

The whole SHC/AHC difference sits in two lines of solve: SHC subtracts $F_{uaf}$ from RAM and freezes X's link; AHC keeps X's column in the domain and lets its link optimise. Clear both and compare.

In [17]:

Copied!





def summarise(n: pypsa.Network) -> pd.Series:
    netpos = n.links_t.p0.iloc[0]
    price = n.buses_t.marginal_price.iloc[0]
    dual = n.model.constraints["fb_domain"].dual.sel(snapshot="now").to_pandas()
    return pd.Series(
        {
            **{f"NP {z}": netpos[z] for z in [*zones, "X"]},
            **{f"price {z}": price[z] for z in [*zones, "X"]},
            "shadow AB(+)": -dual.get("AB (+)", 0.0),
            "total cost": n.objective,
        }
    )


nets = {mode: solve(mode) for mode in ["shc", "ahc"]}
compare = pd.DataFrame(
    {mode.upper(): summarise(nn) for mode, nn in nets.items()}
).round(1)
compare
def summarise(n: pypsa.Network) -> pd.Series:
    netpos = n.links_t.p0.iloc[0]
    price = n.buses_t.marginal_price.iloc[0]
    dual = n.model.constraints["fb_domain"].dual.sel(snapshot="now").to_pandas()
    return pd.Series(
        {
            **{f"NP {z}": netpos[z] for z in [*zones, "X"]},
            **{f"price {z}": price[z] for z in [*zones, "X"]},
            "shadow AB(+)": -dual.get("AB (+)", 0.0),
            "total cost": n.objective,
        }
    )


nets = {mode: solve(mode) for mode in ["shc", "ahc"]}
compare = pd.DataFrame(
    {mode.upper(): summarise(nn) for mode, nn in nets.items()}
).round(1)
compare

INFO:pypsa.consistency:Sanitizing network...

INFO:pypsa.components._types.carriers:Adding 1 missing carriers: ['AC']

INFO:pypsa.components._types.carriers:Assigned colors to 1 carriers using 'tab10' palette.

INFO:pypsa.consistency:Network sanitization complete.

/opt/hostedtoolcache/Python/3.13.16/x64/lib/python3.13/contextlib.py:85: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
  return func(*args, **kwds)

INFO:linopy.model: Solve problem using Gurobi solver

INFO:linopy.model:Solver options:
 - log_to_console: False

INFO:linopy.io: Writing time: 0.02s

Set parameter WLSAccessID

Set parameter WLSSecret

Set parameter LicenseID to value 2537914

Academic license 2537914 - for non-commercial use only - registered to l.___@tu-berlin.de

Read LP format model from file /tmp/linopy-problem-8y1o3esz.lp

Reading time = 0.00 seconds

obj: 28 rows, 8 columns, 41 nonzeros

Set parameter LogToConsole to value 0

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 8 primals, 31 duals
Objective: 4.70e+04
Solver: gurobi
Runtime: 0.00s
Dual bound: 4.70e+04
Solver model: available
Solver message: 2

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, Link-p_set, fb_domain were not assigned to the network.

INFO:pypsa.consistency:Sanitizing network...

INFO:pypsa.components._types.carriers:Adding 1 missing carriers: ['AC']

INFO:pypsa.components._types.carriers:Assigned colors to 1 carriers using 'tab10' palette.

INFO:pypsa.consistency:Network sanitization complete.

/opt/hostedtoolcache/Python/3.13.16/x64/lib/python3.13/contextlib.py:85: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
  return func(*args, **kwds)

INFO:linopy.model: Solve problem using Gurobi solver

INFO:linopy.model:Solver options:
 - log_to_console: False

INFO:linopy.io: Writing time: 0.02s

Set parameter WLSAccessID

Set parameter WLSSecret

Set parameter LicenseID to value 2537914

Academic license 2537914 - for non-commercial use only - registered to l.___@tu-berlin.de

Read LP format model from file /tmp/linopy-problem-x3p7guwr.lp

Reading time = 0.00 seconds

obj: 27 rows, 8 columns, 46 nonzeros

Set parameter LogToConsole to value 0

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 8 primals, 27 duals
Objective: 3.50e+04
Solver: gurobi
Runtime: 0.00s
Dual bound: 3.50e+04
Solver model: available
Solver message: 2

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Link-fix-p-lower, Link-fix-p-upper, fb_domain were not assigned to the network.

Out[17]:

	SHC	AHC
NP A	2000.0	2000.0
NP B	-1300.0	-1500.0
NP C	-1000.0	-1000.0
NP X	300.0	500.0
price A	10.0	10.0
price B	80.0	20.0
price C	45.0	15.0
price X	20.0	20.0
shadow AB(+)	105.0	15.0
total cost	47000.0	35000.0

Reading the columns:

X net position: frozen at 300 MW under SHC; the market pulls it to 500 MW under AHC.
Congestion: AB (+) binds in both, but its shadow price collapses from 105 €/MW to 15 €/MW. The extra cheap import relieves the very border that limited trade.
Prices: under SHC blocked zone B pays its own 80 €/MWh; under AHC the import reaches it and B falls to 20 €/MWh, compressing the spread.
Welfare: total cost drops from 47000 € to 35000 €, a 12000 € gain from 200 MW of cheap X displacing expensive B generation.

In [18]:

Copied!





ent = [*zones, "X"]
x = np.arange(len(ent))
fig, ax = plt.subplots(figsize=(6, 4), layout="constrained")
for i, (mode, col) in enumerate(
    zip(["SHC", "AHC"], ["#5ec962", "#3b528b"], strict=True)
):
    bars = ax.bar(
        x + (i - 0.5) * 0.4,
        compare.loc[[f"price {z}" for z in ent], mode],
        width=0.4,
        label=mode,
        color=col,
    )
    ax.bar_label(bars, fmt="%.0f", fontsize=8, padding=2)
ax.set_xticks(x, ent)
ax.set(
    xlabel="zone",
    ylabel="price (EUR/MWh)",
)
ax.legend(frameon=False);
ent = [*zones, "X"]
x = np.arange(len(ent))
fig, ax = plt.subplots(figsize=(6, 4), layout="constrained")
for i, (mode, col) in enumerate(
    zip(["SHC", "AHC"], ["#5ec962", "#3b528b"], strict=True)
):
    bars = ax.bar(
        x + (i - 0.5) * 0.4,
        compare.loc[[f"price {z}" for z in ent], mode],
        width=0.4,
        label=mode,
        color=col,
    )
    ax.bar_label(bars, fmt="%.0f", fontsize=8, padding=2)
ax.set_xticks(x, ent)
ax.set(
    xlabel="zone",
    ylabel="price (EUR/MWh)",
)
ax.legend(frameon=False);

Appendix: FBMC without auxiliary components¶

The pool bus and the zone -> pool links are a modelling device, not part of the market. Their only jobs are to give each zone an import/export outlet and to enforce sum(NP) = 0. Both can be done instead by modifying the underlying linopy model: add one net-position variable per zone directly inside the nodal balance,

$$\text{gen}_z - \text{load}_z - NP_z = 0 ,$$

keep the balance (so the zonal prices still come from its dual), add a global sum(NP) = 0 constraint, and put the flow-based domain on the NP variables. No auxiliary components needed.

In [19]:

Copied!





na = pypsa.Network()
na.add("Bus", zones)
na.add("Load", zones, bus=zones, p_set=loads)
na.add("Generator", zones, bus=zones, p_nom=4000, marginal_cost=cost)
na.sanitize()

m = na.optimize.create_model()
balance = m.constraints["Bus-nodal_balance"]  # dim "name" = the zone buses
netpos = m.add_variables(
    coords=[na.snapshots, pd.Index(zones, name="name")], name="net_position"
)
balance.update(lhs=balance.lhs - netpos)  # gen - load - NP = 0  ->  NP = gen - load
m.add_constraints(netpos.sum("name") == 0, name="global_balance")
m.add_constraints(netpos @ domain[zones] <= domain["RAM"], name="fb_domain")
na.optimize.solve_model()

pd.DataFrame(
    {
        "net position (MW)": (na.generators_t.p.iloc[0] - loads)[zones].round(0),
        "price (EUR/MWh)": na.buses_t.marginal_price.iloc[0][zones].round(1),
    }
)
na = pypsa.Network()
na.add("Bus", zones)
na.add("Load", zones, bus=zones, p_set=loads)
na.add("Generator", zones, bus=zones, p_nom=4000, marginal_cost=cost)
na.sanitize()

m = na.optimize.create_model()
balance = m.constraints["Bus-nodal_balance"]  # dim "name" = the zone buses
netpos = m.add_variables(
    coords=[na.snapshots, pd.Index(zones, name="name")], name="net_position"
)
balance.update(lhs=balance.lhs - netpos)  # gen - load - NP = 0  ->  NP = gen - load
m.add_constraints(netpos.sum("name") == 0, name="global_balance")
m.add_constraints(netpos @ domain[zones] <= domain["RAM"], name="fb_domain")
na.optimize.solve_model()

pd.DataFrame(
    {
        "net position (MW)": (na.generators_t.p.iloc[0] - loads)[zones].round(0),
        "price (EUR/MWh)": na.buses_t.marginal_price.iloc[0][zones].round(1),
    }
)

INFO:pypsa.consistency:Sanitizing network...

INFO:pypsa.components._types.carriers:Adding 1 missing carriers: ['AC']

INFO:pypsa.components._types.carriers:Assigned colors to 1 carriers using 'tab10' palette.

INFO:pypsa.consistency:Network sanitization complete.

/opt/hostedtoolcache/Python/3.13.16/x64/lib/python3.13/contextlib.py:85: FutureWarning: The default value of `include_objective_constant` will change from True to False in version 2.0. Set `include_objective_constant` explicitly to suppress this warning. Using False improves LP numerical conditioning by not including the objective constant as a variable.
  return func(*args, **kwds)

INFO:linopy.model: Solve problem using Gurobi solver

INFO:linopy.model:Solver options:
 - log_to_console: False

INFO:linopy.io: Writing time: 0.02s

Set parameter WLSAccessID

Set parameter WLSSecret

Set parameter LicenseID to value 2537914

Academic license 2537914 - for non-commercial use only - registered to l.___@tu-berlin.de

Read LP format model from file /tmp/linopy-problem-brsbk_cp.lp

Reading time = 0.00 seconds

obj: 16 rows, 6 columns, 27 nonzeros

Set parameter LogToConsole to value 0

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 6 primals, 16 duals
Objective: 6.50e+04
Solver: gurobi
Runtime: 0.00s
Dual bound: 6.50e+04
Solver model: available
Solver message: 2

INFO:pypsa.optimization.optimize:The variable 'net_position' could not be mapped to the network component because it does not include the symbol '-'.

INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, global_balance, fb_domain were not assigned to the network.

Out[19]:

	net position (MW)	price (EUR/MWh)
name
A	2000.0	10.0
B	-1000.0	80.0
C	-1000.0	45.0

Identical to the pool clearing of Part 2 with NP = (2000, -1000, -1000) and prices A=10, B=80, C=45.