# Machine Learning for Electricity Market Clearing

Laurent Pagnier 1, Robert Ferrando 1, Yury Dvorkin 2, and Michael Chertkov 1
1 Program in Applied Mathematics, University of Arizona, Tucson, AZ, USA
<EMAIL_ADDRESS>2 Electrical and
Computer Engineering, New York University, New York, NY, USA
<EMAIL_ADDRESS>
###### Abstract

This paper seeks to design a machine learning twin of the optimal power flow
(OPF) optimization, which is used in market-clearing procedures by wholesale
electricity markets. The motivation for the proposed approach stems from the
need to obtain the digital twin, which is much faster than the original, while
also being sufficiently accurate and producing consistent generation
dispatches and locational marginal prices (LMPs), which are primal and dual
solutions of the OPF optimization, respectively. Availability of market-
clearing tools based on this approach will enable computationally tractable
evaluation of multiple dispatch scenarios under a given unit commitment.
Rather than direct solution of OPF, the Karush–Kuhn–Tucker (KKT) conditions
for the OPF problem in question may be written, and in parallel the LMPs of
generators and loads may be expressed in terms of the OPF Lagrangian
multipliers. Also, taking advantage of the practical fact that many of the
Lagrangian multipliers associated with lines will be zero (thermal limits are
not binding), we build and train an ML scheme which maps flexible resources
(loads and renewables) to the binding lines, and supplement it with an
efficient power-grid aware linear map to optimal dispatch and LMPs. The scheme
is validated and illustrated on IEEE models. We also report a trade of
analysis between quality of the reconstruction and number of samples needed to
train the model.

###### Index Terms:

Optimal Power Flow, Locational Marginal Prices, Physics-Informed Machine
Learning, Karush–Kuhn–Tucker Conditions

##  I Introduction

This paper applies machine learning for accelerating computations of market-
clearing outcomes in short-term electricity markets, such as hour-ahead or
real-time dispatches and locational marginal prices (LMPs), with the intention
of enabling high-speed and high-accuracy computations for a large number of
net load (load minus renewables) conditions. In particular, we are motivated
by the physics-informed machine learning (PIML) schemes outlined in [1]-[2]
for parameter and state estimation of large-scale power transmission networks,
and seek to extend these approaches to enable accurate and coherent prediction
of power dispatches and LMPs. The need for coherence arises from the
requirement for dispatch and LMP predictions to satisfy market design
properties (e.g. efficiency, cost recovery, and revenue adequacy, [3]), which
is a particular challenge for emerging stochastic market designs [4]. To this
end, we seek to relate primal (dispatch) and dual (LMPs) Optimal Power Flow
(OPF) solutions by internalizing conditions for market efficiency, cost
recovery, and revenue adequacy in the proposed machine learning approach.
Therefore, we develop a variant of the “classification-then-regression” (CTR)
method proposed in [5]. CTR is used to predict generator dispatch, governed by
SCED, through regression, after employing classification to discern which
generators are “saturated,” i.e., functioning at their minimum or maximum
capacities. We intend to follow the CTR framework for power lines, determining
which lines are saturated at their thermal limits. These thermal limits, in
addition to other constraints, can be encapsulated via the Karush-Kuhn-Tucker
(KKT) conditions (see, e.g., [6]), which are necessary for optimality of the
solution (and sufficient under conditions satisfied under DC-OPF assumptions)
and relate both primal and dual variables. By the KKT condition of
complementary slackness, lines which are not saturated should have a
corresponding dual variable of zero. In addition, as will be derived below,
and is discussed in [7], LMPs can be computed directly via a system of
equations relying on the KKT conditions. As such, CTR will enable us to
quickly solve for the variables needed to predict LMPs, representing an
improvement in how LMPs are computed via ML by relating them to primal
(dispatch) solutions using the market efficiency, cost recovery, and revenue
adequacy requirements. The method outlined will allow system operators and
market participants to quickly ascertain risk, and respond with strategic,
cost and profit-effective operational decisions and bids.

In addition to the PIML schemes for power systems in [1]-[2], ML techniques
have become increasingly prominent in enhancing the solution of optimal power
flow problems. The impetus for injecting ML into OPF is almost universally to
improve computational efficiency. To this end, a wide variety of paradigms
have been explored. Some examples include the use of deep neural networks
(NNs) [8]; graphical and convolutional NNs [9], [1]; random forests to resolve
voltage magnitudes and phases [10]; and deep reinforcement learning in systems
with high renewable penetration [11]. To further expedite the process, the
stacked extreme learning machine (SELM) framework has been explored to bypass
the need to learn parameters of a NN [12]. Most akin to our present work are
efforts to accelerate the solution of OPF by identifying sets of active
constraints. The general philosophy behind this method, and an illustration
that the number of “relevant” active sets at the optimal solution of an OPF is
few, is explored in detail in [13]. A more focused explanation and
illustration on DC-OPF is given in [14], and the notion that OPF can be
reduced to solving a system of linear equations once the active sets are
defined is introduced. In particular, a NN-based classifier is developed for
determining the most likely “active sets” corresponding to the DC-OPF
solution. Given a set of manifestations of uncertainty from, e.g., renewables
and demand patterns, the classifier returns a set $\mathcal{A}$ of active
sets, which should contain the “optimal active set” with high probability.
This is, by contrast, more efficient than traditional “ensemble methods,” in
which _all_ possible active sets are considered. This general methodology is
also followed in [5], on identifying sets of active generators. Additionally,
a variant of the active set philosophy was successfully applied by the winners
of the ARPA-E GO competition in 2018–2019 [15].

With the exception of [5], most of the above and related works are posed from
the perspective of a system operator aiming to determine and deploy optimal
dispatch. In addition, they are concerned mostly with the primal OPF problem.
Furthermore, in those works which do explicitly address ML for pricing of
generators, such as [16], a map directly to LMPs (or, to be precise, “LMP
spread”) is proposed, rather than solving for LMPs algebraically after
intermediate quantities are learned. In our work, we consider the dual problem
to DC-OPF, from which locational marginal prices (LMPs) of generators may be
computed. This tactic yields an ML scheme that results in LMPs which are
faithful to the constraints imposed during OPF, and is therefore of use to
energy market participants as well as system operators.

In the remainder of the manuscript, we state two formulations of DC-optimal
power flow (DC-OPF), using the latter to construct a system of equations from
which all primal and dual variables may be solved. We then outline how LMPs
may be computed as a function of these primal and dual variables. The heart of
the mansucript is a deep learning scheme that expedites the solution of OPF
and computation of LMPs, by classifying a priori which lines in the grid are
saturated and incorporating the corresponding thermal constraints, satisfied
with equality, into the aforementioned system of equations. The paper
concludes with numerical evaluation of this classification, and of the scheme
in general.

##  II Problem formulation

An optimal power flow (OPF) is the task of finding the economic dispatch,
i.e., the least-cost way to generate enough electricity to meet demand. In the
following, we first provide a synopsis of the traditional formulation of a DC-
OPF and how it can be rewritten in another equivalent form based on the power
transfer distribution factors (PTDFs) 111PTDF generalizes to a more general
case of linearized AC-OPF, which we plan to discuss in the extended version of
the manuscript.. In the DC approximation [17], the ohmic losses in the system
are neglected, and the power lines are solely characterized by their
susceptances $b_{ij}$. In the present work, we focus on the applications of
the DC-OPF to real-time dispatches. It is assumed that unit commitment has
been determined a priori, e.g., via security-constrained unit commitment
(SCUC) [18]. For the sake of simplicity, we assume that the generators have
quadratic costs, with quadratic coefficients $q_{g}$ and linear coefficients
$c_{g}$. The traditional formulation of the DC-OPF problem is

$\displaystyle\min_{p_{g}}$
$\displaystyle\sum_{g}\big{(}q_{g}p_{g}^{2}+c_{g}p_{g}\big{)}\,,$ (1)
$\displaystyle{\rm s.t.}$ $\displaystyle\forall g:p_{g}^{\min}\leq p_{g}\leq
p_{g}^{\max}\,,$ (2) $\displaystyle\forall i:\sum_{j\sim
i}b_{ij}(\theta_{i}-\theta_{j})=\sum_{g\in G_{i}}p_{g}-\ell_{i}\,,$ (3)
$\displaystyle\forall j\sim i:-f_{ij}^{\max}\leq
b_{ij}(\theta_{i}-\theta_{j})\leq f_{ij}^{\max}\,,$ (4)
$\displaystyle\theta_{\rm slack}=0\,.$ (5)

where $G_{i}$ is the set of generators connected to bus $i$ and $j\sim i$
denotes the neighbours of bus $i$. The objective is to minimize the total cost
of generation under some operational constraints. Constraint (2) enforces the
power flow equations. Constraint (3) ensures that generators do not surpass
their capacities when dispatched. Constraint (4) enforces that the power flows
must be smaller in absolute value than their thermal limits $f_{ij}^{\max}$.
Finally, constraint (5) establishes a reference (or slack bus) to which all
voltage phases are compared.

The PTDF matrix $\bm{\Phi}$ allows us to directly obtain the power flows in
the system from load $\ell_{i}$ and generator output $p_{g}$. For the sake of
completeness, we briefly show how the PTDF matrix is obtained. Let
$\bm{B}=\bm{A}\,{\rm diag}(\bm{b})\,\bm{A}^{\top}$, with $\bm{b}$ being the
vector composed of the line susceptances and $\bm{A}$ the incidence matrix,
defined as

$A_{ik}=\begin{cases}-1,&\text{ if line $k$ starts at bus $i$}\,,\\\
\phantom{-}1,&\text{ if line $k$ ends at bus $i$}\,,\\\ \phantom{-}0,&\text{
otherwise}\,.\end{cases}$ (6)

With these new definitions, Eqs. (3) and (4) can be rewritten in a matrix form
as

$\displaystyle\bm{P}$ $\displaystyle=\bm{B}\bm{\theta}\,,$ (7)
$\displaystyle\bm{F}$ $\displaystyle={\rm
diag}(\bm{b})\bm{A}^{\top}\bm{\theta}\,,$ (8)

where $P_{i}=\sum_{g\in G_{i}}p_{g}-l_{i}$ and with the lines being relabeled
as $ij\rightarrow k$. $\bm{B}$ has one zero eigenvalue associated with the
eigenvector $1_{b}$, consisting of a vector of $1$’s for each bus, and hence
solving Eq. (7) for $\bm{\theta}$ requires the use of a pseudoinverse
$\bm{B}^{\dagger}$. Finally, injecting the result in Eq. 8, we obtain

$\bm{F}={\rm
diag}(\bm{b})\bm{A}^{\top}\bm{B}^{\dagger}\bm{P}\equiv\Phi\bm{P}\,.$ (9)

The pseudoinverse matrix $\bm{B}^{\dagger}$ has also one zero eigenvalue
associated with the eigenvector $1_{b}$, and therefore $\sum_{i}\Phi_{ki}=0$.

Using $\Phi$, it is possible to solve the OPF problem without explicitly
computing $\theta_{j}$. Hence, we reformulate the OPF problem (1)-(5) as

$\displaystyle\min_{p_{g}}$
$\displaystyle\sum_{g}\big{(}q_{g}p_{g}^{2}+c_{g}p_{g}\big{)}\,,$ (10)
$\displaystyle{\rm s.t.}$ $\displaystyle\forall g:\,p_{g}^{\min}\leq p_{g}\leq
p_{g}^{\max}\,,$ (11) $\displaystyle\sum_{g}p_{g}=\sum_{i}\ell_{i}\,,$ (12)
$\displaystyle\forall k:-f_{k}^{\max}\leq\sum_{i}\Phi_{ki}\Big{(}\sum_{g\in
G_{i}}p_{g}-\ell_{i}\Big{)}\leq f_{k}^{\max}.$ (13)

This reformulation allows to bypass the introduction of dual variables
associated with voltage phases, which yields computational benefits to be
detailed in the following sections.

In what follows, we make the following
Assumption: As the unit commitment is assumed to be provided and fixed, we
only consider non-committed active (“free” for simplicity) generators, i.e.
$p_{g}^{\min}<p_{g}<p_{g}^{\max}$, and hence no Lagrangian multipliers were
assigned to enforce these min/max constraints. Furthermore, the dispatch
associated with committed generators is incorporated as “negative load” for
the buses to which they are associated.

From the DC-OPF problem (10) - (13) and the above assumption, we form a
partial Lagrangian function.

$\displaystyle\mathcal{L}(p_{g},$
$\displaystyle\lambda,\mu_{k},\nu_{k})\vcentcolon=\sum_{g}\big{(}q_{g}p_{g}^{2}+c_{g}p_{g}\big{)}+\lambda\Big{(}\sum_{g}p_{g}-\sum_{i}\ell_{i}\Big{)}$
$\displaystyle+\sum_{k}\mu_{k}\Big{[}\sum_{i}\Phi_{ki}\Big{(}\sum_{g\in
G_{i}}p_{g}-\ell_{i}\Big{)}-f_{k}^{\text{max}}\Big{]}$
$\displaystyle-\sum_{k}\nu_{k}\Big{[}\sum_{i}\Phi_{ki}\Big{(}\sum_{g\in
G_{i}}p_{g}-\ell_{i}\Big{)}+f_{k}^{\text{max}}\Big{]}\,.$ (14)

Once the OPF problem is solved, the locational marginal price (LMP) at every
bus is obtained as:

$\forall i:{\rm
LMP}_{i}\equiv\frac{\partial\mathcal{L}}{\partial\ell_{i}}=-\lambda+\sum_{k}\Phi_{ki}\big{(}\nu_{k}-\mu_{k}\big{)}\,,$
(15)

where the first term represents the energy components of LMPs and the second
term represents the congestion component of LMPs. If there is no congestion in
the transmission network, then the second term disappears from (15), as there
are no binding power flow limits and their respective dual variables are zero.

Finally, by construction the solution is _coherent_ , i.e. it provides
rigorous guarantees for the following key (in)equalities:
$\bullet$ Revenue Adequacy: The sum of payments collected by the market from
electricity buyers (i.e., loads) must be equal or exceed the sum of payments
made by the market to the electricity sellers (i.e., generators):

$\sum_{g}{\rm LMP}_{i(g)}p_{g}\leq\sum_{i}{\rm LMP}_{i}\ell_{i}\,,$ (16)

where $i(g)$ is the index of the bus to which the generator $g$ is connected.
$\bullet$ Cost Recovery: Each generator does not produce if their marginal
cost is higher than their respective LMP; that is, every generator must break
even from the participation in the market:

$\forall g\colon\,q_{p}p_{g}^{2}+c_{g}p_{g}\leq{\rm LMP}_{i(g)}p_{g}\,,$ (17)

$\bullet$ Efficiency via Strong Duality: Since DC-OPF, as stated, satisfies
Slater’s conditions, the optimal values of the primal and dual objectives are
equal [6], which indicates the least-cost primal objective value and thus
efficiency. As such,

$\displaystyle\lambda\Big{(}\sum_{g}p_{g}-\\!\sum_{i}\ell_{i}\Big{)}+\sum_{k}\mu_{k}\Big{[}\\!\sum_{i}\Phi_{ki}\Big{(}\sum_{g\in
G_{i}}p_{g}-\ell_{i}\Big{)}-f_{k}^{\text{max}}\Big{]}$
$\displaystyle-\sum_{k}\nu_{k}\Big{[}\sum_{i}\Phi_{ki}\Big{(}\sum_{g\in
G_{i}}p_{g}-\ell_{i}\Big{)}+f_{k}^{\text{max}}\Big{]}=0\,.$ (18)

##  III Fast ML approach to solving the OPF problem

In this section, we present our approach to solving the OPF problem. It
provides a fast way of obtaining generation dispatches and LMPs. We call the
approach identification-then-solving (ITS). We will detail ITS in the
following paragraphs, but in short, the procedure amounts to:

  1. 1.

finding a mapping from loads $\ell_{i}$ to binding line constraints,

  2. 2.

solving a system of linear equations derived from the OPF constraints, and

  3. 3.

recovering the LMPs.

###  III-A Identification of Binding Line Constraints

_Binding constraints_ are those that are active at the optimal DC-OPF
solution. Our approach originates from the empirical finding that there is
typically only a small set of lines which are loaded to their maximal
capacities $f_{k}^{\max}$ [14], and that this set does not vary much from one
solved OPF instance to another, given the same system parameters and similar
load configurations. This is to be expected, as it is well-known that power
systems usually have transmission “bottle-necks”, i.e. portions that tend to
be congested (e.g., as the interface between New York City and Long Island and
the rest of the NYISO transmission network [19]). Therefore, most of the
constraints defined in Eq. (13) are never active. From the complementary
slackness condition [6], if a line constraint is not active, then its dual
variable, $\mu_{k}$ or $\nu_{k}$, is zero.

It should therefore be possible to correctly identify these lines and,
moreover, to exploit their properties to gain information about the state of
the system. In other words, we seek to find a map from an input to the binding
line constraints

$\bm{\ell}\longrightarrow\bigg{[}\\!\\!\begin{array}[]{c}\mathcal{B}(\bm{\mu})\\\
\mathcal{B}(\bm{\nu})\end{array}\\!\\!\bigg{]}\,,$ (19)

where $\mathcal{B}(x)=\\{1:\text{ if }x\neq 0\,,\;0:\text{ otherwise}\\}$. In
the present work, we limit ourselves to inputs consisting only of the loads
$\ell_{i}$, but one may imagine more advanced inputs. Multiple approaches can
be undertaken to find this map. Here, we opt for using a simple deep NN for
simplicity of its implementation.

###  III-B Solving the OPF as a Set of Linear Equations

We now show that once the identification stage is concluded as described
above, the OPF optimization problem is equivalent to solving a system of
linear equations (SLE). As previously mentioned in our assumption, from unit
commitment, we assume that we know which generators are producing at their
full capacity, and we treat them as negative loads. In the following, $p_{g}$
denotes power outputs of the “free” generators, and $\ell_{i}$ the residual
load obtained after the subtraction of committed generators.

From $\mathcal{B}(\bm{\mu})$ and $\mathcal{B}(\bm{\nu})$, we define
$\mathcal{C}_{\nu}=\\{k\,\colon\,\nu_{k}\neq 0\\}$ and
$\mathcal{C}_{\mu}=\\{k,\colon\,\mu_{k},\neq 0\\}$ which gather the indices of
non-zero dual variables associated with binding line constraints to be known.
(Note that there is a one-to-one correspondence between these sets and
$\mathcal{B}(\bm{\mu})$ and $\mathcal{B}(\bm{\nu})$.) Assuming these sets to
be known, the optimization problem defined by Eqs. (10) - (13) becomes
equivalent to solving the following set of linear equations

$\displaystyle\forall\mu_{k}^{*}:\,$
$\displaystyle\sum_{i}\Phi_{ki}\Big{(}\sum_{g\in
G_{i}}p_{g}-\ell_{i}\Big{)}=f_{k}^{\text{max}}\,,$ (20)
$\displaystyle\forall\nu_{k}^{*}:\,$
$\displaystyle\sum_{i}\Phi_{ki}\Big{(}\sum_{g\in
G_{i}}p_{g}-\ell_{i}\Big{)}=-f_{k}^{\text{max}}\,,$ (21) $\displaystyle\forall
g:\,$ $\displaystyle\,2q_{g}p_{g}+c_{g}+\lambda+\\!\\!\sum_{k\in
C_{\mu}}\\!\Phi_{ki}\mu_{k}^{*}-\\!\\!\sum_{k\in
C_{\nu}}\\!\Phi_{ki}\nu_{k}^{*}=0\,,$ (22)
$\displaystyle\sum_{g}p_{g}=\sum_{i}\ell_{i}.$ (23)

Once these equations are solved, LMPs are straightforwardly obtained from
$\lambda$, $\mu_{k}^{*}$, and $\nu_{k}^{*}$ by applying Eq. (15). Importantly,
_this solution process can be performed significantly faster than solving the
optimization problem_. In particular, solving the system of linear equations
is already expedient, and benefits from efficient linear algebra utilities in
most programming environments. As such, the “rate-limiting step” is the
classification of binding line constraints, which can be made similarly
efficient if a suitable ML methodology is chosen. In our case, to be described
below, once the classification NN is trained and performing satisfactorily,
the process can immediately be repeated on an ensemble of samples.
Furthermore, ITS is clear and intuitive, accessible to practitioners with
fundamental knowledge of linear algebra and machine learning.

##  IV Numerical Experiments

In this section, we apply the ITS approach to a standard test case. We detail
how we build our data sets and identifier NN. We use our NN identifier to
learn the binding line constraints for different unit commitment
configurations. This information is then used to obtain the LMPs and
generation dispatch outputs. Most notably, we investigate how the results of
our ITS method deteriorate when it fails to correctly identify the binding
line constraints, highlighting the need for a judicious choice of ML
methodology for this crucial step.

###  IV-A Generating Data Sets

We demonstrate our method by applying it to the standard IEEE-118 test case.
Unit commitment configurations are obtained from different load configurations
corresponding to annual hourly profiles. We extract the 15 configurations that
occur the most often. For each of the configurations, a set of 500 samples is
obtained by adding noise to each of the individual loads. We use different
volatility strengths, with standard deviations of 1%, 5%, and 10% of the
nominal loads. Out of 15 configurations, we select three of them to present
our results. They are chosen to represent various degrees of success: for
Config. #1, our method is accurate in almost every case and performs well for
most samples, but starts to fail when volatility is high in Config. #2; and it
fails in several cases even with a relatively low volatility and a large
training set for Config. #3.

| Free generators
---|---
Config. #1 | $\\{3,5,11,12,18,30,34,40,42,43\\}$
Config. #2 | $\\{2,5,12,26,30,39\\}$
Config. #3 | $\\{5,12,14,20,30,37,39\\}$

###  IV-B Implementing the Identification Process

As we previously mentioned, we decided to use a neural network
$\mathbf{NN}_{\bm{\psi}}$ to identify the binding line constraints

$\left[\\!\\!\begin{array}[]{c}\bm{y}^{\nu}\\\ \bm{y}^{\mu}\\\
\end{array}\\!\\!\right]=\mathbf{NN}_{\bm{\psi}}(\bm{\ell})\,,$ (24)

It is trained to minimize the following loss function

$L_{\psi}=\\!\\!\sum_{s}\sum_{k}\\!\Big{[}L_{\rm
reg}\big{(}y_{k}^{\nu(s)}\\!\\!\\!,\mathcal{B}\big{(}\nu_{k}^{(s)}\big{)}+L_{\rm
reg}\big{(}y_{k}^{\nu(s)}\\!\\!\\!,\mathcal{B}\big{(}\nu_{k}^{(s)}\big{)}\big{)}\\!\Big{]},$
(25)

where $L_{\rm reg}$ is the cost function for a logistic regression, which is
defined as

$L_{\rm reg}(x,y)=\begin{cases}-\log(x),&\text{ if }y=1\,,\\\
-\log(1-x),&\text{ if }y=0\,.\end{cases}$ (26)

We use a 4-layer NN with 500 units per layer. The last layer uses a sigmoid
activation function, which means that the outputs of the neural network lay in
$]0,1[$. Once the NN is trained, we use the threshold of $0.5$ to transform
the outputs $y_{k}^{\nu(s)}$ and $y_{k}^{\mu(s)}$ to binary variables; i.e.,
if $y_{k}^{\nu(s)}>0.5$, set $y_{k}^{\nu(s)}$ to $1$, otherwise set to $0$,
and similar for $y_{k}^{\mu(s)}$. We use Flux.jl [20] to implement the NN and
train it.

###  IV-C Accuracy of the Identification

Table I presents the accuracy of our identification process. As explained
previously, Config. #1 was selected because our method performs well for it.
Indeed, there are only a few misidentifications, always representing less than
1% of the binding line constraints, which leads to excellent predictions for
LMPs and generation dispatch outputs. Config. #2 and #3 have far more
misidentifications. Surprisingly, larger training sets tend to increase their
number. As we will show in the next subsection, the number of
misidentifications is directly related to the mismatch, and binding line
identifications based on larger training sets tend to give better market-
clearing results. This also means that failing to identify certain binding
constraints is more consequential for the ultimate dispatch and LMP
predictions than others.

|  | Config. #1 | Config. #2 | Config. #3
---|---|---|---|---
$\sigma$ | size | training | testing | training | testing | training | testing
1% | 50 | 0(0.00) | 0(0.00) | 0(0.00) | 34(0.02) | 1(0.01) | 66(0.05)
1% | 100 | 0(0.00) | 0(0.00) | 0(0.00) | 48(0.03) | 0(0.00) | 273(0.20)
1% | 200 | 0(0.00) | 6(0.00) | 0(0.00) | 104(0.06) | 0(0.00) | 342(0.25)
5% | 50 | 0(0.00) | 0(0.00) | 0(0.00) | 16(0.01) | 9(0.03) | 50(0.05)
5% | 100 | 0(0.00) | 0(0.00) | 0(0.00) | 43(0.03) | 0(0.00) | 165(0.14)
5% | 200 | 0(0.00) | 6(0.00) | 0(0.00) | 79(0.05) | 0(0.00) | 232(0.19)
10% | 50 | 0(0.00) | 0(0.00) | 4(0.00) | 19(0.02) | 26(0.05) | 36(0.04)
10% | 100 | 0(0.00) | 0(0.00) | 0(0.00) | 24(0.02) | 6(0.01) | 103(0.01)
10% | 200 | 0(0.00) | 2(0.00) | 0(0.00) | 35(0.03) | 1(0.00) | 133(0.15)
Table I: Number of misidentifications, i.e. $y_{k}^{\nu|\mu(s)}=1$ and
$\mathcal{B}\big{(}\nu|\mu_{k}^{(s)}\big{)}=0$ or vice versa, over the
training and testing sets. The ratio of misidentifications to number of
binding line constraints is given in parenthesis.

We have thus far discussed results for three representative configurations.
Table II shows the fraction of misidentifications over the testing set for the
remaining 12 unit commitment configurations under inspection. On average,
their difficulty in terms of training is somewhere between that of Config. #1
and Config. #2. Namely, the fraction of misidentifications over the testing
set is generally less than 10%, suggesting that our method is tractable and
yields accurate classification results in the vast majority of typical cases.
If encountering a limited number of difficult unit commitment configurations,
one is advised to increase the size of the training set or perform manual
computations after training to ensure the fidelity of the classifications
predicted. From a system operations standpoint, the correct course of action
will depend on the particular situation, and indeed, it may be worthwhile to
consider alternative methods for the computation of generation dispatches and
LMPs for such cases.

$\sigma$ | #4 | #5 | #6 | #7 | #8 | #9 | #10 | #11 | #12 | #13 | #14 | #15
---|---|---|---|---|---|---|---|---|---|---|---|---
1% | 0.00 | 0.00 | 0.00 | 0.00 | 0.00 | 0.02 | 0.01 | 0.00 | 0.00 | 0.01 | 0.01 | 0.00
5% | 0.04 | 0.01 | 0.02 | 0.00 | 0.03 | 0.07 | 0.02 | 0.00 | 0.02 | 0.06 | 0.05 | 0.01
10% | 0.05 | 0.03 | 0.04 | 0.02 | 0.06 | 0.08 | 0.06 | 0.02 | 0.06 | 0.09 | 0.05 | 0.05
Table II: Fraction of misidentifications over testing set for the 13 unit
commitment configurations not presented in Table I.

###  IV-D Comparison of the Solutions of the SLE and of the OPF Problem

Figure 1: Mismatch between obtained LMPs (top) and generator outputs (bottom)
for different load volatility standard deviations: 1% (blue), 5% (orange) and
10% (green) of their nominal values, and the training set consists of: 50
(solid), 100 (dashed) are 200 (dotted) samples. Each panel shows the fraction
of the testing set that get maximal mismatches smaller than the ordinate.

To evaluate the performance of ITS below, we consider the maximum absolute
error (MAE) in LMPs and dispatch across each bus. This choice of metric is
natural, given that market participants are often interested in the most
extreme violations of what is expected in normal operations. Such deviations
lead to “price spikes” in the LMPs, which trading agents can profit off of if
predicted accurately. However, if these spikes are incorrectly anticipated,
perhaps as artifacts of the ITS method, it is not as straightforward for
trading agents to identify their maximally profitable strategy.

Fig. 1 compares the results of our ITS methods with the ones obtained by
solving directly the optimization problem (10)-(13). For Config. #1, the
results agree, apart from a very small number of samples when the load
volatility is high. For Config. #2, the higher the load volatility, the higher
the mismatches. Training over a large set seems to provide better results. We
still obtain excellent results in more than 85% of cases with a training of
200 samples. Finally, for Config. #3, our method is off for a significant
fraction of the testing set. We hypothesize that misidentifying a binding line
with a large associated dual variable affects the accuracy of our method more
than the total number of misidentifications. If this hypothesis turns out
correct, one can take this into account by using $\mu_{k}^{*}$ and
$\nu_{k}^{*}$ as weights in the loss function. This effectively introduces a
penalty for each individual misidentification, providing a safeguard against
the sensitivity of dispatch and LMP computations to the success of the
identification step, or lack thereof.

We analyze if the results we obtain with the proposed method satisfy Eqs. (16)
- (18). Table III shows that Revenue Adequacy and Strong Duality, as stated in
Eq. (16) and (18), respectively, are almost always satisfied. Cost Recovery as
stated in Eq. (17) is more prone to be violated. When the load volatility is
low, Cost Recovery is satisfied in most cases. As the volatility increases,
the results start to violate Cost Recovery in a majority of cases. The
magnitude of such violations depend on a chosen unit commitment configuration
and is consistent with the current industry practice, that is; such violations
can be settled in a out-of-market manner using an uplift mechanism [21].

$\sigma$ | Revenue Adequacy | Strong Duality | Cost Recovery
---|---|---|---
1% | 1.000 | 1.000 | 0.808
5% | 0.999 | 0.997 | 0.352
10% | 0.999 | 0.992 | 0.060
Table III: Fraction of the testing samples (over all 15 unit commitment
configurations) that satisfy the key (in)equalities as stated in Eq.
(16)-(18).

##  V Conclusion

In this manuscript, we suggest an identification-then-solving machine learning
scheme, consisting of a neural network augmented by traditional linear
algebra, to substitute for solution of LMPs. The engine of ITS is to exploit
the conditions resulting from balance of load and dispatch, complementary
slackness, and stationarity at the point of optimal dispatch. Our scheme
relies heavily, almost exclusively, on using a neural network to classify
which lines in the transmission grid are binding, and it is observed that the
accuracy of both the primal and dual OPF solutions are highly sensitive to the
accuracy of the classification step. While the classifier itself performs
reasonably well, in order for our proposed scheme to have value, nearly all
binding lines must be correctly anticipated, so that all significantly large
dual variables are incorporated into the constructed system. Nonetheless, it
is encouraging that ITS almost exactly replicates both generator dispatch and
LMPs when the classification scheme achieves its goals. In subsequent works,
we intend to improve the fidelity of the neural network, namely by taking
advantage of the topology of the grid at hand to construct a graphical neural
network. Once we achieve the desired level of performance with a GNN, we will
also conduct numerical experiments on a larger system, to ascertain the extent
to which ITS is robust against various contingencies in power grid operations.
Overall, ITS presents a theoretically and mechanically simple alternative to
computing dispatch and LMPs through OPF, and such methodology will be valuable
to rapid decision-making on the part of both system operators and energy
market participants.

## Acknowledgment

The authors would like to collectively acknowledge with gratitude the ARPA-E
PERFORM grant from which this work arose. In addition, R. Ferrando was
partially supported by the University of Arizona NSF RTG in Data-Driven
Discovery. All authors also thank Daniel Bienstock and Robert Mieth for useful
discussions, guidance, and technical assistance at various stages of this
work.

## Appendix: Notations

Throughout our derivations, matrix $\bm{D}$ distributes the generators to the
buses they are connected to, i.e.

$D_{ig}=\left\\{\begin{array}[]{l}1\,,\text{ if generator $g$ is connected to
bus $i$}\,,\\\ 0\,,\text{ otherwise}\,.\end{array}\right.$ (27)

The operative system of equations (20)-(23) may be more compactly represented
in matrix form. Below, $\bm{Q}={\rm diag}(\\{q_{g}\\})$ is a matrix of
quadratic costs, and $\bm{\Phi}_{\nu}$, $\bm{\Phi}_{\mu}$ are reduced PTDF
matrices, only accounting for those unsaturated lines in which
$\mu_{k},\nu_{k}\neq 0$.

$\\!\\!\\!\begin{bmatrix}\\!\\!2\bm{Q}&\\!\\!\\!\\!\\!\\!\\!\\!\\!(\bm{\Phi}_{\nu}\bm{D})^{\top}&\\!\\!\\!\\!\\!(\bm{\Phi}_{\mu}\bm{D})^{\top}&\\!\\!\\!\\!\bm{1}_{\rm
g}\\\ \bm{\Phi}_{\nu}\bm{D}&\bm{0}&\bm{0}&0\\\
-\bm{\Phi}_{\mu}\bm{D}&\bm{0}&\bm{0}&0\\\ \bm{1}_{\rm
g}^{\top}&\bm{0}&\bm{0}&0\\\ \end{bmatrix}\cdot\begin{bmatrix}\bm{p}\\\
\bm{\nu}^{*}\\\ \bm{\mu}^{*}\\\
\lambda\end{bmatrix}\\!=\\!\begin{bmatrix}-\bm{c}\\\
\bm{f}_{\nu}-\bm{\Phi}_{\nu}\bm{l}\\\ \bm{f}_{\mu}+\bm{\Phi}_{\mu}\bm{l}\\\
\sum_{i}\ell_{i}\\\ \end{bmatrix}.$ (28)

The above matrix may be nearly degenerate if; e.g., two generators share the
same “coordinates” such as costs. As such, when solving the resulting system
of equations, extra care must be taken when solving numerically; for example,
one may left-multiply by pseudo-inverse of the coefficient matrix of the
variables. Finally, it is worth noting that following our chosen DC-OPF
formulation, we do not have to introduce an additional $N_{\rm b}+1$ equations
and unknowns associated with voltage phase angles.

## References

  * [1] L. Pagnier and M. Chertkov, “Embedding power flow into machine learning for parameter and state estimation,” _arXiv:2103.14251_ , 2021.
  * [2] ——, “Physics-informed graphical neural network for parameter & state estimations in power systems,” _arXiv:2102.06349_ , 2021.
  * [3] H. R. Varian, _Microeconomic analysis_ , 3rd ed. New York: Norton, 1992.
  * [4] Y. Dvorkin, “A chance-constrained stochastic electricity market,” _IEEE Transactions on Power Systems_ , vol. 35, no. 4, pp. 2993–3003, 2020.
  * [5] W. Chen, S. Park, M. Tanneau, and P. V. Hentenryck, “Learning optimization proxies for large-scale security-constrained economic dispatch,” _arXiv:2112.13469_ , 2021.
  * [6] A. Beck, _Introduction to nonlinear optimization: Theory, algorithms, and applications with MATLAB_. SIAM, 2014.
  * [7] E. Litvinov, T. Zheng, G. Rosenwald, and P. Shamsollahi, “Marginal loss modeling in lmp calculation,” _IEEE Transactions on Power Systems_ , vol. 19, no. 2, pp. 880–888, 2004.
  * [8] X. Pan, T. Zhao, and M. Chen, “Deepopf: Deep neural network for dc optimal power flow,” 2020.
  * [9] T. Falconer and L. Mones, “Leveraging power grid topology in machine learning assisted optimal power flow,” _CoRR_ , vol. abs/2110.00306, 2021. [Online]. Available: https://arxiv.org/abs/2110.00306
  * [10] J. Rahman, C. Feng, and J. Zhang, “Machine learning-aided security constrained optimal power flow,” in _2020 IEEE Power Energy Society General Meeting (PESGM)_ , 2020, pp. 1–5.
  * [11] J. Li, R. Zhang, H. Wang, Z. Liu, H. Lai, and Y. Zhang, “Deep reinforcement learning for optimal power flow with renewables using spatial-temporal graph information,” 2021.
  * [12] X. Lei, Z. Yang, J. Yu, J. Zhao, Q. Gao, and H. Yu, “Data-driven optimal power flow: A physics-informed machine learning approach,” 2020.
  * [13] S. Misra, L. Roald, and Y. Ng, “Learning for constrained optimization: Identifying optimal active constraint sets,” 2019.
  * [14] D. Deka and S. Misra, “Learning for DC-OPF: Classifying active sets using neural nets,” in _2019 IEEE Milan PowerTech, arXiv/1902.05607_ , pp. 1–6.
  * [15] F. E. Curtis, D. K. Molzahn, S. Tu, A. Wächter, E. Wei, and E. Wong, “A decomposition algorithm for large-scale security-constrained ac optimal power flow,” _arXiv:2110.01737_ , 2021.
  * [16] Y. Li, N. Yu, and W. Wang, “Machine learning-driven virtual bidding with electricity market efficiency analysis,” 2021.
  * [17] W. Stevenson Jr and J. Grainger, _Power System Analysis_. McGraw-Hill, 1994.
  * [18] R. Fernández-Blanco, Y. Dvorkin, and M. A. Ortega-Vazquez, “Probabilistic security-constrained unit commitment with generation and transmission contingencies,” _IEEE Transactions on Power Systems_ , vol. 32, no. 1, pp. 228–239, 2017.
  * [19] N. S. Services, “Guide 01: Market participants user’s guide,” 2021. [Online]. Available: https://www.nyiso.com/documents/20142/3625950/mpug.pdf/
  * [20] M. Innes, “Flux: Elegant machine learning with julia,” _Journal of Open Source Software_ , 2018.
  * [21] Y. M. Al-Abdullah, M. Abdi-Khorsand, and K. W. Hedman, “The role of out-of-market corrections in day-ahead scheduling,” _IEEE Transactions on Power Systems_ , vol. 30, no. 4, pp. 1937–1946, 2015.