Modern power grids are vast, interconnected systems composed of thousands of generators, transmission assets, and load centers operating in real time. System operators are tasked with maintaining a delicate balance between electricity supply and demand at every moment, since even small imbalances can cause deviations in system frequency, risking equipment damage or large-scale blackouts. To prevent such instability, operators must continuously clear electricity markets and dispatch resources in a way that satisfies physical, economic, and regulatory constraints. As renewables, battery storage, and electric vehicle charging play larger roles in the grid, managing resource dispatch also becomes an increasingly sophisticated challenge.

A naive first attempt at solving the problem might involve getting a generating capacity and price (in $/MW) for every power plant, and assigning loads to plants in a way which minimizes total cost. However, this neglects several important constraints. One is that transmission lines have capacity limits, so there’s no point generating power if it can’t be brought to the consumer. Another is that power plants don’t work like the dimmer switch on a lightbulb — turning them on and off have costs, and power plants can’t ramp between generation loads instantly. Telling plant operators when they need to turn on and how fast to ramp ahead of time is known as the unit commitment problem. This is the challenge we’re focused on solving with machine learning techniques.

Unit Commitment

But first, what makes the unit commitment problem so difficult to solve efficiently? The unit commitment constraint effectively functions as a binary mask on top of the base generation capacity problem, so that for every commitment mask (of what power plants are on/off), the optimal generation for each plant must be solved as a subproblem. This problem formulation is nonconvex, which means that local minima are not necessarily global minima. When the grid size grows to hundreds of potential plants to turn on or off, the number of possible masks to solve grows combinatorially, making exhaustive search impossible. Instead, one of two options is typically used: mixed-integer solvers or heuristic commitment.

Mixed integer solvers include the binary commitment mask as an optimizable variable, which must be treated as an integer constraint. This problem formulation is solved by many commercially available codes, which typically use a “branch and bound” algorithm. Describing branch and bound is a bit beyond the scope of a single post, so refer to a resource such as this for more information. [1] Practically, mixed integer solvers offer optimality guarantees and much better performance compared to exhaustive search, but still scale poorly with increasing decision variable count and are significantly slower than solving the equivalent problem without the binary constraint. A sample problem formulation for the unit commitment problem as a mixed integer program is shown below:

The mixed-integer formulation is typically suitable for grid operators in real-time operation, but the computational cost begins to grow when running forward-looking simulations that are essential for grid planning and determining the economic value of future projects. Currently, one of the industry-standard methods for solving that problem is to use a heuristic commitment map. Unlike in the mixed integer program, the unit commitment is solved prior to the optimization with a heuristic that looks forward 24–72 hours. With the unit commitment constraint set by the heuristic, the problem becomes convex and can be solved rapidly.

Given the relative simplicity of the second method combining a heuristic with a convex linear program solve, it is useful to look at the heuristic in more detail. A good heuristic should take in the relevant problem constraints and weights (such as capacity limits, startup costs, projected demand), and output a binary commitment map which produces both a feasible and near-optimal solution (compared to the more expensive mixed-integer program). For applications which require a large number of evaluations such as grid planning, the heuristic must also evaluate rapidly and scale well with increasing grid complexity.

Graph Neural Networks

For this, we turned to graph neural networks (GNNs). The structure of the transmission lines, generators, and loads within a power grid take the form of a graph. Loads and generators are located at buses which serve as nodes in the graph, and transmission lines, represented by edges, connect the buses. Graph neural networks are machine learning architectures uniquely suited to take advantage of the connectivity of graph structures to perform learned inference tasks at the level of individual nodes and edges or the graph as a whole. The complexity of a GNN can be linear in the number of edges and inductive over changes in graph structure. These properties allow for training data to be drawn from a variety of sample power grid architectures with favorable scaling to realistic large power grids.

Graph neural networks have had past successes in power transmission problems. Optimal power flow (OPF) is a classic problem in transmission grid operation, describing where to optimally generate power given a realistic transmission grid and load distributed over that grid. Solving the OPF problem is prohibitively expensive in real time, so grid operators typically solve simplified versions with associated costs in accuracy. CANOS (Constraint-Augmented Neural OPF Solver) is a graph neural network architecture which generates near-optimal solutions while scaling well to utility-scale grid sizes. [2] Additionally, GNNs have been applied to graph-level tasks over power grids, such as predicting whether a given grid operating condition is susceptible to cascading failures. [3]

In the case of the unit commitment heuristic, we want the model to predict whether a given generator is committed at a specific time, given the same information which is provided to the mixed integer optimization problem. This is a node-level inference problem, where the GNN predicts a binary classification denoting the commitment status of a generator node. Before we start looking at GNN architectures, let’s first check out the training data we’ll use for the problem.

UnitCommitment.jl and PSTCA

The data we’ll be using for training will come from the UnitCommitment.jl package. [4] The package compiles benchmark instances of the Unit Commitment problem, as well as a formulation for solving them with mixed integer programming techniques. Specifically, the package uses the security constrained unit commitment formulation, which attempts to find the optimal unit commitment which are robust to a set of contingencies, which dictate changes to the network architecture such as downed transmission lines. These contingencies will be useful later for augmenting our dataset. For the ground truth unit commitment solutions, we will use the built-in tools in UnitCommitment.jl.

The specific dataset we are using is the Power Systems Test Case Archive (PSTCA), compiled by the University of Washington and translated into the standard UnitCommitment.jl json data format. [5] These cases represent portions of the grid maintained by the American Electric Power System in the 1960s, and range in size from networks with 14 buses up to 300 buses. There is one variation for each day of a sample year, with demand projected at each bus over the next 36 hours starting from midnight of the given day. Take a look below for a snippet of a json file describing a simple 14-bus case.

We see in the header some information about the case setup:

"SOURCE": "- University of Washington, Dept. of Electrical Engineering. \"Power Systems Test Case Archive\". Available at: http://www.ee.washington.edu/research/pstca/ (Accessed: Nov 14, 2020)\n\n- D. Zimmerman, C. E. Murillo-Sandnchez and R. J. Thomas. \"Matpower: Steady-state operations, planning, and analysis tools forpower systems research and education\", IEEE Transactions on PowerSystems, vol. 26, no. 1, pp. 12 \u201319, Feb. 2011. DOI: 10.1109/TPWRS.2010.2051168\n\n- A. S. Xavier, A. M. Kazachkov, O. Yurdakul, F. Qiu. \"UnitCommitment.jl: A Julia/JuMP Optimization Package for Security-Constrained Unit Commitment (Version 0.3)\". Zenodo (2022). DOI: 10.5281/zenodo.4269874\n",
"Parameters": {
  "Version": "0.3",
  "Power balance penalty ($/MW)": 1000.0,
  "Time horizon (h)": 36
}

Then, the generators in the case are listed:

"Generators": {
    "g1": {
      "Bus": "b1",
      "Production cost curve (MW)": [
        36.75123365792352,
        110.10634590601819,
        183.46145815411288,
        256.8165704022076,
        330.1716826503023
      ],
      "Production cost curve ($)": [
        1160.1370560288465,
        4051.592909653724,
        7310.413973593571,
        10614.711162105745,
        14293.941140664052
      ],
      "Startup costs ($)": [
        25056.51,
        29059.15,
        34238.45
      ],
      "Startup delays (h)": [
        1,
        2,
        4
      ],
      "Ramp up limit (MW)": 231.12,
      "Ramp down limit (MW)": 231.12,
      "Startup limit (MW)": 231.12,
      "Shutdown limit (MW)": 231.12,
      "Minimum uptime (h)": 1,
      "Minimum downtime (h)": 1,
      "Reserve eligibility": [
        "r1"
      ],
      "Initial status (h)": 24,
      "Initial power (MW)": 230.74887999999999
    },
# etc.
}

A piecewise linear cost curve is used to interpolate between electricity production amount and generation cost. Similarly, leaving a power plant shut off for longer periods could incur larger costs when restarting, so another piecewise linear cost curve is used to represent the cost of starting the plant back up after being shut down for a given number of hours.

Finally, we can look at the transmission lines, which connect pairs of buses, and projected loads on buses over the next 36 hours:

"Transmission lines": {
  "l1": {
    "Source bus": "b1",
    "Target bus": "b2",
    "Reactance (ohms)": 0.05917,
    "Susceptance (S)": 29.49686
  },
#etc.
}
"Buses": {
    "b1": {
      "Load (MW)": 0.0
    },
    "b2": {
      "Load (MW)": [
        19.33301,
        18.57311,
        18.06675,
        17.76778,
        17.95852,
        18.33015,
        19.37086,
        20.13001,
        #etc.
      ]
    },
#etc.
}

Next up, we need to transform this data into a format usable for our GNN.

Generating Training Data

To get the ground truth solution from a given UnitCommitment.jl case, we can solve the problem with the built-in `optimize()` method. This generates a binary mask for each of the 36 lookahead hours given projected demand as the solution, which we will then use to train the GNN. For those interested, we are using the HiGHS mixed integer solver within the JuMP framework, and the solution strategy is described in a 2019 paper by Xavier et. al. [6]

using HiGHS
using UnitCommitment

# load a benchmark case
instance = UnitCommitment.read_benchmark("matpower/case14/2017-01-01")

# build optimization model
model = UnitCommitment.build_model(instance = instance, optimizer = HiGHS.Optimizer)

# solve the built model
UnitCommitment.optimize!(model)

# extract the solution
solution = UnitCommitment.solution(model)

The solution is printed out to a json file, which contains the variable “is on” that we care about , which is our binary unit commitment.

"Is on": {
  "g1": [
    1.0,
    1.0,
    1.0,
    1.0,
    1.0,
# etc.
  ],
}

Let’s take a look at a sample input and solution graph. In the figure below, the generator buses are represented in orange and the load buses in blue. Each bus is connected to others by transmission lines which are represented by a lumped element resistance and susceptance. Additional node-level information is not depicted, such as startup costs and ramping rates for the generation assets.

A solution graph is depicted below, corresponding to the unit commitment for the case above. This serves as a binary mask over all of the generator nodes, which becomes the optimization objective for the GNN. Note that only buses containing generator nodes are affected by the commitment status.

However, so far we’ve only looked at the 4 available power grid architectures within the benchmark instances (14, 57, 118, and 300 bus networks) of which we only have 365 examples each. To make sure our GNN can adapt to unseen transmission grid structures, we take advantage of the contingencies. UnitCommitment.jl contains a set of available contingencies for each transmission network case, which consists of all of the single transmission line removals which can occur without making the transmission network a disjoint graph. We sample five random contingencies for each data point, solve them, and add them to the dataset for training.

Now that we have solved cases, these must be converted into a format compatible with PyTorch Geometric, the machine learning package we plan on using to implement the GNNs. We use the following code to do so. First, we encode generator nodes and connect them to buses:

generators = data["Generators"]
buses = data["Buses"]
lines = data["Transmission lines"]

data_hetero = HeteroData()

# Build generator features
gen_features = []
gen_bus_idx = []
for g_name, g in generators.items():
    if len(g.get("Startup costs ($)", [])) < 3:
        g["Startup costs ($)"] = prev_startup_costs.get(g_name, [0.0]*3)
    else:
        prev_startup_costs[g_name] = g["Startup costs ($)"]

    if len(g.get("Startup delays (h)", [])) < 3:
        g["Startup delays (h)"] = prev_startup_delays.get(g_name, [0.0]*3)
    else:
        prev_startup_delays[g_name] = g["Startup delays (h)"]

    features = (
        g["Production cost curve (MW)"]
        + g["Production cost curve ($)"]
        + g["Startup costs ($)"]
        + g["Startup delays (h)"]
        + [ 
            g["Ramp up limit (MW)"],
            g["Ramp down limit (MW)"],
            g["Startup limit (MW)"],
            g["Shutdown limit (MW)"],
            g["Minimum uptime (h)"],
            g["Minimum downtime (h)"],
            g["Initial status (h)"],
            g["Initial power (MW)"]
        ]
    )
    gen_features.append(features)
    gen_bus_idx.append(bus2idx[g["Bus"]])

gen_features = torch.tensor(gen_features, dtype=torch.float)
gen_features = (gen_features - gen_features.mean(dim=0)) / (gen_features.std(dim=0) + 1e-6)

data_hetero["generator"].x = gen_features
src = torch.arange(len(gen_bus_idx))
dst = torch.tensor(gen_bus_idx)
data_hetero["generator", "produces_at", "bus"].edge_index = torch.stack([src, dst], dim=0)
data_hetero["bus", "served_by", "generator"].edge_index = torch.stack([dst, src], dim=0)

Then, we add in bus features for load and transmission:

bus_features = []
for bus_data in buses.values():
    load_seq = bus_data["Load (MW)"]
    if isinstance(load_seq, (int, float)):
        load_seq = [load_seq] * lookahead_hours
    elif len(load_seq) < lookahead_hours:
        load_seq = list(load_seq) + [0.0]*(lookahead_hours - len(load_seq))
    bus_features.append(load_seq)

bus_features = torch.tensor(bus_features, dtype=torch.float)
data_hetero["bus"].x = bus_features

# Build transmission edges, removing failed lines if specified
bus_edge_index = []
bus_edge_attr = []
failed_line_set = set(failed_lines) if failed_lines else set()

for line_id, props in lines.items():
    # Skip failed lines in contingency scenarios
    if line_id in failed_line_set:
        continue
        
    src = bus2idx[props["Source bus"]]
    dst = bus2idx[props["Target bus"]]
    bus_edge_index.append([src, dst])
    bus_edge_index.append([dst, src])
    bus_edge_attr.append([props["Reactance (ohms)"], props["Susceptance (S)"]])
    bus_edge_attr.append([props["Reactance (ohms)"], props["Susceptance (S)"]])

if bus_edge_index:  # Only add if edges exist
    bus_edge_index = torch.tensor(bus_edge_index, dtype=torch.long).T
    bus_edge_attr = torch.tensor(bus_edge_attr, dtype=torch.float)
    data_hetero["bus", "transmission", "bus"].edge_index = bus_edge_index
    data_hetero["bus", "transmission", "bus"].edge_attr = bus_edge_attr

Finally, we add in reserve nodes. We haven’t talked about them much, but they essentially serve as an additional requirement for generation capacity across the entire grid which is ready to go, forecasted out over the next 36 hours.

reserve_features = []
reserve_names = list(data.get("Reserves", {}).keys())
reserve_name2idx = {r: i for i, r in enumerate(reserve_names)}

for r_name in reserve_names:
    r_data = data["Reserves"][r_name]
    amount_seq = r_data["Amount (MW)"]
    if len(amount_seq) < lookahead_hours:
        amount_seq = list(amount_seq) + [0.0]*(lookahead_hours - len(amount_seq))
    reserve_features.append(amount_seq)

if reserve_features:
    reserve_features = torch.tensor(reserve_features, dtype=torch.float)
    data_hetero["reserve"].x = reserve_features

reserve_src, reserve_dst = [], []
for r_idx, r_name in enumerate(reserve_names):
    for g_idx, g in enumerate(generators.values()):
        if r_name in g.get("Reserve eligibility", []):
            reserve_src.append(r_idx)
            reserve_dst.append(g_idx)

if reserve_src:
    data_hetero["reserve", "backed_by", "generator"].edge_index = torch.tensor(
        [reserve_src, reserve_dst], dtype=torch.long
    )

Model Architecture

Our model can be decomposed into four distinct parts:

1. Load in MLP

First, we use two MLP layers to project the 24 generator features to a 128 dimensional space. Both the reserve and load node types have no feature information and thus ignore the load-in step. Instead they are initialized to random 128-dimensional vectors.

2. Structural GNN

The structural encoder is a two-layer HEATConv network. This model augments the Graph Attention Transformer model with edge features and heterogenous node types, all of which are necessary to realize the interactions between generators, nodes, and transmission lines. [7] Each HEAT layer performs heterogeneous, edge-aware 4-head attention over the power-system graph. The update rule in each layer resembles (skip-connection not included here but tacked on to the second layer):

Importantly, line parameters such as reactance and susceptance enter directly through edge-attribute embeddings.

Since power networks are highly interconnected, we discovered that deep GNNs (more than 3 layers) suffer from oversmoothing, where all node embeddings become indistinguishable. Using two HEAT layers performed best in the example systems, but additional layers may be acceptable in larger systems that are less interconnected. Below is a little code snippet showing an individual HEAT layer and the forward pass of the GNN.

self.heat1 = HEATConv(
            in_channels=hidden_dim, out_channels=hidden_dim,
            num_node_types=self.num_node_types, num_edge_types=self.num_edge_types,
            edge_type_emb_dim=edge_emb_dim, edge_dim=2, edge_attr_emb_dim=edge_emb_dim,
            heads=num_heads, concat=False
        )

x = self.heat1(x, edge_index, node_type, edge_type, edge_attr)
x1 = torch.relu(x)
x = self.heat2(x1, edge_index, node_type, edge_type, edge_attr)
x_out = torch.relu(x) + x1
x_out = self.dropout(x_out)

3. Temporal Layer

After structural message passing is completed by the two-layer HEAT encoder, temporal dynamics are modeled outside of the GNN using a lightweight convolutional sequence model (GRU was also used, but didn’t show any advantages). This design strictly separates structural learning from temporal learning, preventing any form of time leakage through message passing.

Once the GNN runs, load and reserve 36-hour trajectories are introduced through a learned linear projection. These projected sequences are concatenated to the corresponding static load and reserve embeddings to form temporal node representations. Generator nodes do not receive direct temporal inputs. Instead, their static embeddings are broadcast across time, reflecting the fact that their physical characteristics do not change over the horizon. The resulting temporal node sequences are processed by a two-layer 1D convolutional network. Below is the construction of the temporal layers and how we implement them:

if temporal_method == '1d_conv':
    self.temporal = nn.Sequential(
        nn.Conv1d(hidden_dim, hidden_dim, kernel_size=5, padding=2),
        nn.BatchNorm1d(hidden_dim),
        nn.ReLU(),
        nn.Dropout(dropout),
        nn.Conv1d(hidden_dim, hidden_dim, kernel_size=3, padding=1),
        nn.BatchNorm1d(hidden_dim),
        nn.ReLU()
    )

# Get structural embeddings from GNN
z_graph, _, offsets = self.encoder(x_dict_static, graph.edge_index_dict, graph.edge_attr_dict)

# Extract generator structural embeddings (only computed once!)
gen_emb_static = z_graph[offsets['generator']:offsets['generator']+n_gen]  # [n_gen, hidden_dim]
bus_emb_static = z_graph[offsets['bus']:offsets['bus']+n_bus]  # [n_bus, hidden_dim]
res_emb_static = z_graph[offsets['reserve']:offsets['reserve']+n_reserve]  # [n_reserve, hidden_dim]

# Build temporal sequences efficiently (vectorized over time)
bus_temporal = bus_feats.unsqueeze(-1)  # [n_bus, T, 1]
res_temporal = res_feats.unsqueeze(-1)  # [n_reserve, T, 1]

# Project temporal features
bus_temporal_emb = self.bus_temporal_proj(bus_temporal)  # [n_bus, T, hidden_dim]
res_temporal_emb = self.reserve_temporal_proj(res_temporal)  # [n_reserve, T, hidden_dim]

# Combine structural + temporal for bus and reserve
bus_emb_seq = bus_emb_static.unsqueeze(1) + bus_temporal_emb  # [n_bus, T, hidden_dim]
res_emb_seq = res_emb_static.unsqueeze(1) + res_temporal_emb  # [n_reserve, T, hidden_dim]

# Generators: expand static embedding across time (no temporal features for generators)
gen_emb_seq = gen_emb_static.unsqueeze(1).expand(-1, self.T, -1)  # [n_gen, T, hidden_dim]

# Concatenate all nodes
z_seq = torch.cat([gen_emb_seq, bus_emb_seq, res_emb_seq], dim=0)  # [total_nodes, T, hidden_dim]

# Apply temporal model
if self.temporal_method == '1d_conv':
    # Conv1d expects [batch, channels, length]
    z_seq_t = z_seq.permute(0, 2, 1)  # [total_nodes, hidden_dim, T]
    z_temporal = self.temporal(z_seq_t)  # [total_nodes, hidden_dim, T]
    z_temporal = z_temporal.permute(0, 2, 1)  # [total_nodes, T, hidden_dim]

4. Decoder

Finally, each temporal node embedding is decoded through a small MLP to produce 36-step unit-commitment. Only generator predictions are retained as outputs:

# Decoder
self.decoder = nn.Sequential(
    nn.Linear(hidden_dim, hidden_dim // 2),
    nn.ReLU(),
    nn.Dropout(dropout),
    nn.Linear(hidden_dim // 2, output_dim)
)

# Decode to predictions
out_seq = self.decoder(z_temporal)  # [total_nodes, T, 1]
out_seq = out_seq.permute(1, 0, 2)  # [T, total_nodes, 1]

# Extract only generator predictions
gen_predictions = out_seq[:, :n_gen, :]  # [T, n_gen, 1]

We use binary cross entropy loss as the metric to backpropagate for the accuracy of the unit commitment predictions. The full model flow is shown in the diagram below:

Results

Our main benchmark was a simple ridge regression. To do this we simply flattened the inputs (generator features, 36-hour load and reserve lookahead) and used them to predict the flattened unit commitments. In an effort to speed up training we reduced the input features to one hundred dimensions using an efficient PCA algorithm. However, this ridge model needed to be trained separately on each different case as it had abysmal accuracy between different grid set ups. As a result, we could not use it to provide a reasonable benchmark against a mix of systems, but this inability to adapt to new systems highlights the main advantage of our model: flexibility.

While these results may look lackluster at first glance, they are in fact extremely strong. It is important to remember that not only does the ridge model completely fail to adapt to new grid structures, it also cannot handle contingencies. Contingencies are when something goes wrong in the grid (a line or a generator goes out) and they happen all the time. A strong power model needs to be able to handle these issues and still produce valid unit commitment instructions. While the ridge can be reformulated to accommodate line contingencies by zero-padding, the resulting performance (not shown in the table above) hovers between 70–85% depending on the active sets of contingencies in the dataset. We expect that on larger grids with more unique outage combinations and complex dependencies, ML methods that are structure-agnostic (like ridge) will simply not be able to cut it.

Back to the table above, it is impressive that our generalized model–which can work on a variety of different grid structures and accommodate any combination of contingency/outage situations–performs almost on par with the ridge model which can only perform under strict fixed settings.

Next Steps

We believe that with a few small (albeit time consuming) improvements, our model could be used as a much faster alternative to the MIP or an ugly heuristic that is currently employed to solve the unit commitment problem. First, we need to get our code factorized to run on GPUs to save time. After factorizing the code properly, we expect training time to be drastically reduced. With this reduction, we could train the model on considerably larger grids, with many years of data and a wide variety of outage contingencies. By doing so, we believe that we could create a single model that could generalize to all sorts of different power grids. Lastly, we plan to compare our unit commitment output to real MIP unit commitment output on dispatch. What exactly does this mean? Once units are committed, they are then dispatched in real time so that supply meets demand. Only units that were committed in advance are allowed to be dispatched, hence the importance of getting unit commitment correct. Conveniently, this dispatch model is a simple LP and allows us to quickly compare whether our unit commitment model leads to system failures and energy costs on par with the traditional MIP and heuristic unit commitment models.

As highlighted throughout this post, our end goal is speed. Even if our model produces slightly worse results than the traditional MIP or heuristic models, we expect it to run orders of magnitude faster. Why is speed so important? While real-time operations place more emphasis on accuracy, speed is much more important for forward-looking studies. Each study contains hundreds of different weather, load, and generation scenarios over thousands of hours. That’s not even taking into account that each of these scenario-hours rely on random probabilistic outage draws and done via Monte Carlo simulation! As a result, the electrical-grid study space is largely bottlenecked by the MIP/heuristic unit commitment problem and is in desperate need of a faster alternative. Solving it via the model we proposed (if properly factored for GPUs) compared to the heuristic could conservatively increase the speed of each iteration by 50 times from 500 ms to somewhere around 10 ms. We hope that by continuing to build out this model together over the next quarter, we can provide the electricity industry with a much-needed speedy alternative.

Learn More

If you want to learn more about the project and play around with the tools and datasets, we have a Colab demonstrating the workflow here. We also host the code and training data on a Github repo here.

References

[1] Boyd, S., Ghosh, A., & Magnani, A. (2003, November 1). Branch and bound methods (Notes for EE392o, Stanford University). Stanford University. https://web.stanford.edu/class/ee392o/bb.pdf

[2] Piloto, L., Liguori, S., Madjiheurem, S., Zgubic, M., Lovett, S., Tomlinson, H., Elster, S., Apps, C., & Witherspoon, S. (2024). CANOS: A fast and scalable neural AC-OPF solver robust to N-1 perturbations (arXiv:2403.17660) [Preprint]. arXiv. https://arxiv.org/abs/2403.17660

[3] Varbella, A., Amara, K., Gjorgiev, B., El-Assady, M., & Sansavini, G. (2024). PowerGraph: A power grid benchmark dataset for graph neural networks (arXiv:2402.02827v2) [Preprint]. arXiv. https://arxiv.org/abs/2402.02827v2

[4] Xavier, A. S., Kazachkov, A. M., Yurdakul, O., He, J., & Qiu, F. (2025). UnitCommitment.jl: A Julia/JuMP optimization package for security-constrained unit commitment (Version 0.4.2) [Computer software]. Zenodo. https://doi.org/10.5281/zenodo.17737778

[5] University of Washington Electrical & Computer Engineering. (n.d.). Power Systems Test Case Archive [Data repository]. University of Washington. https://labs.ece.uw.edu/pstca/

[6] Xavier, Á. S., Qiu, F., Wang, F., & Thimmapuram, P. R. (2019). Transmission constraint filtering in large-scale security-constrained unit commitment. IEEE Transactions on Power Systems, 34(3), 2457–2460. https://doi.org/10.1109/TPWRS.2019.2892620

[7] Mo, X., Xing, Y., & Lv, C. (2021). Heterogeneous edge-enhanced graph attention network for multi-agent trajectory prediction (arXiv:2106.07161) [Preprint]. arXiv. https://arxiv.org/abs/2106.07161

Using Graph Neural Networks to Solve the Unit Commitment Problem was originally published in Stanford CS224W: Machine Learning with Graphs on Medium, where people are continuing the conversation by highlighting and responding to this story.