Line Loading Limits¶
The thermal rating s_nom is not the only reason to cap the flow on a transmission line. Operators derate lines to keep a security margin, adapt ratings to the weather, and restrict voltage angle differences to preserve stability. This example shows how each of these limits is expressed in a PyPSA linear optimal power flow (LOPF) on a small meshed network, and what changes once line capacities become an investment decision.
It covers:
- Static N-1 approximation via
Line.s_max_pu - Dynamic line rating via a time-varying
n.lines_t.s_max_pu - Voltage angle difference limits, first as a manual
s_max_puconversion and then natively viaLine.v_ang_max - Angle limits with extendable capacities, where an absolute angle cap does not scale with
s_nom - Temporary admissible loading (TATL) by relaxing
line2'ss_max_pulimit insideextra_functionality, rationed by a rolling energy budget
Related user guide pages: Dispatch Limits (section "Voltage Angle Limits"), Linearised Power Flow and Contingencies.
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pypsa
A Small Meshed Network¶
Three buses form a triangle. A cheap generator sits at bus0, an expensive one at bus1, and the load at bus2. All lines have a rating of 100 MVA, but line2 (the direct path from bus0 to bus2) has the lowest reactance. In the linearised power flow, the power from gen0 splits in inverse proportion to the reactance of the two paths: 35/45 of it takes line2, the remaining 10/45 the detour via bus1. Hence line2 is the line that binds first whenever gen0 runs at full output.
We wrap the network in a function so that we can start from a clean copy later on, and define two small helpers that we reuse throughout the example: one plots the line loading on the network map, the other returns the voltage angle difference across each line in degrees. The angles are available in n.buses_t.v_ang right after n.optimize(), no separate power flow is required.
def build_network(load: float = 120) -> pypsa.Network:
n = pypsa.Network()
n.set_snapshots(range(4))
n.add("Carrier", ["AC", "coal", "gas"])
n.add(
"Bus",
["bus0", "bus1", "bus2"],
v_nom=380,
carrier="AC",
x=[0, 2, 1],
y=[0, 0, 1.7],
)
n.add(
"Line", "line0", bus0="bus0", bus1="bus1", x=15, r=1.5, s_nom=100, carrier="AC"
)
n.add(
"Line", "line1", bus0="bus1", bus1="bus2", x=20, r=2.0, s_nom=100, carrier="AC"
)
n.add(
"Line", "line2", bus0="bus0", bus1="bus2", x=10, r=1.0, s_nom=100, carrier="AC"
)
n.add("Generator", "gen0", bus="bus0", p_nom=300, marginal_cost=10, carrier="coal")
n.add("Generator", "gen1", bus="bus1", p_nom=300, marginal_cost=30, carrier="gas")
n.add("Load", "load0", bus="bus2", p_set=load)
n.calculate_dependent_values()
return n
def plot_loading(n: pypsa.Network, ax: plt.Axes, title: str, snapshot: int = 0) -> None:
flow = n.lines_t.p0.loc[snapshot]
loading = flow.abs() / n.lines.s_nom_opt
collections = n.plot(
ax=ax,
geomap=False,
margin=0.3,
bus_size=0.012,
line_width=4,
line_color=loading,
line_cmap="viridis",
line_cmap_norm=plt.Normalize(0, 1),
line_flow=flow / 100,
)
ax.figure.colorbar(
collections["branches"]["Line"], ax=ax, label="Line loading [p.u.]"
)
centroid = n.buses[["x", "y"]].mean()
for bus, row in n.buses.iterrows():
ax.annotate(
bus, (row.x, row.y), textcoords="offset points", xytext=(0, 14), ha="center"
)
for line, row in n.lines.iterrows():
mid = (
n.buses.loc[row.bus0, ["x", "y"]] + n.buses.loc[row.bus1, ["x", "y"]]
) / 2
shift = 0.2 * (mid - centroid) / np.hypot(*(mid - centroid))
ax.annotate(
f"{line}\n{flow[line]:.1f} MW",
(mid.x + shift.x, mid.y + shift.y),
ha="center",
va="center",
fontsize=9,
)
ax.set_title(title)
def angle_differences(n: pypsa.Network) -> pd.DataFrame:
v_ang = np.rad2deg(n.buses_t.v_ang)
return pd.DataFrame(
{line: v_ang[row.bus0] - v_ang[row.bus1] for line, row in n.lines.iterrows()}
)
n = build_network()
n.optimize(solver_name="highs", output_flag=False)
fig, ax = plt.subplots(figsize=(7, 5))
plot_loading(n, ax, "Unconstrained dispatch")
fig.tight_layout()
/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) /tmp/ipykernel_3321/848909555.py:74: 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. n.optimize(solver_name="highs", output_flag=False)
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options: - output_flag: False - log_to_console: False
INFO:linopy.io: Writing time: 0.03s
INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 20 primals, 56 duals Objective: 4.80e+03 Solver: highs Runtime: 0.00s MIP gap: inf Dual bound: 0.00e+00 Solver model: available Solver message: Optimal
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Kirchhoff-Voltage-Law were not assigned to the network.
Without any additional limit the cheap gen0 covers the whole load and line2 runs at 93% of its rating. This is the reference dispatch against which the following limits are measured.
1. Static N-1 Approximation via s_max_pu¶
The simplest security margin is to derate every line below its thermal limit, e.g. to 70% of s_nom, by setting a static Line.s_max_pu. The headroom makes it less likely that the remaining lines overload when a parallel circuit trips.
This is only a cheap approximation of N-1 security. It does not check any specific outage, and a well-loaded network can still violate limits after a real contingency. The exact method, which enforces the line ratings for every single branch outage, is shown in the Security-Constrained LOPF example and described in the Contingencies user guide.
n.lines["s_max_pu"] = 0.7
n.optimize(solver_name="highs", output_flag=False)
fig, ax = plt.subplots(figsize=(7, 5))
plot_loading(n, ax, "Static derating to 70% of s_nom")
fig.tight_layout()
loading = n.lines_t.p0.loc[0].abs() / n.lines.s_nom
loading.rename("loading [p.u.]").to_frame().T
/tmp/ipykernel_3321/2052764756.py:2: 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. n.optimize(solver_name="highs", output_flag=False)
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options: - output_flag: False - log_to_console: False
INFO:linopy.io: Writing time: 0.02s
INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 20 primals, 56 duals Objective: 1.04e+04 Solver: highs Runtime: 0.00s MIP gap: inf Dual bound: 0.00e+00 Solver model: available Solver message: Optimal
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Kirchhoff-Voltage-Law were not assigned to the network.
| name | line0 | line1 | line2 |
|---|---|---|---|
| loading [p.u.] | 0.2 | 0.5 | 0.7 |
line2 is now loaded at exactly 70%. To respect the cap the optimiser shifts 70 MW of generation to the expensive gen1, 20 MW of which flow backwards over line0 to bus0. The static derating scales with s_nom: a line with twice the rating gets twice the headroom.
2. Dynamic Line Rating via Time-Varying s_max_pu¶
Overhead line ratings depend on ambient conditions. A well-cooled line can carry more current than its static rating suggests, while a hot, still day requires a lower limit. Such a dynamic line rating (DLR) is modelled by supplying a per-snapshot series in n.lines_t.s_max_pu, which takes precedence over the static Line.s_max_pu for the lines it covers.
We give line2 a rating that varies over the four snapshots, mimicking changing weather.
n.lines["s_max_pu"] = 1.0
n.lines_t.s_max_pu["line2"] = pd.Series([0.9, 0.8, 0.7, 0.85], index=n.snapshots)
n.optimize(solver_name="highs", output_flag=False)
flow = n.lines_t.p0["line2"]
rating = n.lines_t.s_max_pu["line2"] * n.lines.at["line2", "s_nom"]
fig, ax = plt.subplots(figsize=(7, 4))
ax.bar(n.snapshots, flow, color="tab:blue", alpha=0.7, label="Flow on line2")
ax.step(
n.snapshots,
rating,
where="mid",
color="tab:red",
linewidth=2,
label="Dynamic rating",
)
ax.axhline(n.lines.at["line2", "s_nom"], color="grey", linestyle="--", label="s_nom")
ax.set_xlabel("Snapshot")
ax.set_ylabel("Power [MW]")
ax.set_xticks(n.snapshots)
ax.set_title("Dynamic line rating on line2")
ax.legend()
ax.grid(alpha=0.3)
fig.tight_layout()
pd.DataFrame({"flow [MW]": flow, "rating [MW]": rating})
/tmp/ipykernel_3321/1373415494.py:3: 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. n.optimize(solver_name="highs", output_flag=False)
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options: - output_flag: False - log_to_console: False
INFO:linopy.io: Writing time: 0.02s
INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 20 primals, 56 duals Objective: 7.70e+03 Solver: highs Runtime: 0.00s MIP gap: inf Dual bound: 0.00e+00 Solver model: available Solver message: Optimal
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Kirchhoff-Voltage-Law were not assigned to the network.
| flow [MW] | rating [MW] | |
|---|---|---|
| snapshot | ||
| 0 | 90.0 | 90.0 |
| 1 | 80.0 | 80.0 |
| 2 | 70.0 | 70.0 |
| 3 | 85.0 | 85.0 |
In every snapshot the dispatch is rerouted just enough that the flow on line2 tracks its time-varying rating. The line's nominal s_nom is untouched, the rating only tells the optimiser how much of it is usable at each point in time.
3. Voltage Angle Difference Limits¶
In the linearised power flow, the voltage angle difference across a line is proportional to its flow:
$$\theta_{\text{bus0}} - \theta_{\text{bus1}} = x_{\text{pu,eff}} \cdot s$$
Here $x_{\text{pu,eff}}$ is the effective per-unit reactance (n.lines.x_pu_eff, computed by n.calculate_dependent_values()) and $s$ is the flow in MW, since PyPSA uses a base power of 1 MVA. A maximum angle difference $\bar\theta$ in degrees therefore translates into a flow cap, which we can express as an equivalent s_max_pu:
$$s_{\text{max,pu}} = \frac{\mathrm{deg2rad}(\bar\theta)}{x_{\text{pu,eff}} \cdot s_{\text{nom}}}$$
To see the mechanism, we first apply this conversion by hand and cap line2 at an angle difference of 0.3 degrees. The time-varying rating from the previous section is removed again.
n.lines_t.s_max_pu = n.lines_t.s_max_pu.drop(columns="line2")
v_ang_max = 0.3
x_pu_eff = n.lines.at["line2", "x_pu_eff"]
angle_cap = np.deg2rad(v_ang_max) / x_pu_eff
n.lines.loc["line2", "s_max_pu"] = angle_cap / n.lines.at["line2", "s_nom"]
print(f"x_pu_eff = {x_pu_eff:.3e} p.u.")
print(
f"flow cap = {angle_cap:.1f} MW, equivalent s_max_pu = {n.lines.at['line2', 's_max_pu']:.3f}"
)
n.optimize(solver_name="highs", output_flag=False)
manual = angle_differences(n).loc[0]
manual.rename("angle difference [deg]").to_frame().T
/tmp/ipykernel_3321/340905905.py:13: 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. n.optimize(solver_name="highs", output_flag=False)
x_pu_eff = 6.925e-05 p.u. flow cap = 75.6 MW, equivalent s_max_pu = 0.756
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options: - output_flag: False - log_to_console: False
INFO:linopy.io: Writing time: 0.02s
INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 20 primals, 56 duals Objective: 9.05e+03 Solver: highs Runtime: 0.00s MIP gap: inf Dual bound: 0.00e+00 Solver model: available Solver message: Optimal
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Kirchhoff-Voltage-Law were not assigned to the network.
| line0 | line1 | line2 | |
|---|---|---|---|
| angle difference [deg] | -0.052284 | 0.352284 | 0.3 |
The angle difference across line2 is exactly 0.3 degrees, so the manual conversion enforces the desired limit. The flow cap is only 75.6 MW here, i.e. the angle limit is stricter than the thermal rating of 100 MVA.
Native v_ang_max¶
The same limit can be set directly via Line.v_ang_max (in degrees). The optimisation then adds the constraints Line-v_ang-lower and Line-v_ang-upper, which bound the flow $s$ by $\pm\,\mathrm{deg2rad}(\bar\theta) / x_{\text{pu,eff}}$, exactly the cap derived above.
Like s_max_pu, the limit is symmetric: v_ang_max bounds the magnitude of $\theta_{\text{bus0}} - \theta_{\text{bus1}}$, so both flow directions are capped.
n.lines.loc["line2", "s_max_pu"] = 1.0
n.lines.loc["line2", "v_ang_max"] = v_ang_max
n.optimize(solver_name="highs", output_flag=False)
n.model.constraints["Line-v_ang-upper"]
/tmp/ipykernel_3321/1490400602.py:3: 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. n.optimize(solver_name="highs", output_flag=False)
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options: - output_flag: False - log_to_console: False
INFO:linopy.io: Writing time: 0.03s
INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 20 primals, 64 duals Objective: 9.05e+03 Solver: highs Runtime: 0.00s MIP gap: inf Dual bound: 0.00e+00 Solver model: available Solver message: Optimal
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Kirchhoff-Voltage-Law, Line-v_ang-lower, Line-v_ang-upper were not assigned to the network.
Constraint `Line-v_ang-upper` [snapshot: 4, name: 1]: ----------------------------------------------------- [0, line2]: +1 Line-s[0, line2] ≤ 75.60766319639434 [1, line2]: +1 Line-s[1, line2] ≤ 75.60766319639434 [2, line2]: +1 Line-s[2, line2] ≤ 75.60766319639434 [3, line2]: +1 Line-s[3, line2] ≤ 75.60766319639434
native = angle_differences(n).loc[0]
fig, ax = plt.subplots(figsize=(7, 4))
pd.DataFrame({"manual s_max_pu": manual, "native v_ang_max": native}).plot.bar(
ax=ax, color=["tab:blue", "tab:orange"], rot=0
)
ax.axhline(v_ang_max, color="tab:red", linestyle="--", label="v_ang_max")
ax.axhline(-v_ang_max, color="tab:red", linestyle=":", label="-v_ang_max")
ax.axhline(0, color="black", linewidth=0.8)
ax.set_ylim(-0.45, 0.45)
ax.set_ylabel("Voltage angle difference [deg]")
ax.set_xlabel("")
ax.set_title("Angle difference across lines (snapshot 0)")
ax.legend(loc="upper left", ncol=1, bbox_to_anchor=(1, 1), frameon=False)
ax.grid(alpha=0.3, axis="y")
fig.tight_layout()
The native constraint reproduces the manual result: line2 sits exactly at its 0.3 degree limit, while the other two lines are unconstrained and only follow from the network's Kirchhoff laws.
!!! info "Native voltage angle difference limits"
Enforcement of `Line.v_ang_max` and `Transformer.v_ang_max` in the optimisation is implemented natively (see [GitHub issue #1481](https://github.com/PyPSA/PyPSA/issues/1481)). The manual `s_max_pu` conversion above remains a useful way to understand the mechanism. For a transformer the angle difference additionally couples with its optimisable `phase_shift`.4. Angle Limits with Extendable Capacities¶
The two kinds of limit behave very differently once s_nom becomes a decision variable. A derating s_max_pu is a fraction of the optimised capacity s_nom_opt, so it grows with expansion. The angle limit does not: it is an absolute cap $\mathrm{deg2rad}(\bar\theta) / x_{\text{pu,eff}}$ in MW, and the reactance is a fixed parameter in the LP. The linear optimal power flow assumes that the impedance does not change with s_nom. Note that the manual conversion into s_max_pu from the previous section is no longer valid here: the product s_max_pu * s_nom would have to shrink as s_nom grows, which is not linear. The native constraint bounds the flow directly and contains no s_nom, so it stays linear.
Physically, adding parallel circuits does lower the reactance and relaxes the angle cap. Within a single LP the reactance stays fixed, so the cap cannot be relaxed by investment. To let the reactance, and with it the angle cap, follow the expanded capacity, use n.optimize.optimize_transmission_expansion_iteratively(), which rescales x by s_nom_prev / s_nom_opt after each iteration (for extendable lines without a type).
For the experiment we raise the load to 200 MW so that expansion pays off, make line2 extendable at a small capital cost and give it a symmetric angle limit of 0.5 degrees. Without any angle limit the optimiser would expand line2 to 155.6 MW and serve the whole load from gen0.
v_ang_max = 0.5
def expansion_network() -> pypsa.Network:
n = build_network(load=200)
n.lines.loc["line2", ["s_nom_extendable", "s_nom_min", "capital_cost"]] = [
True,
100,
5,
]
n.lines.loc["line2", "v_ang_max"] = v_ang_max
return n
n = expansion_network()
n.optimize(solver_name="highs", output_flag=False)
angle_cap = np.deg2rad(v_ang_max) / n.lines.at["line2", "x_pu_eff"]
print(f"angle cap on line2: {angle_cap:.1f} MW")
pd.DataFrame(
{
"s_nom [MW]": n.lines.s_nom,
"s_nom_opt [MW]": n.lines.s_nom_opt,
"flow [MW]": n.lines_t.p0.loc[0],
"angle difference [deg]": angle_differences(n).loc[0],
}
)
/tmp/ipykernel_3321/3741719373.py:16: 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. n.optimize(solver_name="highs", output_flag=False)
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options: - output_flag: False - log_to_console: False
INFO:linopy.io: Writing time: 0.04s
INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 22 primals, 65 duals Objective: 1.52e+04 Solver: highs Runtime: 0.00s MIP gap: inf Dual bound: 0.00e+00 Solver model: available Solver message: Optimal
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Line-ext-s-lower, Line-ext-s-upper, Kirchhoff-Voltage-Law, Line-v_ang-lower, Line-v_ang-upper were not assigned to the network.
angle cap on line2: 126.0 MW
| s_nom [MW] | s_nom_opt [MW] | flow [MW] | angle difference [deg] | |
|---|---|---|---|---|
| line0 | 100.0 | 100.000000 | -14.641123 | -0.087141 |
| line1 | 100.0 | 100.000000 | 73.987228 | 0.587141 |
| line2 | 100.0 | 126.012772 | 126.012772 | 0.500000 |
n.generators_t.p.loc[0].rename("dispatch [MW]").to_frame().T
| name | gen0 | gen1 |
|---|---|---|
| dispatch [MW] | 111.371649 | 88.628351 |
The optimiser expands line2 to 126.0 MW, exactly the angle cap, and not a single MW further. Extra capacity would cost money without allowing extra flow, because the cap is fixed by the reactance and not by s_nom_opt. As a consequence 88.6 MW of the load still have to come from the expensive gen1, whose injection at bus1 reaches bus2 over line1.
The shadow price of the binding constraint shows how much the angle limit costs per MW and snapshot:
n.model.constraints["Line-v_ang-upper"].dual.to_series().unstack()
| name | line2 |
|---|---|
| snapshot | |
| 0 | -55.0 |
| 1 | -60.0 |
| 2 | -60.0 |
| 3 | -60.0 |
5. Temporary Admissible Loading (TATL)¶
The limits so far are hard caps that hold in every snapshot. Operators also work with two ratings at once: a PATL (permanently admissible transmission loading) that the line may carry indefinitely, and a higher TATL (temporarily admissible transmission loading) that it may carry only for a short time. Loading in the band between the two is still admissible, not an overload, but it heats the conductor, so the loading above the PATL is rationed by an energy budget; sustained loading in the band would drive the conductor temperature past its design limit (see arXiv:2112.06667).
PyPSA has no native attribute for this. We add it through the extra_functionality hook, the general way to attach custom limits to the optimisation. For line2 we introduce a non-negative variable $e_{l,t}$ for the loading above the PATL in each snapshot and
- relax the permanent cap by this excess, $|s_{l,t}| \leq \text{PATL}_l + e_{l,t}$,
- cap the excess at $\text{TATL}_l - \text{PATL}_l = 30$ MW so the flow never exceeds the 100 MW TATL,
- bound the excess energy by a budget $E_l$ over a rolling window of $W$ snapshots, $\sum_{\tau=t-W+1}^{t} w_\tau \, e_{l,\tau} \leq E_l$, with $w_t$ the snapshot weightings in hours.
$$e_{l,t} \geq 0, \qquad -\text{PATL}_l - e_{l,t} \leq s_{l,t} \leq \text{PATL}_l + e_{l,t}, \qquad \sum_{\tau=t-W+1}^{t} w_\tau \, e_{l,\tau} \leq E_l$$
How s_max_pu enters. Here s_max_pu is the PATL. line2 keeps the s_max_pu = 0.7 of the reference solve, so PyPSA already builds its flow limits $-0.7\,s_\text{nom} \leq s_{l,t} \leq 0.7\,s_\text{nom}$ as the constraints Line-fix-s-upper and Line-fix-s-lower. Instead of adding new constraints, we relax these two in place: subtracting $e_{l,t}$ from the upper limit and adding it to the lower one lets the flow surpass the s_max_pu rating by up to $e_{l,t}$. Capping $e_{l,t}$ at $\text{TATL} - \text{PATL} = 30$ MW reinstates the 100 MW TATL as the hard ceiling, and the rolling budget rations how long the band is used. The permanent rating thus stays a genuine s_max_pu that is briefly surpassed, not a separate soft limit. The constraint is named Line-fix-s-... for non-extendable lines; extendable lines use Line-ext-s-..., and lines with transmission losses carry extra loss terms, so adapt the name to your case.
How the budget E is chosen. $E_l$ is the excess loading above the PATL, integrated over time, that the conductor can absorb before its temperature reaches the design limit. It follows from a thermal model of the conductor (its heat capacity and the temperature margin between the PATL and TATL steady states) or, in practice, from an operational rule such as "the TATL may be used for 20 minutes". We use a small illustrative value below.
We give line2 a permanent rating of 70 MW, a temporary ceiling of 100 MW, and a peaky load so that using the band pays off. As a reference we first solve with the PATL alone, set as a static s_max_pu of 0.7.
n = build_network()
n.loads_t.p_set["load0"] = pd.Series([90.0, 125.0, 95.0, 90.0], index=n.snapshots)
n.lines["s_max_pu"] = 1.0
n.lines.loc["line2", "s_max_pu"] = 0.7
n.optimize(solver_name="highs", output_flag=False)
patl_cost = n.objective
pd.DataFrame(
{
"load [MW]": n.loads_t.p_set["load0"],
"line2 flow [MW]": n.lines_t.p0["line2"],
"gen1 [MW]": n.generators_t.p["gen1"],
}
).round(1)
/tmp/ipykernel_3321/1654114155.py:5: 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. n.optimize(solver_name="highs", output_flag=False)
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options: - output_flag: False - log_to_console: False
INFO:linopy.io: Writing time: 0.02s
INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 20 primals, 56 duals Objective: 5.87e+03 Solver: highs Runtime: 0.00s MIP gap: inf Dual bound: 0.00e+00 Solver model: available Solver message: Optimal
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Kirchhoff-Voltage-Law were not assigned to the network.
| load [MW] | line2 flow [MW] | gen1 [MW] | |
|---|---|---|---|
| snapshot | |||
| 0 | 90.0 | 70.0 | -0.0 |
| 1 | 125.0 | 70.0 | 81.7 |
| 2 | 95.0 | 70.0 | 11.7 |
| 3 | 90.0 | 70.0 | -0.0 |
line2 sits at its 70 MW PATL in every snapshot. During the load peak in snapshot 1 the cheap gen0 cannot push more power through the line, so the expensive gen1 has to cover the difference. This is the reference cost against which the temporary loading into the PATL–TATL band is measured.
PATL = 70.0
budget = 20.0
window = 2
def add_tatl(n: pypsa.Network, snapshots: pd.Index) -> None:
m = n.model
TATL = n.lines.at["line2", "s_nom"]
excess = m.add_variables(
lower=0,
upper=TATL - PATL,
coords=m["Line-s"].sel(name=["line2"]).coords,
name="Line-excess",
)
for name, sign in [("Line-fix-s-upper", -1), ("Line-fix-s-lower", 1)]:
con = m.constraints[name]
term = (sign * excess).reindex(name=con.lhs.indexes["name"], fill_value=0)
con.update(lhs=con.lhs + term)
w = n.snapshot_weightings.objective
windowed = (excess.sel(name="line2", drop=True) * w).rolling(snapshot=window).sum()
m.add_constraints(windowed <= budget, name="tatl-budget")
n.optimize(solver_name="highs", output_flag=False, extra_functionality=add_tatl)
excess = n.model.solution["Line-excess"].sel(name="line2", drop=True).to_series()
pd.DataFrame(
{
"load [MW]": n.loads_t.p_set["load0"],
"line2 flow [MW]": n.lines_t.p0["line2"],
"excess above PATL [MW]": excess,
"gen1 [MW]": n.generators_t.p["gen1"],
}
).round(1)
/tmp/ipykernel_3321/1016986648.py:24: 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. n.optimize(solver_name="highs", output_flag=False, extra_functionality=add_tatl)
INFO:linopy.model: Solve problem using Highs solver
INFO:linopy.model:Solver options: - output_flag: False - log_to_console: False
INFO:linopy.io: Writing time: 0.03s
INFO:linopy.constants: Optimization successful: Status: ok Termination condition: optimal Solution: 24 primals, 60 duals Objective: 4.67e+03 Solver: highs Runtime: 0.00s MIP gap: inf Dual bound: 0.00e+00 Solver model: available Solver message: Optimal
INFO:pypsa.optimization.optimize:The shadow-prices of the constraints Generator-fix-p-lower, Generator-fix-p-upper, Line-fix-s-lower, Line-fix-s-upper, Kirchhoff-Voltage-Law, tatl-budget were not assigned to the network.
/home/runner/work/PyPSA/PyPSA/.venv/lib/python3.13/site-packages/linopy/common.py:306: UserWarning: Coordinates across variables not equal. Perform outer join. warn(
| load [MW] | line2 flow [MW] | excess above PATL [MW] | gen1 [MW] | |
|---|---|---|---|---|
| snapshot | ||||
| 0 | 90.0 | 70.0 | 0.0 | -0.0 |
| 1 | 125.0 | 86.1 | 16.1 | 33.3 |
| 2 | 95.0 | 73.9 | 3.9 | -0.0 |
| 3 | 90.0 | 70.0 | 0.0 | -0.0 |
fig, ax = plt.subplots(figsize=(7, 4))
flow = n.lines_t.p0["line2"]
ax.bar(n.snapshots, flow, color="tab:blue", alpha=0.7, label="Flow on line2")
ax.axhline(PATL, color="tab:green", linestyle="--", label="PATL (permanent)")
ax.axhline(
n.lines.at["line2", "s_nom"],
color="tab:red",
linestyle="--",
label="TATL (ceiling)",
)
ax.set_xlabel("Snapshot")
ax.set_ylabel("Power [MW]")
ax.set_xticks(n.snapshots)
ax.set_title(
f"line2 loading with TATL budget ({budget:.0f} MWh per {window}-snapshot window)"
)
ax.legend()
ax.grid(alpha=0.3)
fig.tight_layout()
The optimiser uses the band where it saves the most: about 16 MW of excess at the load peak in snapshot 1 and 4 MW in snapshot 2, so line2 carries roughly 86 MW at the peak. There the cheap gen0 displaces the expensive gen1, while in the two slack snapshots the flow stays at the 70 MW PATL.
The budget applies over a rolling window of window = 2 snapshots rather than the whole horizon: the weighted excess in any two consecutive snapshots must stay below the 20 MWh budget. This mirrors the physics, the conductor cools between snapshots, so the relevant limit is the recent excess energy, not the total over the year; it is the same rolling-window pattern used for maintenance-scheduling constraints. Here the window covering snapshots 1 and 2 binds.
The dual of the budget constraint prices the excess allowance: relaxing the binding window's budget by one MWh lowers the operating cost by that shadow price.
print(f"cost with PATL only: {patl_cost:,.0f}")
print(f"cost with TATL budget: {n.objective:,.0f}")
n.model.constraints["tatl-budget"].dual
cost with PATL only: 5,867 cost with TATL budget: 4,667
<xarray.DataArray 'dual' (snapshot: 4)> Size: 32B array([ -0., -0., -60., -0.]) Coordinates: * snapshot (snapshot) int64 32B 0 1 2 3