Graph Convolutional Neural Networks for Optimal Power Flow Locational
  Marginal Price
The real-time electricity market with the integration of renewable energies
and electric vehicles have been receiving significant attention recently. So
far most of the literature addresses the optimal power flow (OPF) problem in
the real-time electricity market context by iterative methods. However, solving
OPF problems in real-time is challenging due to the high computational
complexity by the iterative methods. Motivated by this fact, in this paper, we
propose a Chebyshev Graph Convolutional Neural Networks (ChebGCN) to improve
the efficiency of integrating low-carbon energy sources into power grids and to
address scalability and adaptivity of end-to-end existing OPF solutions. The
proposed GCN method is capable to predict the optimal energy market marginal
prices in real time. Numerical analysis is used to benchmark the results and
validate the improvement.
Adrian-Petru Surani, Rahul Sahetiya
# Graph Convolutional Neural Networks for Optimal Power Flow Locational Marginal Price

###### Abstract

The real-time electricity market with the integration of renewable energies and electric vehicles have been receiving significant attention recently. So far most of the literature addresses the optimal power flow (OPF) problem in the real-time electricity market context by iterative methods. However, solving OPF problems in real-time is challenging due to the high computational complexity by the iterative methods. Motivated by this fact, in this paper, we propose a Chebyshev Graph Convolutional Neural Networks (ChebGCN) to improve the efficiency of integrating low-carbon energy sources into power grids and to address scalability and adaptivity of end-to-end existing OPF solutions. The proposed GCN method is capable to predict the optimal energy market marginal prices in real time. Numerical analysis is used to benchmark the results and validate the improvement.

## I Introduction

Electricity market pricing is one of the core tasks of operating power-grids, as the real-time market determines the incremental adjustment for the day-ahead by solving the optimal-power flow (OPF) problem [1]. The real-time OPF for market pricing also ensures high efficiency and reliability of grid operations, especially in the case of low-carbon energy sources integration, including renewable energy generations and electric vehicle integration. As the amount of smart, flexible resources' participation in distribution system operations keeps increasing, market-based approaches are needed for economic efficiency. [2] A key component to enable power distribution markets and provide economic incentives to market participants is the distribution locational marginal price (DLMP).

Typically, solving OPF problems by iterative methods incurs excessive computational complexity, limiting its applicability in large-scale power networks and real-time implementations. The locational marginal price (LMP) represents the marginal cost to serve one incremental unit of demand at a specific location in electric power networks, [3, 4]. Typically the LMP is calculated by solving the direct current optimal power flow (DC OPF) problem. Due to high voltage volatility and high power losses due to high R/X ratio, reactive power modeling is needed. As the DC OPF approach is unable to capture power losses and reactive powers in distribution systems, it is unsuitable for calculating the DLMP. A different approach has been proposed to calculate DLMPs in [2, 5, 6] which uses the single-phase alternating current optimal power flow (AC OPF) problem for the balanced system. At the same time, as LMPs relate to the duality analysis for OPF, their dependence on grid topology has been recognized in [7, 8].

To improve the high computation complexity of the accurate AC OPF problem, which is non-convex and non-linear, machine learning (ML) techniques have been proposed throughout existing literature: identifying active constraints [9], finding a warm start for iterative OPF solutions [10], or addressing the feasibility issue [11]. These models require an extensive off-line training of neural network (NN) models and as they rely on end-to-end NNs, also incur a high model and computation complexity for large-scale power grids. An additional cost is represented by the re-training of the models whenever the system inputs change as a result of frequently varying topologies.

To address these problems, we are motivated to utilize graph convolutional neural networks (GCN) to approximate the optimal marginal prices. The proposed method considers the power system measurements as the low-pass graph signals, and derive the suitable Graph Shift Operator (GSO) to design GCN. The proposed method also designs the regulation terms for the feasibility of power flow constraints.

The rest of the paper is organized as follows: The AC OPF formulation is introduced initially in Section II. In Section III, the proposed approach is discussed including various adaptations and optimizations such as convolution in Fourier domain to speed up computations. The case study results on a transmission system are given in Section IV, while Section V concludes the paper and discusses future research directions and exploratory questions.

## II Problem Formulation

The LMP at a certain location represents the incremental cost to supply an extra unit of load at this location. The typical approach is to use the DC OPF problem to determine the LMPs in the transmission level. The formulation of the considered AC OPF problem is:
\[\begin{array}{ll}\min\limits_{P_{G,k}}&\sum\limits_{k\in\mathcal{N}_{G}}C_{i}(P_ {G,k})\\ \text{s.t.}&P_{G,i}-P_{L,i}=\Re\{V_{i}\cdot(I_{i})^{*}\},i\in\mathcal{N}\\ &Q_{G,i}-Q_{L,i}=\Im\left\{V_{i}\cdot(I_{i})^{*}\right\},i\in\mathcal{N}\\ &\underline{P}_{G,k}\leq P_{G,k}\leq\overline{P}_{G,k},k\in\mathcal{N}_{G}\\ &\underline{Q}_{G,k}\leq Q_{G,k}\leq\overline{Q}_{G,k},k\in\mathcal{N}_{G}\\ &\underline{V_{i}}\leq|V_{i}|\leq\overline{V_{i}},i\in\mathcal{N}\end{array} \tag{1}\]

where \(\mathcal{N}_{G}\) denotes the set of buses which have generators. \(\Re\{\cdot\}\) and \(\Im\{\cdot\}\) denote the real and the imaginary part of a complex number, respectively. \(P_{G,k}\) and \(Q_{G,k}\) represent the active and reactive power output of the generator at bus k. Similarly, \(\underline{P}_{G,k},\underline{Q}_{G,k}\) and \(\overline{P}_{G,k},\overline{Q}_{G,k}\) correspond to the lower bounds and upper bounds for the active and reactive power generation. \(|V_{i}|\) corresponds to the voltage magnitude at bus \(i\), and \(\underline{V_{i}},\overline{V_{i}}\) the associated lower and upper bounds.

For simplicity, let \(u\) denote the vector of all control variables and \(x\) the vector of all state variables, including the voltage magnitude and angle at every bus and phase. We define \(h(x,u)\leq 0\) to represent the inequality constraints in (1). The Lagrangian function of the formulated AC OPF problem can be written as:

\[\begin{array}{ll}&L(x,u,\lambda,\nu,\mu)=\sum\limits_{k\in\mathcal{N}_{G}}C _{i}(P_{G,k})\\ &-\sum_{i\in\mathcal{N}}\lambda_{i}\left(P_{G,i}-P_{L,i}-\Re\{V_{i}\cdot(I_{i })^{*}\}\right)\\ &-\sum_{i\in\mathcal{N}}\nu_{i}\cdot(Q_{G,i}-Q_{L,i}-\Im\{V_{i}\cdot(I_{i}) ^{*}\}^{*})\\ &+\sum\limits_{m\in\mathcal{H}}\mu_{m}\cdot h_{m}(x,u)\end{array} \tag{2}\]

where \(\mathcal{N}\) represents the set of all buses in the system, \(\mathcal{P}_{i}\) the set of phases at bus i, and \(\mathcal{H}\) the set of inequality constraints. \(\lambda_{i}\) and \(\nu_{i}\) are the Lagrange multipliers corresponding to the active power balance equation, and the reactive power balance equation. \(u_{m}\) is the Lagrange multiplier associated with the inequality constraint \(h_{m}(x,u)\leq 0\).

Assuming the considered AC OPF problem has an optimal solution \((x^{*},u^{*})\), it can be calculated as:

\[DLMP_{i}=\left.\frac{\partial f}{\partial P_{L,i}}\right|_{x^{*},u^{*}}=\left. \frac{\partial L}{\partial P_{L,i}}\right|_{x^{*},u^{*}}=\lambda_{i} \tag{3}\]

Note that each bus has a marginal price \(\lambda_{i}\).

## III Proposed Approach

The GCNs are used to learn a mapping function from voltage measurements to optimal marginal prices.

### _Chebyshev Graph Neural Networks_

In spectral graph analysis, the graph Laplacian is an essential operator for undirected connected graphs, defined as: \(\mathbf{L}=\mathbf{D}-\mathbf{A}\in\mathbb{R}^{N\times N}\), while the normalized Laplacian is \(\mathbf{L}=\mathbf{I}_{N}-\mathbf{D}^{-1/2}\mathbf{A}\mathbf{D}^{-1/2}\), where \(\mathbf{I}_{N}\) is the identity matrix, \(\mathbf{A}\) is the adjacent matrix, and \(\mathbf{D}\in\mathbb{R}^{N\times N}\) is the diagonal degree matrix with \(\mathbf{D}_{ii}=\sum_{j}\mathbf{A}_{ij}\). The Laplacian matrix is diagonalized by the Fourier basis \(U=[u_{0},...,u_{N-1}]\in\mathbb{R}^{N\times N}\), where \(\{u_{l}\}_{l=0}^{\mathbb{N}-1}\in\mathbb{R}^{N}\) are the orthogonal eigenvectors of \(\mathbf{L}\). The eigenvalue decomposition for the Laplacian matrix is \(L=\mathbf{U}\mathbf{\Lambda}U^{T}\), where \(\Lambda=diag([\lambda_{0},...\lambda_{N-1}])\in\mathbb{R}^{N\times N}\) and \(\{\lambda_{l}\}_{l=0}^{\mathcal{N}-1}\in\mathbb{R}^{\mathcal{N}}\) is the ordered list of non-negative real eigenvalues associated with \(\mathbf{L}\), referred to as the graph's frequencies. The voltage phasor at time \(t\) is represented as the signal \(x=\left[|\mathbf{V}|,\mathbf{\theta}\right]\in\mathbb{R}^{2N}\) and the graph Fourier transform of the signal \(x\in\mathbb{R}^{2N}\) is \(\hat{x}=\mathbf{U}^{T}x\in\mathbb{R}^{2N}\), and its inverse Fourier transform as \(x=\mathbf{U}\hat{x}\). In the Fourier domain, the convolution on graph \(\mathcal{G}\) is denoted \(x_{*\mathcal{G}y}=U((U^{T}x)\odot(U^{T}y))\), where the \(\odot\) represents the element-wise product. Filtering the signal \(x\) by \(g_{\theta}\):

\[y=g_{\theta}(\mathbf{L})x=g_{\theta}(\mathbf{U}\Lambda\mathbf{U}^{T})x=\mathbf{U}g_{\theta}( \Lambda)\mathbf{U}^{T}x \tag{4}\]

The ChebNet model uses mainly Chebyshev polynomials to replace the convolution kernel in the spectral domain to represent the filter, as defined in 5, which is the truncated expansion of order \(K-1\):

\[g_{\theta}(\Lambda)=\sum\limits_{k=0}^{K-1}\theta_{k}T_{k}(\tilde{\Lambda}) \tag{5}\]

In (5), the parameter \(\theta\in\mathbb{R}^{K}\) is the vector of Chebyshev coefficients. \(T_{k}(\tilde{\mathbf{\Lambda}})\in\mathbb{R}^{\mathcal{N}\times\mathcal{N}}\) is the Chebyshev polynomial with order \(k\) evaluated at \(\tilde{\mathbf{\Lambda}}=2\lambda/\lambda_{max}-\mathbf{I}_{\mathcal{N}}\), where \(\tilde{\Lambda}_{ii}\in[-1,1]\) and \(\lambda_{max}\) is the maximum eigenvalue.

The calculation of the Chebyshev polynomial, \(T_{k}(x)\) can be obtained through the following recursive equation:

\[\left\{\begin{array}{ll}T_{0}(x)=1\\ T_{1}(x)=x\\ T_{\mathcal{N}+1}(x)=2xT_{\mathcal{N}}(x)-T_{\mathcal{N}-1}(x)\end{array}\right. \tag{6}\]

Therefore, the convolution operation is expressed as:

\[x_{*\mathcal{G}g_{\theta}}=\mathbf{U}(\sum\limits_{i=0}^{K-1}\theta_{k}T_{k}( \tilde{\mathbf{\Lambda}}))\mathbf{U}^{T}x \tag{7}\]

Similar to \(\tilde{\mathbf{\Lambda}}=2\mathbf{\Lambda}/\lambda_{max}-\mathbf{I}_{\mathcal{N}}\) and \(T_{k}(\tilde{\mathbf{L}})\), \(\tilde{\mathbf{L}}\), and \(T_{i}(\tilde{\mathbf{L}})\) can be defined as:

\[\tilde{\mathbf{L}}=2\mathbf{L}/\lambda_{max}-\mathbf{I}_{\mathcal{N}} \tag{8}\]

and

\[\begin{array}{ll}T_{k}(\tilde{\mathbf{L}})&=T_{k}(2\mathbf{L}/\lambda_{max}-\mathbf{I}_{ \mathcal{N}})\\ &=T_{k}(2\mathbf{U}\mathbf{\Lambda}\mathbf{U}^{T}/\lambda_{max}-\mathbf{U}\mathbf{I}_{\mathcal{N}}U ^{T})\\ &=T_{k}(\mathbf{U}(2\mathbf{\Lambda}/\lambda_{max}-\mathbf{I}_{\mathcal{N}})U^{T})\\ &=\mathbf{U}\mathbf{U}T_{k}(\hat{\mathbf{\Lambda}})\mathbf{U}^{T}\end{array} \tag{9}\]

By performing the Fourier transform, doing the convolution and then recovering the \(\mathbf{U}\) through the inverse Fourier transform, we can avoid extra-calculations, and simplify the equation further:

\[\begin{array}{ll}x_{*\mathcal{G}g_{\theta}}&=\mathbf{U}\left(\sum\limits_{i=0}^{K- 1}\theta_{k}T_{k}(\tilde{\mathbf{\Lambda}})\right)\mathbf{U}^{T}x\\ &=\sum\limits_{i=0}^{K-1}\theta_{k}\left(\mathbf{U}T_{k}(\tilde{\mathbf{\Lambda}})\mathbf{U}^ {T}\right)x\\ &=\sum\limits_{i=0}^{K-1}\theta_{k}T_{k}\left(\tilde{\mathbf{L}}\right)x\end{array} \tag{10}\]
As pointed above in (10), the calculation does not require eigenvector decomposition, which improves the computational complexity. Given a sample \(x\), we calculate the output output feature map denoted as \(y\), using the \(\theta_{i,j}\in\mathbb{R}^{K}\) as Chebyshev coefficients.

\[y=\sum_{k=1}^{K}g_{\theta_{k}}(\mathbf{L})x \tag{11}\]

where the output can be multiple channels.

### _GCN for Marginal Price Prediction and Forecasting_

Inspired by the graph signal \(x\) in (10), we propose to utilize Chebyshev GCN to predict the marginal prices, i.e., \(\lambda_{i}\) of (3), instead of solving a nonconvex optimization problem. In our problem, our input graph signals are the voltage magnitude and voltage angles. Our graph shift operator \(\mathbf{L}\) is the absolute value of system matrix \(\mathbf{Y}\), i.e., \(\mathbf{L}=|\mathbf{Y}|\). The process of training and testing processes can be summarized as follows.

1. Utilize the Interior Point Method (IPM) to solve the Lagrangian problem in (2) to obtain the optimal \(\lambda_{i}^{*},i\in\mathcal{N}\).
2. Train the graph neural networks by regarding \(|V_{i}|,i\in\mathcal{N}\) and \(|\theta_{i}|,i\in\mathcal{N}\) as inputs, and \(\lambda_{i},i\in\mathcal{N}\) as the outputs.
3. After training the GCN, we can predict the unseen \(\lambda_{i},i\in\mathcal{N}\) by putting the new-arrival measurements \(|V_{i}|,i\in\mathcal{N}\) and \(|\theta_{i}|,i\in\mathcal{N}\).

This framework can be also implemented for the forecasting problem if the power demands are time-series measurements. The only difference is Step 2, i.e., regarding \(|V_{i}(t)|,i\in\mathcal{N}\) and \(|\theta_{i}(t)|,i\in\mathcal{N}\) as inputs, and the future \(\lambda_{i}(t+1),i\in\mathcal{N}\) as the outputs.

## IV Analysis and numerical results

We will consider the IEEE 118-bus standard testing cases with the real-world renewable energies and EV charging powers for analysis. The 118-bus system has been used to represent the real-world demand, and for each demand vector the real and imaginary parts were separated.

### _System Setting_

The code was run on a machine powered by the Apple M1 chip, using 16 GB RAM. Most of the computations were run under an upper bound of 2 hours. The Python environment was built using the Conda package system and Python version 3.10.8. Pytorch 1.13.0 was used to build the GNN models and to subsequently benchmark them. We consider two benchmarks, i.e., fully connected neural networks (FCNN) and 1st-order approximated graph neural network [12].

### _Prediction of Marginal Prices_

In order to predict marginal prices, three approaches were used based on PyTorch's classic models: Fully-Connected Neural Networks (FCNNs), GNN (Graph Neural Networks) and Chebyshev GCN. Figure 1 compares the performances of the three methods and subsequently the improvements throughout the stages. In particular, Chebyshev GCN predicts the optimal marginal prices with \(4.5717e^{-5}\) MSEs, while GNN has 0.0016 and FCNN has \(2.3304e^{-4}\) MSEs. The GNN performance has been significantly better in both prediction and forecasting. We also illustrate the performance of FCNN, GNN, and Chebyshev GCN in Figs. 2, 3 and 4. It shows that Chebyshev GCN can predict the ground-truth marginal prices with very high accuracy.
show that ChebyGCN performs slightly better than FCNN. This is because forecasting is a time-series problem, while ChebyGCN does not consider the temporal correlations. We also illustrate three examples of FCNN, GNN and Chebyshev GCN for the forecasting marginal prices in Figs. 6-8. The Chebyshev GCN still produces a promising result that approximates the ground-truth marginal prices.
obtained by considering time-series data, which would require Recurrent Neural Networks (RNNs) and Long Short-Term Memory (LSTM).
