Extra material: Energy dispatch problem#
Preamble: Install Pyomo and a solver#
The following cell sets and verifies a global SOLVER for the notebook. If run on Google Colab, the cell installs Pyomo and the HiGHS solver, while, if run elsewhere, it assumes Pyomo and HiGHS have been previously installed. It then sets to use HiGHS as solver via the appsi module and a test is performed to verify that it is available. The solver interface is stored in a global object SOLVER for later use.
import sys
if 'google.colab' in sys.modules:
%pip install pyomo >/dev/null 2>/dev/null
%pip install highspy >/dev/null 2>/dev/null
solver = 'appsi_highs'
import pyomo.environ as pyo
SOLVER = pyo.SolverFactory(solver)
assert SOLVER.available(), f"Solver {solver} is not available."
Problem description#
Optimal Power Flow (OPF) minimizes electricity dispatch costs while meeting demand and respecting generation and transmission limits. We solve 96 independent 15-minute instances, then compare two policies that restrict fossil-fuel generation.
The final section adds a line-switching variant, allowing a limited number of existing lines to be opened while retaining a connected fixed backbone.
Background: Power networks and power flow physics#
We represent the transmission network by a directed graph \(G=(V, E)\), with nodes \(V\) for generators and loads and edges \(E\) for physical transmission lines.
Each node \(i \in V\) has injection \(p_i\) and demand \(d_i\). Generators \(\mathcal{G} \subseteq V\) are nodes \(i \in V\) with \(p_i \geq 0\) and \(d_i = 0\); each \(i \in \mathcal{G}\) has output bounds \(p_i^{\min}\) and \(p_i^{\max}\). Load nodes \(\mathcal{D} \subseteq V\) have \(p_i = 0\) and \(d_i \geq 0\). Supply and demand match when \(\sum_{i \in V}(p_i-d_i)=0\).
Each line \((i, j) \in E\) has signed flow \(f_{ij} \in \mathbb{R}\) and capacity \(f_{ij}^{\max} \geq 0\). Its orientation \((i,j)\) is a reference: \(f_{ij} < 0\) means flow from \(j\) to \(i\), with \(f_{ji} = -f_{ij}\). Capacity applies in either direction: \(|f_{ij}| = |f_{ji}| \leq f_{ij}^{\max}\).
Unlike generic network flows, power flows must also satisfy physical equations. Ignoring capacity limits for the moment, let \(\theta_{i} \in \mathbb{R}\) be the phase angle at node \(i\) and, for each line \((i,j)\), let \(b_{ij} > 0\) be its susceptance. Given balanced injections, \(\sum_{i \in V}(p_i-d_i)=0\), the flows \(\mathbf{f} \in \mathbb{R}^{m}\) and angles \(\mathbf{\theta} \in \mathbb{R}^{n}\) solve
The first set of constraints ensures flow conservation, and the second relates flows to susceptances and angle differences. For a connected physical network with positive susceptances, matched injections and demand give a unique flow vector. Phase angles are unique up to an additive constant, which can be fixed by choosing one reference angle.
The weighted network Laplacian with edge weights \(b_{ij}\) has rank \(|V|-1\) in a connected network. If the network is disconnected, supply and demand must balance separately in each connected component, and one reference angle is needed per component. Global balance alone is then insufficient.
Optimal power flow problem#
We now choose the injections \(\mathbf{p}\) instead of treating them as given. OPF finds a minimum-cost solution \((\mathbf{p}, \mathbf{f}, \mathbf{\theta})\) satisfying demand, line capacities and generator limits.
Let \(c_i \geq 0\) be the unit production cost of generator \(i\). The model is
Load nodes have \(p_i = p_i^{\min} = p_i^{\max} = 0\) for all \(i \in \mathcal{D}\), so we can define \(p_i\) for every node. Generators have \(d_i = 0\) for all \(i \in \mathcal{G}\).
The decision variables are:
\(p_i\): injections
\(\theta_i\): phase angles
\(f_{ij}\): line flows
All other quantities are instance parameters. Demand and renewable generation vary between time intervals.
Data#
We use an IEEE-118 instance with 118 nodes and 179 lines. The data are loaded as pandas.DataFrames and then converted to a nested dictionary for the Pyomo models.
Setup code and data import#
The repository includes data/nodes.csv and data/edges.csv alongside this notebook. The cell below reads these local files, whether you run from the notebook’s folder or the repository root. If neither local file is present, as when opening the notebook on Colab without the data folder, it downloads both files from a fixed repository snapshot. Keep the two local CSV files together; an incomplete local pair raises an error rather than mixing local and downloaded data.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import networkx as nx
from ast import literal_eval as make_tuple
from pathlib import Path
# Use the bundled data when available; otherwise download the same snapshot.
data_files = ("nodes.csv", "edges.csv")
data_dir = next(
(folder for folder in (Path("data"), Path("notebooks/04/data"))
if any((folder / name).is_file() for name in data_files)),
None,
)
if data_dir is not None:
if not all((data_dir / name).is_file() for name in data_files):
raise FileNotFoundError(f"Keep nodes.csv and edges.csv together in {data_dir}.")
data_source = str(data_dir)
else:
data_source = (
"https://raw.githubusercontent.com/mobook/MO-book/"
"58832cc0eb3c59c1c013451eaf4ccc0e5e12aae9/notebooks/04/data"
)
nodes_df = pd.read_csv(f"{data_source}/nodes.csv", index_col=0)
edges_df = pd.read_csv(f"{data_source}/edges.csv", index_col=0)
# Read all instances
nodes = nodes_df.groupby("instance").apply(
lambda data: data.set_index("node_id").T.to_dict()
)
edges = edges_df.set_index(edges_df["edge_id"].apply(make_tuple)).T.to_dict()
I = [{"nodes": n, "edges": edges} for n in nodes]
# Initialize a network for demonstration purposes
network = I[0]
example_nodes = nodes_df[nodes_df.instance == 0]
Network data#
def visualize_network(network, edge_flows=None, ax=None):
"""Visualize a network with technology-dependent generator colors and gray loads."""
plt.figure(figsize=(12, 10))
g = nx.DiGraph(network["edges"].keys())
pos = nx.layout.kamada_kawai_layout(g, weight=None)
color_mapping = {
"solar": "#ffcb36",
"wind": "white",
"hydro": "#a5efff",
"coal": "#686868",
"gas": "#00ab4e",
np.nan: "#b6b6b6",
}
vertex2color = {
i: color_mapping[data["energy_type"]] for i, data in network["nodes"].items()
}
v2c_list = [vertex2color[i] for i in g.nodes] # Order based on networkx
nodes = nx.draw_networkx_nodes(
g,
pos,
node_size=250,
node_color=v2c_list,
linewidths=2,
)
edges = nx.draw_networkx_edges(
g,
pos,
width=2,
edge_color="#595959",
)
# Gives node colors
ax = plt.gca()
ax.collections[0].set_edgecolor("#595959")
ax.set_axis_off()
visualize_network(network)
The network has 18 generators: 2 hydro (cyan), 4 coal (dark gray), 4 solar (yellow), 2 wind (white) and 6 gas (green), plus 100 load nodes (gray).
Each node has these attributes:
node_id: identifierd: demandp_min,p_max: production boundsc_var: variable production costis_generator: whether the node is a generatorenergy_type: generation technology
Filter generators using is_generator. They have \(d_i=0\). Hydro, solar and wind have zero c_var; solar and wind have p_min = p_max, so their available production is fixed.
example_nodes[example_nodes.is_generator]
| node_id | d | p_min | p_max | c_var | is_generator | energy_type | instance | |
|---|---|---|---|---|---|---|---|---|
| 9 | 9 | 0.0 | 0.000000 | 400.000000 | 0.000000 | True | hydro | 0 |
| 11 | 11 | 0.0 | 0.000000 | 200.000000 | 0.000000 | True | hydro | 0 |
| 24 | 24 | 0.0 | 0.000000 | 397.800000 | 28.948321 | True | coal | 0 |
| 25 | 25 | 0.0 | 0.000000 | 873.000000 | 22.220980 | True | coal | 0 |
| 30 | 30 | 0.0 | 0.000000 | 612.000000 | 25.993982 | True | coal | 0 |
| 45 | 45 | 0.0 | 0.000000 | 720.000000 | 24.202306 | True | coal | 0 |
| 48 | 48 | 0.0 | 0.000000 | 0.000000 | 0.000000 | True | solar | 0 |
| 53 | 53 | 0.0 | 0.000000 | 0.000000 | 0.000000 | True | solar | 0 |
| 58 | 58 | 0.0 | 0.000000 | 0.000000 | 0.000000 | True | solar | 0 |
| 60 | 60 | 0.0 | 0.000000 | 0.000000 | 0.000000 | True | solar | 0 |
| 64 | 64 | 0.0 | 38.484861 | 38.484861 | 0.000000 | True | wind | 0 |
| 65 | 65 | 0.0 | 180.602933 | 180.602933 | 0.000000 | True | wind | 0 |
| 79 | 79 | 0.0 | 0.000000 | 916.200000 | 50.330000 | True | gas | 0 |
| 86 | 86 | 0.0 | 0.000000 | 360.000000 | 50.330000 | True | gas | 0 |
| 88 | 88 | 0.0 | 0.000000 | 1146.600000 | 50.330000 | True | gas | 0 |
| 99 | 99 | 0.0 | 0.000000 | 1175.400000 | 50.330000 | True | gas | 0 |
| 102 | 102 | 0.0 | 0.000000 | 194.400000 | 50.330000 | True | gas | 0 |
| 110 | 110 | 0.0 | 0.000000 | 142.200000 | 50.330000 | True | gas | 0 |
For the load nodes, the only important parameter is the demand \(d_i \geq 0\). All other parameters are either zero, False, or NaN.
example_nodes[~example_nodes.is_generator]
| node_id | d | p_min | p_max | c_var | is_generator | energy_type | instance | |
|---|---|---|---|---|---|---|---|---|
| 0 | 0 | 28.186145 | 0.0 | 0.0 | 0.0 | False | NaN | 0 |
| 1 | 1 | 10.921628 | 0.0 | 0.0 | 0.0 | False | NaN | 0 |
| 2 | 2 | 22.884795 | 0.0 | 0.0 | 0.0 | False | NaN | 0 |
| 3 | 3 | 21.961682 | 0.0 | 0.0 | 0.0 | False | NaN | 0 |
| 4 | 4 | 0.000000 | 0.0 | 0.0 | 0.0 | False | NaN | 0 |
| ... | ... | ... | ... | ... | ... | ... | ... | ... |
| 113 | 113 | 4.396000 | 0.0 | 0.0 | 0.0 | False | NaN | 0 |
| 114 | 114 | 11.864042 | 0.0 | 0.0 | 0.0 | False | NaN | 0 |
| 115 | 115 | 108.311328 | 0.0 | 0.0 | 0.0 | False | NaN | 0 |
| 116 | 116 | 10.703998 | 0.0 | 0.0 | 0.0 | False | NaN | 0 |
| 117 | 117 | 18.242759 | 0.0 | 0.0 | 0.0 | False | NaN | 0 |
100 rows × 8 columns
Edges#
Each of the 179 lines has edge_id (its node pair), b (susceptance) and f_max (capacity).
edges_df
| edge_id | b | f_max | |
|---|---|---|---|
| 0 | (0, 1) | 10.0100 | 271 |
| 1 | (0, 2) | 23.5849 | 271 |
| 2 | (3, 4) | 125.3133 | 316 |
| 3 | (2, 4) | 9.2593 | 315 |
| 4 | (4, 5) | 18.5185 | 316 |
| ... | ... | ... | ... |
| 174 | (64, 65) | 28.9059 | 1427 |
| 175 | (67, 68) | 28.9059 | 1427 |
| 176 | (80, 79) | 28.9059 | 1427 |
| 177 | (86, 85) | 4.8216 | 253 |
| 178 | (67, 115) | 246.9136 | 12992 |
179 rows × 3 columns
Network data structure#
The models access attributes through nested dictionaries:
network = {"nodes": {i: {"attr1": 1, ...},
...
},
"edges": {(i, j): {"attr1": 4, "attr2": 1, ...},
...
}
}
network["nodes"][i] and network["edges"][i, j] contain the attributes of a node and line. For example,
network["nodes"][1]["d"]
10.921627803299574
returns demand at node 1. To list generators:
[i for i, data in network["nodes"].items() if data["is_generator"]]
[9, 11, 24, 25, 30, 45, 48, 53, 58, 60, 64, 65, 79, 86, 88, 99, 102, 110]
Instances#
The list I contains 96 independent 15-minute instances spanning one day, with varying demand and renewable generation. I[0] covers [00:00, 00:15], and I[95] is the last interval. Decisions are not coupled across intervals.
Solving the OPF problem#
We split each flow into nonnegative parts \(f^+_{ij}\) and \(f^-_{ij}\) to express the absolute-flow limits linearly.
We implement this model in Pyomo, solve all 96 instances and report their average objective value.
def OPF1(network):
"""
Input:
- network: a dictionary containing:
- a dictionary of nodes with a dictionary of attributes
- a dictionary of edges with a dictionary of attributes
Output:
- the optimal value of the objective function
"""
# Define a model
model = pyo.ConcreteModel("Q1")
# Declare decision variables
model.p = pyo.Var(network["nodes"], domain=pyo.NonNegativeReals)
model.theta = pyo.Var(network["nodes"], domain=pyo.Reals)
model.fabs = pyo.Var(network["edges"], domain=pyo.NonNegativeReals)
model.fp = pyo.Var(network["edges"], domain=pyo.NonNegativeReals)
model.fm = pyo.Var(network["edges"], domain=pyo.NonNegativeReals)
# Declare objective value
@model.Objective(sense=pyo.minimize)
def objective(m):
return sum(
data["c_var"] * m.p[i]
for i, data in network["nodes"].items()
if data["is_generator"]
)
@model.Expression(network["nodes"])
def outgoing_flow(m, i):
return sum(
m.fp[i, j] - m.fm[i, j]
for j in network["nodes"]
if (i, j) in network["edges"]
)
@model.Expression(network["nodes"])
def incoming_flow(m, i):
return sum(
m.fp[j, i] - m.fm[j, i]
for j in network["nodes"]
if (j, i) in network["edges"]
)
@model.Constraint(network["nodes"])
def flow_conservation(m, i):
return (
m.outgoing_flow[i] - m.incoming_flow[i] == m.p[i] - network["nodes"][i]["d"]
)
@model.Constraint(network["edges"])
def susceptance(m, i, j):
return m.fp[i, j] - m.fm[i, j] == network["edges"][(i, j)]["b"] * (
m.theta[i] - m.theta[j]
)
@model.Constraint(network["edges"])
def abs_flow(m, *e):
return m.fabs[e] == m.fp[e] + m.fm[e]
@model.Constraint(network["edges"])
def flows_upper_bound(m, *e):
return m.fabs[e] <= network["edges"][e]["f_max"]
@model.Constraint(network["nodes"])
def generation_upper_bound(m, i):
return m.p[i] <= network["nodes"][i]["p_max"]
@model.Constraint(network["nodes"])
def generation_lower_bound(m, i):
return network["nodes"][i]["p_min"] <= m.p[i]
# Solve
result = SOLVER.solve(model)
pyo.assert_optimal_termination(result)
return model.objective()
sol1 = [OPF1(instance) for instance in I]
print(f"The average objective value over all instances is: {np.mean(sol1):.2f}")
The average objective value over all instances is: 38507.23
Strict fossil fuel policy, part 1#
In this policy scenario, at most 2 gas plants and 1 coal plant may be active in each interval; inactive plants cannot produce electricity. The variables \(x_i\), \(i\in V\), indicate whether generator \(i\) is active:
As before, we solve all 96 instances and report the average objective value.
def OPF2(network):
"""
Input:
- network: a dictionary containing:
- a dictionary of nodes with a dictionary of attributes
- a dictionary of edges with a dictionary of attributes
Output:
- the optimal value of the objective function
"""
# Define a model
model = pyo.ConcreteModel("Q2")
# Declare decision variables
model.p = pyo.Var(network["nodes"], domain=pyo.NonNegativeReals)
model.theta = pyo.Var(network["nodes"], domain=pyo.Reals)
model.fabs = pyo.Var(network["edges"], domain=pyo.NonNegativeReals)
model.fp = pyo.Var(network["edges"], domain=pyo.NonNegativeReals)
model.fm = pyo.Var(network["edges"], domain=pyo.NonNegativeReals)
model.x = pyo.Var(network["nodes"], domain=pyo.Binary)
# Declare objective value
@model.Objective(sense=pyo.minimize)
def objective(m):
return sum(
data["c_var"] * m.p[i]
for i, data in network["nodes"].items()
if data["is_generator"]
)
# Declare constraints
@model.Expression(network["nodes"])
def outgoing_flow(m, i):
return sum(
m.fp[i, j] - m.fm[i, j]
for j in network["nodes"]
if (i, j) in network["edges"]
)
@model.Expression(network["nodes"])
def incoming_flow(m, i):
return sum(
m.fp[j, i] - m.fm[j, i]
for j in network["nodes"]
if (j, i) in network["edges"]
)
@model.Constraint(network["nodes"])
def flow_conservation(m, i):
return (
m.outgoing_flow[i] - m.incoming_flow[i] == m.p[i] - network["nodes"][i]["d"]
)
@model.Constraint(network["edges"])
def susceptance(m, i, j):
return m.fp[i, j] - m.fm[i, j] == network["edges"][(i, j)]["b"] * (
m.theta[i] - m.theta[j]
)
@model.Constraint(network["edges"])
def abs_flows(m, *e):
return m.fabs[e] == m.fp[e] + m.fm[e]
@model.Constraint(network["edges"])
def flows_upper_bound(m, *e):
return m.fabs[e] <= network["edges"][e]["f_max"]
@model.Constraint(network["nodes"])
def generation_upper_bound(m, i):
return m.p[i] <= network["nodes"][i]["p_max"] * m.x[i]
@model.Constraint(network["nodes"])
def generation_lower_bound(m, i):
return network["nodes"][i]["p_min"] * m.x[i] <= m.p[i]
@model.Constraint()
def max_gas_plants(model):
return (
sum(
model.x[i]
for i, data in network["nodes"].items()
if data["energy_type"] == "gas"
)
<= 2
)
@model.Constraint()
def max_coal_plants(model):
return (
sum(
model.x[i]
for i, data in network["nodes"].items()
if data["energy_type"] == "coal"
)
<= 1
)
# Solve
result = SOLVER.solve(model)
pyo.assert_optimal_termination(result)
return model.objective()
sol2 = [OPF2(instance) for instance in I]
print(f"The average objective value over all instances is: {np.mean(sol2):.2f}")
The average objective value over all instances is: 53120.81
Strict fossil fuel policy, part 2#
We relax the policy when renewable generation is low: the limits of 2 gas plants and 1 coal plant apply only when total solar, wind and hydro production exceeds 1000. Otherwise, the generator-count limits are lifted. We represent these alternatives with big-\(M\) constraints:
We again solve all 96 instances and report the average cost.
def OPF3(network):
"""
Input:
- network: a dictionary containing:
- a dictionary of nodes with a dictionary of attributes
- a dictionary of edges with a dictionary of attributes
Output:
- the optimal value of the objective function
"""
# Define a model
model = pyo.ConcreteModel("Q3")
# Declare decision variables
model.p = pyo.Var(network["nodes"], domain=pyo.NonNegativeReals)
model.theta = pyo.Var(network["nodes"], domain=pyo.Reals)
model.fabs = pyo.Var(network["edges"], domain=pyo.NonNegativeReals)
model.fp = pyo.Var(network["edges"], domain=pyo.NonNegativeReals)
model.fm = pyo.Var(network["edges"], domain=pyo.NonNegativeReals)
model.x = pyo.Var(network["nodes"], domain=pyo.Binary)
model.y = pyo.Var(domain=pyo.Binary)
# Big-Ms
max_total_renewable_production = sum(
data["p_max"]
for i, data in network["nodes"].items()
if data["energy_type"] in ["solar", "wind", "hydro"]
)
M = [4, 3, max_total_renewable_production - 1000]
# Declare objective value
@model.Objective(sense=pyo.minimize)
def objective(m):
return sum(
data["c_var"] * m.p[i]
for i, data in network["nodes"].items()
if data["is_generator"]
)
# Declare constraints
@model.Expression(network["nodes"])
def outgoing_flow(m, i):
return sum(
m.fp[i, j] - m.fm[i, j]
for j in network["nodes"]
if (i, j) in network["edges"]
)
@model.Expression(network["nodes"])
def incoming_flow(m, i):
return sum(
m.fp[j, i] - m.fm[j, i]
for j in network["nodes"]
if (j, i) in network["edges"]
)
@model.Constraint(network["nodes"])
def flow_conservation(m, i):
return (
m.outgoing_flow[i] - m.incoming_flow[i] == m.p[i] - network["nodes"][i]["d"]
)
@model.Constraint(network["edges"])
def susceptance(m, i, j):
return m.fp[i, j] - m.fm[i, j] == network["edges"][(i, j)]["b"] * (
m.theta[i] - m.theta[j]
)
@model.Constraint(network["edges"])
def abs_flows(m, *e):
return m.fabs[e] == m.fp[e] + m.fm[e]
@model.Constraint(network["edges"])
def flows_upper_bound(m, *e):
return m.fabs[e] <= network["edges"][e]["f_max"]
@model.Constraint(network["nodes"])
def generation_upper_bound(m, i):
return m.p[i] <= network["nodes"][i]["p_max"] * m.x[i]
@model.Constraint(network["nodes"])
def generation_lower_bound(m, i):
return network["nodes"][i]["p_min"] * m.x[i] <= m.p[i]
@model.Constraint()
def max_gas_plants(m):
return (
sum(
m.x[i]
for i, data in network["nodes"].items()
if data["energy_type"] == "gas"
)
<= 2 + (1 - m.y) * M[0]
)
@model.Constraint()
def max_coal_plants(m):
return (
sum(
m.x[i]
for i, data in network["nodes"].items()
if data["energy_type"] == "coal"
)
<= 1 + (1 - m.y) * M[1]
)
@model.Constraint()
def renewable_energy_production(m):
return (
sum(
m.p[i]
for i, data in network["nodes"].items()
if data["energy_type"] in ["solar", "wind", "hydro"]
)
<= 1000 + m.y * M[2]
)
# Solve
result = SOLVER.solve(model)
pyo.assert_optimal_termination(result)
return model.objective()
sol3 = [OPF3(instance) for instance in I]
print(f"The average objective value over all instances is: {np.mean(sol3):.2f}")
The average objective value over all instances is: 41608.73
Comparing the three models#
We plot the 96 objective values for each model and interpret their differences below.
objs = [sol1, sol2, sol3]
plt.plot(objs[0], color="blue", label="OPF1")
plt.plot(objs[1], color="green", label="OPF2")
plt.plot(objs[2], color="orange", label="OPF3")
plt.title("Optimal objective values over all instances and models")
plt.xlabel("Instance number")
plt.ylabel("Optimal objective value")
plt.legend()
plt.show()
The goal is to minimize energy-dispatch costs. OPF1 is the base model without plant-activation decisions. OPF2 limits the number of active gas and coal plants, while OPF3 applies those limits only when renewable production exceeds 1000.
Let \(F_2,F_3\) denote the feasible regions projected onto the common production and flow variables. Any OPF2 solution is also feasible for OPF3 by selecting its restricted-generator alternative. Thus \(F_2\subseteq F_3\), and the minimum costs satisfy \(v_3\leq v_2\).
OPF3 can instead keep renewable production at or below 1000 and lift the generator-count limits. Its optimum need not equal either of the other two optima: this renewable-production cap makes its second alternative different from OPF1.
In the 96 instances plotted above, the costs satisfy \(v_1\leq v_3\leq v_2\). The comparison with OPF1 is an observation for these data, rather than a general feasible-region inclusion: OPF2 and OPF3 can deactivate wind or solar plants through \(x_i\), whereas OPF1 enforces their positive minimum-output bounds.
A line-switching variant#
The three models above keep the transmission network fixed. We now allow a limited number of existing lines to be opened. Opening a line sets its flow to zero and removes its phase-angle equation; keeping that equation with zero flow would incorrectly force its endpoint angles to coincide. The signed-flow variables in this section replace the split flows used above.
Let \(S\subseteq E\) contain the switchable lines and let \(z_{ij}=1\) mean that line \((i,j)\in S\) is active. All lines in \(E\setminus S\) remain active. We retain OPF1’s generation bounds, costs and nodal balances, and replace the line constraints by
For fixed lines we retain \(|f_{ij}|\leq f^{\max}_{ij}\) and \(f_{ij}=b_{ij}(\theta_i-\theta_j)\). We fix one reference angle to zero. This is a mixed-integer linear optimization problem, which HiGHS can solve.
For this teaching variant the permanently active lines must form a connected backbone. Along a fixed line, \(|\theta_i-\theta_j|\leq f^{\max}_{ij}/b_{ij}\). Adding these inequalities along any backbone path between the endpoints of a switchable line gives a valid bound on its angle difference. The code uses the shortest such path, weighted by \(f^{\max}/b\), and sets \(M_{ij}\) to its length times \(b_{ij}\). Thus opening a line releases its physical equation without imposing an arbitrary numerical bound. The connected backbone also prevents islanding; choosing all lines freely would require a different formulation.
The function below builds the model. max_open is \(K\), the maximum number of open lines; the initial comparison has all lines active. This represents opening and reclosing existing transmission lines, not building new infrastructure. The plant-activation policies of OPF2 and OPF3 are not included in this extension.
The fixed-backbone path bound is described by Fattahi, Lavaei and Atamtürk in A Bound Strengthening Method for Optimal Transmission Switching in Power Systems, Section IV (2019).
def OPF_switching(network, switchable=(), max_open=1):
"""Build DC optimal transmission switching with a connected fixed backbone."""
nodes, edges = network["nodes"], network["edges"]
switchable = list(switchable)
if len(set(switchable)) != len(switchable) or not set(switchable) <= set(edges):
raise ValueError("Switchable lines must be distinct existing edge keys.")
if not isinstance(max_open, int) or max_open < 0:
raise ValueError("max_open must be a nonnegative integer.")
backbone = nx.Graph()
backbone.add_nodes_from(nodes)
for (i, j), data in edges.items():
if data["b"] <= 0 or data["f_max"] < 0:
raise ValueError("Lines need positive susceptance and nonnegative capacity.")
if (i, j) not in switchable:
backbone.add_edge(i, j, weight=data["f_max"] / data["b"])
if not nx.is_connected(backbone):
raise ValueError("The permanently active lines must connect all nodes.")
big_m = {
(i, j): edges[i, j]["b"] * nx.shortest_path_length(
backbone, i, j, weight="weight"
)
for i, j in switchable
}
model = pyo.ConcreteModel("Optimal transmission switching")
model.V = pyo.Set(initialize=list(nodes))
model.E = pyo.Set(initialize=list(edges), dimen=2)
model.S = pyo.Set(initialize=switchable, dimen=2)
model.p = pyo.Var(model.V, bounds=lambda m, i: (
nodes[i]["p_min"], nodes[i]["p_max"]
))
model.theta = pyo.Var(model.V, domain=pyo.Reals)
model.f = pyo.Var(model.E, domain=pyo.Reals)
model.z = pyo.Var(model.S, domain=pyo.Binary)
model.theta[next(iter(nodes))].fix(0)
model.M = pyo.Param(model.S, initialize=big_m)
@model.Objective(sense=pyo.minimize)
def objective(m):
return sum(nodes[i]["c_var"] * m.p[i] for i in m.V
if nodes[i]["is_generator"])
@model.Constraint(model.V)
def balance(m, i):
return (
sum(m.f[a, b] for a, b in m.E if a == i)
- sum(m.f[a, b] for a, b in m.E if b == i)
== m.p[i] - nodes[i]["d"]
)
@model.Expression(model.E)
def line_status(m, i, j):
return m.z[i, j] if (i, j) in m.S else 1
@model.Constraint(model.E)
def capacity_upper(m, i, j):
return m.f[i, j] <= edges[i, j]["f_max"] * m.line_status[i, j]
@model.Constraint(model.E)
def capacity_lower(m, i, j):
return -edges[i, j]["f_max"] * m.line_status[i, j] <= m.f[i, j]
@model.Constraint(model.E)
def physics_upper(m, i, j):
residual = m.f[i, j] - edges[i, j]["b"] * (m.theta[i] - m.theta[j])
return residual <= (m.M[i, j] * (1 - m.z[i, j]) if (i, j) in m.S else 0)
@model.Constraint(model.E)
def physics_lower(m, i, j):
residual = m.f[i, j] - edges[i, j]["b"] * (m.theta[i] - m.theta[j])
return -residual <= (m.M[i, j] * (1 - m.z[i, j]) if (i, j) in m.S else 0)
@model.Constraint()
def switching_budget(m):
if not switchable:
return pyo.Constraint.Skip
return sum(1 - m.z[e] for e in m.S) <= max_open
return model
# Keep the tighter switching settings separate from the earlier policy solves.
SWITCHING_SOLVER = pyo.SolverFactory("appsi_highs")
SWITCHING_SOLVER.options["mip_rel_gap"] = 0.0
SWITCHING_SOLVER.options["mip_abs_gap"] = 0.0
def solve_switching(model):
"""Request zero MIP gaps before comparing the switching costs."""
results = SWITCHING_SOLVER.solve(model)
pyo.assert_optimal_termination(results)
return pyo.value(model.objective)
A small network#
# A three-bus example with cheap generation at bus 0 and costly generation at bus 2.
triangle = {
"nodes": {
0: {"d": 0, "p_min": 0, "p_max": 100, "c_var": 10, "is_generator": True},
1: {"d": 0, "p_min": 0, "p_max": 0, "c_var": 0, "is_generator": False},
2: {"d": 100, "p_min": 0, "p_max": 100, "c_var": 50, "is_generator": True},
},
"edges": {
(0, 1): {"b": 1000, "f_max": 100},
(1, 2): {"b": 1000, "f_max": 100},
(0, 2): {"b": 1000, "f_max": 30},
},
}
triangle_fixed = OPF_switching(triangle, [(0, 2)], max_open=0)
triangle_switched = OPF_switching(triangle, [(0, 2)], max_open=1)
triangle_costs = [solve_switching(m) for m in [triangle_fixed, triangle_switched]]
pd.DataFrame({
"case": ["All lines active", "At most one open line"],
"dispatch cost": triangle_costs,
"cheap generation": [m.p[0]() for m in [triangle_fixed, triangle_switched]],
"costly generation": [m.p[2]() for m in [triangle_fixed, triangle_switched]],
"line (0, 2) active": [round(m.z[0, 2]()) for m in [triangle_fixed, triangle_switched]],
})
| case | dispatch cost | cheap generation | costly generation | line (0, 2) active | |
|---|---|---|---|---|---|
| 0 | All lines active | 3200.0 | 45.0 | 55.0 | 1 |
| 1 | At most one open line | 1000.0 | 100.0 | 0.0 | 0 |
Why opening a line can help#
With all three lines active, the direct line \((0,2)\) carries two-thirds of the cheap generator’s output. Its capacity of 30 therefore limits that output to 45. The remaining 55 units come from the costly generator, giving a dispatch cost of 3,200. In this small example, bus 2 hosts both demand and a generator; the same nodal balance handles their net injection. Opening the direct line sends all 100 units through the path \(0\to1\to2\), reducing the cost to 1,000. Its endpoints then have different angles: deleting the line’s angle equation is essential.
This is a stylized DC-flow example of the interaction between network topology and congestion. It does not establish that opening lines is generally beneficial, nor assess reliability under outages. It uses the same cost convention as the preceding objective: no duration factor or switching charge has been added.
Using the 118-bus data#
We apply the variant to the first 15-minute instance. A spanning tree stays active, and the first three non-tree lines in sorted order are eligible for switching. This small deterministic candidate set makes the example quick to run and independently check; it is not a search for the best switching candidates across all 179 lines. We permit at most one of these three lines to be open.
With zero allowed openings this formulation reproduces OPF1 for the same instance. Allowing openings includes the all-active solution, so its optimal cost cannot increase. Any saving reported here applies only to this instance and candidate set. The earlier 96-instance policy comparison remains unchanged.
The saving here is only about 0.42 objective units (0.00165% of the all-active cost), so it is negligible for this candidate set. This difference is smaller than HiGHS’s default relative MIP-gap allowance. The separate SWITCHING_SOLVER therefore requests zero relative and absolute MIP gaps; ordinary floating-point feasibility tolerances still apply.
# Keep a spanning tree active and allow three deterministic non-tree candidates.
grid = nx.Graph()
grid.add_nodes_from(network["nodes"])
grid.add_edges_from(network["edges"])
tree_edges = {frozenset(e) for e in nx.minimum_spanning_tree(grid).edges}
switchable_lines = sorted(
e for e in network["edges"] if frozenset(e) not in tree_edges
)[:3]
fixed_network = OPF_switching(network, switchable_lines, max_open=0)
switched_network = OPF_switching(network, switchable_lines, max_open=1)
fixed_cost = solve_switching(fixed_network)
switched_cost = solve_switching(switched_network)
assert abs(fixed_cost - sol1[0]) <= 1e-5
assert switched_cost <= fixed_cost + 1e-5
print("Switchable lines:", switchable_lines)
print(f"All lines active: {fixed_cost:.2f}")
print(f"At most one candidate open: {switched_cost:.2f}")
print("Open lines:", [e for e in switched_network.S if switched_network.z[e]() < 0.5])
saving = fixed_cost - switched_cost
print(f"Dispatch-cost saving: {saving:.4f} ({100 * saving / fixed_cost:.5f}%)")
Switchable lines: [(2, 11), (4, 10), (6, 11)]
All lines active: 25556.39
At most one candidate open: 25555.97
Open lines: [(4, 10)]
Dispatch-cost saving: 0.4205 (0.00165%)