Determining Flow-Based Domains¶

The previous example took the flow-based domain (the zonal PTDF matrix and the RAM vector) as given input. That is how a market modeller usually receives it. But someone has to compute it first. This example opens that black box on a small, stylised grid and explains how physical network data can be turned into a zonal PTDF . NP <= RAM flow-based domain.

The recipe, one step per section:

1. Nodal PTDFs - how nodal injections load each line (linearized grid physics).
2. Generation Shift Keys (GSK) - how a zone's net position is shared out over its nodes.
3. Zonal PTDFs - fold the nodal PTDF through the GSK to get one column per zone.
4. CNEC selection - keep only the lines that trade actually stresses.
5. Base case and RAM - subtract reference flows and safety margins from the line limits.
6. Clear the market - feed the result into the flow-based constraints.
7. N-1 security - add contingency constraints that make the domain physically robust.

The stylised grid¶

Three market zones. Zone A contains two nodes (A1, A2); zones B and C are a single node each. Five lines form a meshed grid with two loops, so power splits over parallel paths, i.e. the loop flows that make flow-based coupling necessary. Reactances are equal (x = 1); line ratings s_nom differ so different lines can bind.

In [2]:

import matplotlib.pyplot as plt
import networkx as nx
import numpy as np
import pandas as pd
from scipy.spatial import ConvexHull, HalfspaceIntersection

import pypsa

pypsa.options.params.optimize.log_to_console = False

zones = ["A", "B", "C"]
node_zone = pd.Series({"A1": "A", "A2": "A", "B1": "B", "C1": "C"}, name="zone")
p_nom_node = pd.Series({"A1": 2000.0, "A2": 1000.0, "B1": 4000.0, "C1": 4000.0})

grid = {
    "L_A1A2": ("A1", "A2", 3000.0),  # internal to zone A
    "L_A2B1": ("A2", "B1", 1000.0),  # A-B border
    "L_B1C1": ("B1", "C1", 1500.0),  # B-C border
    "L_C1A1": ("C1", "A1", 2000.0),  # C-A border
    "L_A1B1": ("A1", "B1", 1200.0),  # second A-B path (closes the mesh)
}
n = pypsa.Network()
n.add("Bus", list(node_zone.index))
n.add(
    "Line",
    list(grid),
    bus0=[v[0] for v in grid.values()],
    bus1=[v[1] for v in grid.values()],
    x=1.0,
    r=0.0,
    s_nom=[v[2] for v in grid.values()],
);

import matplotlib.pyplot as plt import networkx as nx import numpy as np import pandas as pd from scipy.spatial import ConvexHull, HalfspaceIntersection import pypsa pypsa.options.params.optimize.log_to_console = False zones = ["A", "B", "C"] node_zone = pd.Series({"A1": "A", "A2": "A", "B1": "B", "C1": "C"}, name="zone") p_nom_node = pd.Series({"A1": 2000.0, "A2": 1000.0, "B1": 4000.0, "C1": 4000.0}) grid = { "L_A1A2": ("A1", "A2", 3000.0), # internal to zone A "L_A2B1": ("A2", "B1", 1000.0), # A-B border "L_B1C1": ("B1", "C1", 1500.0), # B-C border "L_C1A1": ("C1", "A1", 2000.0), # C-A border "L_A1B1": ("A1", "B1", 1200.0), # second A-B path (closes the mesh) } n = pypsa.Network() n.add("Bus", list(node_zone.index)) n.add( "Line", list(grid), bus0=[v[0] for v in grid.values()], bus1=[v[1] for v in grid.values()], x=1.0, r=0.0, s_nom=[v[2] for v in grid.values()], );

/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)

In [3]:

layout = {"A1": (0, 1), "A2": (1.4, 1), "B1": (1.4, 0), "C1": (0, 0)}
zone_colour = {"A": "#0072B2", "B": "#D55E00", "C": "#009E73"}

g = nx.Graph()
g.add_edges_from((v[0], v[1]) for v in grid.values())

fig, ax = plt.subplots(figsize=(4.5, 4.2), layout="constrained")
nx.draw_networkx_edges(g, layout, ax=ax, edge_color="0.6")
nx.draw_networkx_nodes(
    g,
    layout,
    ax=ax,
    node_size=1400,
    node_color=[zone_colour[node_zone[b]] for b in g.nodes],
)
nx.draw_networkx_labels(g, layout, ax=ax, font_color="white", font_weight="bold")
nx.draw_networkx_edge_labels(
    g,
    layout,
    ax=ax,
    font_size=8,
    edge_labels={(v[0], v[1]): f"{v[2]:.0f} MW" for v in grid.values()},
)
ax.set_axis_off();

layout = {"A1": (0, 1), "A2": (1.4, 1), "B1": (1.4, 0), "C1": (0, 0)} zone_colour = {"A": "#0072B2", "B": "#D55E00", "C": "#009E73"} g = nx.Graph() g.add_edges_from((v[0], v[1]) for v in grid.values()) fig, ax = plt.subplots(figsize=(4.5, 4.2), layout="constrained") nx.draw_networkx_edges(g, layout, ax=ax, edge_color="0.6") nx.draw_networkx_nodes( g, layout, ax=ax, node_size=1400, node_color=[zone_colour[node_zone[b]] for b in g.nodes], ) nx.draw_networkx_labels(g, layout, ax=ax, font_color="white", font_weight="bold") nx.draw_networkx_edge_labels( g, layout, ax=ax, font_size=8, edge_labels={(v[0], v[1]): f"{v[2]:.0f} MW" for v in grid.values()}, ) ax.set_axis_off();

Step 1 - Nodal PTDF¶

The nodal PTDF answers a purely physical question: if we inject 1 MW at node $n$ and withdraw it at the slack node, how much extra flow appears on line $l$? It follows from a DC power-flow linearisation of the grid, and PyPSA computes it for us:

$$F_l = \sum_n \text{PTDF}^{\text{nodal}}_{l,n}\, P_n .$$

The slack node's column is zero (its injection is the balancing one). Here PyPSA picks A1 as slack.

In [4]:

n.determine_network_topology()
sub = n.sub_networks.obj.iloc[0]
sub.calculate_PTDF()

H = pd.DataFrame(sub.PTDF, index=n.lines.index, columns=sub.buses_o)
H.index.name, H.columns.name = "line", "node"
H.round(3)

n.determine_network_topology() sub = n.sub_networks.obj.iloc[0] sub.calculate_PTDF() H = pd.DataFrame(sub.PTDF, index=n.lines.index, columns=sub.buses_o) H.index.name, H.columns.name = "line", "node" H.round(3)

Out[4]:

node	A1	A2	B1	C1
line
L_A1A2	0.0	-0.625	-0.25	-0.125
L_A2B1	0.0	0.375	-0.25	-0.125
L_B1C1	0.0	0.125	0.25	-0.375
L_C1A1	0.0	0.125	0.25	0.625
L_A1B1	0.0	-0.250	-0.50	-0.250

Step 2 - Generation Shift Keys (GSK)¶

Zones, not nodes, trade. A zone's net position $NP_z$ (net export) is a single number, but physically it is an injection spread over the zone's nodes. The GSK encodes that split:

$$\text{GSK}_{n,z} = \frac{\partial P_n}{\partial NP_z},\qquad \sum_{n \in z}\text{GSK}_{n,z}=1 .$$

A GSK of 0.3 means node $n$ takes 0.3 MW of every extra MW that zone $z$ exports. Only zone A has a choice to make. We weight its two nodes by installed capacity (A1 twice A2), a common rule; others (flat, or pro rata to the forecast dispatch) would give different weights and a different domain.

In [5]:

gsk = p_nom_node.groupby(node_zone).transform(lambda s: s / s.sum())
gsk

gsk = p_nom_node.groupby(node_zone).transform(lambda s: s / s.sum()) gsk

Out[5]:

A1    0.666667
A2    0.333333
B1    1.000000
C1    1.000000
dtype: float64

Step 3 - Zonal PTDF¶

Fold the nodal PTDF through the GSK to get the sensitivity of each line to a zonal net position:

$$\text{PTDF}^{\text{zonal}}_{l,z} = \sum_n \text{PTDF}^{\text{nodal}}_{l,n}\,\text{GSK}_{n,z} .$$

Weight each nodal column by its GSK, then sum the columns within each zone.

In [6]:

def zonal(Hn: pd.DataFrame) -> pd.DataFrame:
    return Hn.multiply(gsk).T.groupby(node_zone).sum().T


Z = zonal(H)
Z.round(3)

def zonal(Hn: pd.DataFrame) -> pd.DataFrame: return Hn.multiply(gsk).T.groupby(node_zone).sum().T Z = zonal(H) Z.round(3)

Out[6]:

zone	A	B	C
line
L_A1A2	-0.208	-0.25	-0.125
L_A2B1	0.125	-0.25	-0.125
L_B1C1	0.042	0.25	-0.375
L_C1A1	0.042	0.25	0.625
L_A1B1	-0.083	-0.50	-0.250

The B and C columns come straight from those single nodes; zone A's column is the capacity-weighted blend of A1 and A2.

Step 4 - Critical Network Element selection¶

Not every line needs a constraint. A line is a CNEC (critical network element) if trade loads it meaningfully. The rule of thumb is a zonal PTDF above 5% for at least one zone. Note this catches internal lines too: L_A1A2 sits inside zone A yet is strongly loaded by A's trades, so it is critical.

In [7]:

thresh = 0.05
sens = Z.abs().max(axis=1)
cnecs = Z.index[sens > thresh]
pd.DataFrame({"max |PTDF|": sens.round(3), "is CNEC": sens > thresh})

thresh = 0.05 sens = Z.abs().max(axis=1) cnecs = Z.index[sens > thresh] pd.DataFrame({"max |PTDF|": sens.round(3), "is CNEC": sens > thresh})

Out[7]:

	max |PTDF|	is CNEC
line
L_A1A2	0.250	True
L_A2B1	0.250	True
L_B1C1	0.375	True
L_C1A1	0.625	True
L_A1B1	0.500	True

In this small mesh all five lines qualify, the internal one included. On a real grid the 5% filter discards the many lines that trade barely touches, which is what keeps the domain small.

Step 5 - Base case and Remaining Available Margin¶

A line's thermal rating $F^{\max}$ is not all available to the day-ahead market. Some bites are taken out of it, e.g.:

• the base-case flow $F_0$ - the loop/reference flow already on the line from trades settled elsewhere (forward, bilateral, internal). We get it by running a reference nodal injection $P^{bc}$ through the nodal PTDF, $F_0 = \text{PTDF}^{\text{nodal}} P^{bc}$;
• the Flow Reliability Margin $F^{RM}$ - a safety buffer for the approximations. The GSK assumes a forecast of how each zone dispatches; when the real dispatch differs, the true flow differs from the zonal-PTDF estimate, and this margin absorbs that gap. Here, a flat 10% of $F^{\max}$.

$$\text{RAM} = F^{\max} - F^{RM} - F_0 .$$

Because $F_0$ has a direction, the positive and negative margins differ.

In [8]:

p_bc = pd.Series({"A1": 300.0, "A2": 0.0, "B1": -500.0, "C1": 200.0})  # sums to 0
F0 = pd.Series(H.values @ p_bc.reindex(H.columns).values, index=H.index)

Fmax = n.lines.s_nom
frm = 0.10 * Fmax
ram_pos = Fmax - frm - F0  # positive flow direction
ram_neg = Fmax - frm + F0  # negative flow direction

pd.DataFrame(
    {
        "Fmax": Fmax,
        "FRM": frm,
        "F0": F0.round(0),
        "RAM+": ram_pos.round(0),
        "RAM-": ram_neg.round(0),
    }
).loc[cnecs]

p_bc = pd.Series({"A1": 300.0, "A2": 0.0, "B1": -500.0, "C1": 200.0}) # sums to 0 F0 = pd.Series(H.values @ p_bc.reindex(H.columns).values, index=H.index) Fmax = n.lines.s_nom frm = 0.10 * Fmax ram_pos = Fmax - frm - F0 # positive flow direction ram_neg = Fmax - frm + F0 # negative flow direction pd.DataFrame( { "Fmax": Fmax, "FRM": frm, "F0": F0.round(0), "RAM+": ram_pos.round(0), "RAM-": ram_neg.round(0), } ).loc[cnecs]

Out[8]:

	Fmax	FRM	F0	RAM+	RAM-
line
L_A1A2	3000.0	300.0	100.0	2600.0	2800.0
L_A2B1	1000.0	100.0	100.0	800.0	1000.0
L_B1C1	1500.0	150.0	-200.0	1550.0	1150.0
L_C1A1	2000.0	200.0	0.0	1800.0	1800.0
L_A1B1	1200.0	120.0	200.0	880.0	1280.0

Step 6 - Assemble the domain and clear the market¶

Each CNEC becomes two signed rows (+ and - direction), giving the PTDF . NP <= RAM domain.

In [9]:

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


def domain_of(
    Zsel: pd.DataFrame, ram_pos: pd.Series, ram_neg: pd.Series
) -> pd.DataFrame:
    """Signed +/- rows PTDF . NP <= RAM for the given CNECs."""
    pos = Zsel.copy()
    pos["RAM"] = ram_pos.values
    neg = -Zsel.copy()
    neg["RAM"] = ram_neg.values
    d = pd.concat([pos.rename(lambda s: f"{s} (+)"), neg.rename(lambda s: f"{s} (-)")])
    d.index.name = "cnec"
    return d


def clear(domain: pd.DataFrame) -> tuple[pd.Series, pd.Series, pd.Series]:
    """Clear the toy market subject to the flow-based domain; return NP, prices, shadow prices."""
    mk = pypsa.Network()
    mk.add("Bus", [*zones, "pool"])
    mk.add("Load", zones, bus=zones, p_set=loads.values)
    mk.add("Generator", zones, bus=zones, p_nom=4000, marginal_cost=cost.values)
    mk.add("Link", zones, bus0=zones, bus1="pool", p_nom=1e5, p_min_pu=-1)
    mk.sanitize()
    m = mk.optimize.create_model()
    ptdf = domain[zones].rename_axis(columns="name")
    m.add_constraints(
        m["Link-p"].sel(name=zones) @ ptdf <= domain["RAM"], name="fb_domain"
    )
    mk.optimize.solve_model()
    NP = mk.links_t.p0.iloc[0][zones]
    prices = mk.buses_t.marginal_price.iloc[0][zones]
    mu = m.constraints["fb_domain"].dual.sel(snapshot="now").to_pandas()
    return NP, prices, mu[mu.abs() > 1e-4]


domain_N = domain_of(Z.loc[cnecs], ram_pos[cnecs], ram_neg[cnecs])
NP_N, price_N, mu_N = clear(domain_N)
pd.DataFrame({"net position (MW)": NP_N.round(0), "price (EUR/MWh)": price_N.round(1)})

loads = pd.Series({"A": 500.0, "B": 1500.0, "C": 1000.0}) cost = pd.Series({"A": 10.0, "B": 80.0, "C": 50.0}) def domain_of( Zsel: pd.DataFrame, ram_pos: pd.Series, ram_neg: pd.Series ) -> pd.DataFrame: """Signed +/- rows PTDF . NP <= RAM for the given CNECs.""" pos = Zsel.copy() pos["RAM"] = ram_pos.values neg = -Zsel.copy() neg["RAM"] = ram_neg.values d = pd.concat([pos.rename(lambda s: f"{s} (+)"), neg.rename(lambda s: f"{s} (-)")]) d.index.name = "cnec" return d def clear(domain: pd.DataFrame) -> tuple[pd.Series, pd.Series, pd.Series]: """Clear the toy market subject to the flow-based domain; return NP, prices, shadow prices.""" mk = pypsa.Network() mk.add("Bus", [*zones, "pool"]) mk.add("Load", zones, bus=zones, p_set=loads.values) mk.add("Generator", zones, bus=zones, p_nom=4000, marginal_cost=cost.values) mk.add("Link", zones, bus0=zones, bus1="pool", p_nom=1e5, p_min_pu=-1) mk.sanitize() m = mk.optimize.create_model() ptdf = domain[zones].rename_axis(columns="name") m.add_constraints( m["Link-p"].sel(name=zones) @ ptdf <= domain["RAM"], name="fb_domain" ) mk.optimize.solve_model() NP = mk.links_t.p0.iloc[0][zones] prices = mk.buses_t.marginal_price.iloc[0][zones] mu = m.constraints["fb_domain"].dual.sel(snapshot="now").to_pandas() return NP, prices, mu[mu.abs() > 1e-4] domain_N = domain_of(Z.loc[cnecs], ram_pos[cnecs], ram_neg[cnecs]) NP_N, price_N, mu_N = clear(domain_N) pd.DataFrame({"net position (MW)": NP_N.round(0), "price (EUR/MWh)": price_N.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.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-qwxj0jm9.lp

Reading time = 0.00 seconds

obj: 26 rows, 6 columns, 51 nonzeros

Set parameter LogToConsole to value 0

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 6 primals, 26 duals
Objective: 3.20e+04
Solver: gurobi
Runtime: 0.00s
Dual bound: 3.20e+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]:

	net position (MW)	price (EUR/MWh)
name
A	2450.0	10.0
B	-1500.0	70.0
C	-950.0	50.0

In [10]:

mu_N.round(1).rename("shadow price (EUR/MW)").to_frame()

mu_N.round(1).rename("shadow price (EUR/MW)").to_frame()

Out[10]:

	shadow price (EUR/MW)
cnec
L_A2B1 (+)	-160.0

A exports its cheap power and B imports, but the A-B line L_A2B1 (+) binds first (RAM+ = 800 MW). Its shadow price sets the price spread.

Step 7 - N-1 security (critical branches under critical outages)¶

So far every constraint assumed the full grid. But the market outcome must also survive the loss of any single element. For each outage of a line $k$, the flow redistributes onto the surviving lines by the line outage distribution factors (LODF), which PyPSA computes as calculate_BODF. The post-outage nodal PTDF of a monitored line $l$ is

$$\text{PTDF}^{(k)}_{l,n} = \text{PTDF}_{l,n} + \text{LODF}_{l,k}\,\text{PTDF}_{k,n},$$

and folding it through the GSK gives a post-contingency zonal PTDF. Each surviving line under each outage is an extra CNEC - a critical branch under a critical outage (CBCO) with its own RAM. We sweep every single outage and add these rows to the domain.

In [11]:

sub.calculate_BODF()
BODF = pd.DataFrame(sub.BODF, index=n.lines.index, columns=n.lines.index)

rows = []
for k in n.lines.index:  # outaged line
    Hk = H.add(np.outer(BODF[k], H.loc[k]), axis=0)  # post-outage nodal PTDF
    Zk = zonal(Hk)
    F0k = pd.Series(Hk.values @ p_bc.reindex(H.columns).values, index=H.index)
    for l in n.lines.index:  # monitored line
        if l == k or Zk.loc[l].abs().max() <= thresh:
            continue
        rows.append(
            [
                f"{l} | out {k}",
                *Zk.loc[l, zones],
                Fmax[l] - frm[l] - F0k[l],
                Fmax[l] - frm[l] + F0k[l],
            ]
        )

cbco = pd.DataFrame(rows, columns=["cnec", *zones, "RAM+", "RAM-"]).set_index("cnec")
f"{len(cbco)} critical branch-outage pairs (CBCOs)"

sub.calculate_BODF() BODF = pd.DataFrame(sub.BODF, index=n.lines.index, columns=n.lines.index) rows = [] for k in n.lines.index: # outaged line Hk = H.add(np.outer(BODF[k], H.loc[k]), axis=0) # post-outage nodal PTDF Zk = zonal(Hk) F0k = pd.Series(Hk.values @ p_bc.reindex(H.columns).values, index=H.index) for l in n.lines.index: # monitored line if l == k or Zk.loc[l].abs().max() <= thresh: continue rows.append( [ f"{l} | out {k}", *Zk.loc[l, zones], Fmax[l] - frm[l] - F0k[l], Fmax[l] - frm[l] + F0k[l], ] ) cbco = pd.DataFrame(rows, columns=["cnec", *zones, "RAM+", "RAM-"]).set_index("cnec") f"{len(cbco)} critical branch-outage pairs (CBCOs)"

Out[11]:

'20 critical branch-outage pairs (CBCOs)'

In [12]:

domain_cbco = domain_of(cbco[zones], cbco["RAM+"], cbco["RAM-"])
domain_N1 = pd.concat([domain_N, domain_cbco])
NP_N1, price_N1, mu_N1 = clear(domain_N1)

pd.concat(
    {
        "N only": pd.DataFrame({"NP (MW)": NP_N.round(0), "price": price_N.round(1)}),
        "N and N-1": pd.DataFrame(
            {"NP (MW)": NP_N1.round(0), "price": price_N1.round(1)}
        ),
    },
    axis=1,
)

domain_cbco = domain_of(cbco[zones], cbco["RAM+"], cbco["RAM-"]) domain_N1 = pd.concat([domain_N, domain_cbco]) NP_N1, price_N1, mu_N1 = clear(domain_N1) pd.concat( { "N only": pd.DataFrame({"NP (MW)": NP_N.round(0), "price": price_N.round(1)}), "N and N-1": pd.DataFrame( {"NP (MW)": NP_N1.round(0), "price": price_N1.round(1)} ), }, axis=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-ea7evk7o.lp

Reading time = 0.00 seconds

obj: 66 rows, 6 columns, 143 nonzeros

Set parameter LogToConsole to value 0

INFO:linopy.constants: Optimization successful: 
Status: ok
Termination condition: optimal
Solution: 6 primals, 66 duals
Objective: 9.10e+04
Solver: gurobi
Runtime: 0.00s
Dual bound: 9.10e+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[12]:

	N only		N and N-1
	NP (MW)	price	NP (MW)	price
name
A	2450.0	10.0	1120.0	10.0
B	-1500.0	70.0	-1307.0	80.0
C	-950.0	50.0	187.0	50.0

In [13]:

mu_N1.round(1).rename("shadow price (EUR/MW)").to_frame()

mu_N1.round(1).rename("shadow price (EUR/MW)").to_frame()

Out[13]:

	shadow price (EUR/MW)
cnec
L_A2B1 | out L_A1B1 (+)	-120.0

Adding the contingencies tightens the domain: A's export falls from ~2450 to ~1120 MW, and the binding element is now L_A2B1 | out L_A1B1, a critical branch under a critical outage. The market is curtailed so that trade stays secure even if the meshed second A-B path trips.

Step 8 - Domain Plot¶

With three zones and sum(NP) = 0 the domain is two-dimensional: fix NP_C = -(NP_A + NP_B) and every CNEC PTDF . NP <= RAM becomes a straight cut in the (NP_A, NP_B) plane. Their intersection is the feasible polygon. The grey lines are the individual N-state CNEC cuts (one per PTDF row); the blue polygon is their intersection, and the cuts that never touch it are simply slack.

Overlaying the N and N-1 domains shows the effect of including contingency constraints (their cuts omitted for clarity). They shave corners off the N polygon, and the cleared operating point moves from the edge of the large domain onto the boundary of the smaller one.

In [14]:

def polygon(domain: pd.DataFrame) -> np.ndarray:
    """Ordered vertices of {PTDF . NP <= RAM} in the (NP_A, NP_B) plane, NP_C = -(NP_A + NP_B)."""
    a = (domain["A"] - domain["C"]).values
    b = (domain["B"] - domain["C"]).values
    halfspaces = np.column_stack([a, b, -domain["RAM"].values])
    v = HalfspaceIntersection(halfspaces, np.zeros(2)).intersections
    return v[ConvexHull(v).vertices]


vN, vN1 = polygon(domain_N), polygon(domain_N1)

fig, ax = plt.subplots(figsize=(5.5, 5.5), layout="constrained")
span = np.abs(vN).max() * 2
t = np.array([-span, span])
for (
    c1,
    c2,
    r,
) in zip(  # each N-state CNEC = one grey cut; the polygon is their intersection
    (domain_N["A"] - domain_N["C"]).values,
    (domain_N["B"] - domain_N["C"]).values,
    domain_N["RAM"].values,
    strict=True,
):
    xy = (t, (r - c1 * t) / c2) if abs(c2) > abs(c1) else ((r - c2 * t) / c1, t)
    ax.plot(*xy, color="0.8", lw=0.6, zorder=0)
for verts, colour, label in [
    (vN, "#0072B2", "N domain"),
    (vN1, "#D55E00", "N-1 domain"),
]:
    ax.fill(
        verts[:, 0],
        verts[:, 1],
        facecolor=colour,
        edgecolor=colour,
        alpha=0.25,
        lw=1.5,
        label=label,
    )
for netpos, colour, label, off, ha in [
    (NP_N, "#0072B2", "N optimum", (10, 0), "left"),
    (NP_N1, "#D55E00", "N-1 optimum", (-10, 0), "right"),
]:
    ax.plot(netpos["A"], netpos["B"], "o", color=colour, ms=8)
    ax.annotate(
        label,
        (netpos["A"], netpos["B"]),
        textcoords="offset points",
        xytext=off,
        ha=ha,
        va="center",
    )
ax.axhline(0, color="k", linestyle="--", lw=0.6)
ax.axvline(0, color="k", linestyle="--", lw=0.6)
ax.set_aspect("equal")
ax.set_xlim(vN[:, 0].min() - 500, vN[:, 0].max() + 1700)
ax.set_ylim(vN[:, 1].min() - 500, vN[:, 1].max() + 500)
ax.set_xlabel("net position A (MW)")
ax.set_ylabel("net position B (MW)")
ax.legend(loc="lower left", frameon=False)
ax.set_title("N-1 contracts the flow-based domain");

def polygon(domain: pd.DataFrame) -> np.ndarray: """Ordered vertices of {PTDF . NP <= RAM} in the (NP_A, NP_B) plane, NP_C = -(NP_A + NP_B).""" a = (domain["A"] - domain["C"]).values b = (domain["B"] - domain["C"]).values halfspaces = np.column_stack([a, b, -domain["RAM"].values]) v = HalfspaceIntersection(halfspaces, np.zeros(2)).intersections return v[ConvexHull(v).vertices] vN, vN1 = polygon(domain_N), polygon(domain_N1) fig, ax = plt.subplots(figsize=(5.5, 5.5), layout="constrained") span = np.abs(vN).max() * 2 t = np.array([-span, span]) for ( c1, c2, r, ) in zip( # each N-state CNEC = one grey cut; the polygon is their intersection (domain_N["A"] - domain_N["C"]).values, (domain_N["B"] - domain_N["C"]).values, domain_N["RAM"].values, strict=True, ): xy = (t, (r - c1 * t) / c2) if abs(c2) > abs(c1) else ((r - c2 * t) / c1, t) ax.plot(*xy, color="0.8", lw=0.6, zorder=0) for verts, colour, label in [ (vN, "#0072B2", "N domain"), (vN1, "#D55E00", "N-1 domain"), ]: ax.fill( verts[:, 0], verts[:, 1], facecolor=colour, edgecolor=colour, alpha=0.25, lw=1.5, label=label, ) for netpos, colour, label, off, ha in [ (NP_N, "#0072B2", "N optimum", (10, 0), "left"), (NP_N1, "#D55E00", "N-1 optimum", (-10, 0), "right"), ]: ax.plot(netpos["A"], netpos["B"], "o", color=colour, ms=8) ax.annotate( label, (netpos["A"], netpos["B"]), textcoords="offset points", xytext=off, ha=ha, va="center", ) ax.axhline(0, color="k", linestyle="--", lw=0.6) ax.axvline(0, color="k", linestyle="--", lw=0.6) ax.set_aspect("equal") ax.set_xlim(vN[:, 0].min() - 500, vN[:, 0].max() + 1700) ax.set_ylim(vN[:, 1].min() - 500, vN[:, 1].max() + 500) ax.set_xlabel("net position A (MW)") ax.set_ylabel("net position B (MW)") ax.legend(loc="lower left", frameon=False) ax.set_title("N-1 contracts the flow-based domain");

References¶

• Van den Bergh, Boury & Delarue (2016)
• Schönheit et al. (2020/2021)