Learning to Solve Large-Scale Security-Constrained Unit Commitment Problems

TL;DR

This paper introduces machine learning techniques to improve the efficiency of solving large-scale security-constrained unit commitment problems in power systems, achieving significant speedups while maintaining solution quality.

Contribution

The paper presents novel ML methods to predict redundant constraints and feasible regions, significantly accelerating MIP-based SCUC solutions with guarantees.

Findings

Average 4.3x faster solutions with optimality guarantees

Average 10.2x faster solutions without guarantees

Method shows robustness against dataset shift

Abstract

Security-Constrained Unit Commitment (SCUC) is a fundamental problem in power systems and electricity markets. In practical settings, SCUC is repeatedly solved via Mixed-Integer Linear Programming, sometimes multiple times per day, with only minor changes in input data. In this work, we propose a number of machine learning (ML) techniques to effectively extract information from previously solved instances in order to significantly improve the computational performance of MIP solvers when solving similar instances in the future. Based on statistical data, we predict redundant constraints in the formulation, good initial feasible solutions and affine subspaces where the optimal solution is likely to lie, leading to significant reduction in problem size. Computational results on a diverse set of realistic and large-scale instances show that, using the proposed techniques, SCUC can be solved on average 4.3x faster with optimality guarantees, and 10.2x faster without optimality guarantees, but with no observed reduction in solution quality. Out-of-distribution experiments provides evidence that the method is somewhat robust against dataset shift.

Tables

Table 1. Table 1: Hyper-parameters for affine subspace predictor.

Hyperplane	$z^{\text{FIX}}$	$z^{\text{MIN}}$	$z^{\text{MAX}}$	${\alpha }^{\text{R}}$	${\alpha }^{\text{P}}$
$x_{gt}=0$	1.000	0.250	0.750	0.90	0.90
$x_{gt}=1$	1.000	0.250	0.750	0.75	0.75
$x_{gt}=x_{g,t+1}$	0.975	0.025	0.975	0.50	0.75

Table 2. Table 2: Size of selected instances.

Instance	Buses	Units	Lines
case1888rte	1,888	297	2,531
case1951rte	1,951	391	2,596
case2848rte	2,848	547	3,776
case3012wp	3,012	502	3,572
case3375wp	3,374	596	4,161
case6468rte	6,468	1,295	9,000
case6470rte	6,470	1,330	9,005
case6495rte	6,495	1,372	9,019
case6515rte	6,515	1,388	9,037

Table 3. Table 3: Impact of different transmission predictors on number of constraints added per time period, number of iterations and running time.

Method	Iterations	Violations	Time (s)	Speedup

Equations

$$\sum _{g\in G}{c}_g(x_{g\bullet },y_{g\bullet })$$

$$(x_{g\bullet },y_{g\bullet })\in {\mathcal{G}}_g$$

$$\sum _{g\in G}{y}_{gt}=\sum _{b\in B}{d}_{bt}$$

$$-F_l^0\le \sum _{b\in B}{\delta }_{lb}^0\left(\sum _{g\in G_b}{y}_{gt}-d_{bt}\right)\le F_l^0$$

$$-F_l^c\le \sum _{b\in B}{\delta }_{lb}^c\left(\sum _{g\in G_b}{y}_{gt}-d_{bt}\right)\le F_l^c$$

$$x_{gt}\in \{0,1\}$$

$$y_{gt}\ge 0$$

$$\underset{\pi ,{\pi }_0,\alpha }{\mathrm{minimize}}$$

$$h_i\left({\pi }^T{p}^i-{\pi }_0\right)\ge 1-{\alpha }_i$$

$${\alpha }_i\ge 0,$$

$$\sum _{i=1}^k\sum _{(p,h)\in {\mathcal{D}}_i}\frac{f\left({\varphi }_i(p),h\right)}{s},$$

$$f_{\text{precision}}(\tilde{h},h)=\frac{\left|\{j:{\tilde{h}}_j=1\wedge h_j=1\}\right|}{\left|\{j:{\tilde{h}}_j=1\}\right|},$$

$$f_{\text{recall}}(\tilde{h},h)=\frac{\left|\{j:{\tilde{h}}_j=1\wedge h_j=1\}\right|}{\left|\{j:h_j=1\}\right|}.$$

$$-F_l^c\le \sum _{b\in B}{\delta }_{lb}^c\left(\sum _{g\in G_b}{y}_{gt}-d_{bt}\right)\le F_l^c$$

$$h_{l,c,t}\in \{\text{ENFORCE},\text{RELAX}\}$$

$$x_{gt}$$

$$\max \{\bar{z},1-\bar{z}\}\cdot (1-{\alpha }^{\text{P}})+{\alpha }^{\text{P}},$$

Full text

1 Introduction

Security-Constrained Unit Commitment (SCUC) is one of the most fundamental optimization problems in power systems and electricity markets, being solved a number of times daily by major reliability coordinators, independent system operators (ISOs), and regional transmission organizations (RTOs) to schedule 3.8 trillion kWh of energy production and clear a $400 billion electricity market annually in the United States (EiaWebsite). The problem asks for the most cost-efficient schedule for power generation units under a large number of physical, operational and economic constraints. Among other requirements, a feasible generation schedule must not only guarantee that enough power is being produced during each hour to satisfy the demand, but also that no transmission lines are overloaded during the delivery of the electric power. Furthermore, the schedule must be robust against small variations in demands and against the unexpected loss of a small number of network elements, such as generation units and transmission lines.

The computational performance of SCUC is an extremely important practical issue, given the very short clearing windows enforced in most electricity markets. For instance, in Midcontinent Independent System Operator (or MISO, one of the largest electricity markets in North America), systems operators only have 3 or 4 hours, after receiving all bids and offers, to produce a low-cost and fair generation schedule. During this short time window, SCUC must be solved multiple times, under a number of different scenarios. Using state-of-art software, realistic large-scale instances can be usually solved to 0.1% optimality within approximately 20 minutes of running time (ChenCastoWangWangWangWan2016). Improvements in computational performance would allow markets to implement numerous enhancements that could bring significant economic benefits and improved market efficiency, such as more accurate modelling of enhanced combined cycle units, higher resolution for production cost curves, sub-hourly dispatch, longer planning horizons, among others. Alternatively, a simple reduction of the optimality gap to 0.05% without sacrificing computational running times could lead to significant cost savings, considering the $400 billion market size in U.S. alone.

In the past, SCUC has been solved using various optimization techniques including priority lists (BurnsGibson1975), dynamic programming (Lowery1966) and Lagrangian relaxation (MerlinSandrin1983). Nowadays, SCUC is most commonly formulated as a Mixed-Integer Linear Programming (MIP). The introduction of MIP solvers allowed system operators to easily model more realistic constraints that had been traditionally hard to enforce with previous methods, and have benefited from the constantly advancing capabilities of modern solvers. Since the first MIP formulations, proposed in the 1960s (Garver1962), a number of alternative formulations and polyhedral results have been described (AtakanLulliSen2018; OstrowskiAnjosVannelli2012; RajanTakriti2005; GentileMoralesEspanaRamos2017; DamciKurtKucukyavuzRajanAtamturk2016; espana2015tight; lee2004min; pan2016polyhedral; Knueven2018). There has also been active research on more practical techniques for improving the computational performance of the MIP approach, such as fast identification of redundant constraints (ZhaiGuanChengWu2010; ArdakaniBouffard2015; XavierFengWangThimmapuram2018), decomposition methods (FeizollahiCostleyAhmedGrijalva2015; Kim18TemporalDecomp), among others.

One aspect of SCUC that has been overlooked in previous research is that, in most practical settings, this problem is repeatedly solved with only minor changes in input data, often multiple times per day. Most of the system characteristics, such as the parameters describing the generation units or the topology of the transmission network, are kept almost entirely unchanged from solve to solve. While modern MIP solvers have the ability to reuse previous solutions as warm starts, it is well known that, for SCUC, providing the previous-day optimal solution brings only minor performance improvements (ChenCastoWangWangWangWan2016). In this work, we propose the use of machine learning to effectively extract information from previously solved instances, and show how to leverage this information to significantly accelerate the solution of similar instances in the future.

Machine learning (ML) is a discipline that studies how algorithms can automatically learn from large amounts of data and subsequently use the acquired knowledge to make decisions quickly and accurately. In the 1990s and early 2000s, researchers experimented using artificial neural networks (ANN) to solve simpler variants of SCUC on small-scale instances (SasakiWatanabeKubokawaYorinoYokoyama1992; WangShahidehpour1993; WalshOMalley1997; LiangKang2000). Even when considering very simplified versions of the problem, obtaining sufficiently high-quality solutions to SCUC via ANN proved very challenging, and the approach failed to replace existing methods. Similar explorations with evolutionary algorithms proved equally challenging (JusteKitaTanakaHasegawa1999).

In more recent years, there has been a growing interest, within the broader mathematical optimization community, in applying ML techniques to enhance the performance of current MIP solvers. ML has been used, for example, to automatically construct branching strategies (AlvarezLouveauxWehenkel2017; KhalilLeBodicSongNemhauserDilkina2016; LodiGiulia2017), to better parallelize branch-and-bound methods (AlvarezWehenkelLouveaux2014), to decide when to run primal heuristics (KhalilDilkinaNemhauserAhmedShao2017), to predict resolution outcome (FischettiLodiZarpellon2019), to decide solver parameters (BonamiLodiZarpellon2018) and to automatically construct decompositions (BassoCeselliTettamanzi2017). More broadly, there has been a renewed interest in using ML to tackle hard combinatorial problems. To give a few examples, deep and reinforcement learning have been used to produce heuristic solutions for optimization problems over graphs (KhalilDaiZhangDilkina2017; kool2018attention; MichelPierreAlexandreYossiriLouisMartin2018) and for stochastic integer programs (nair2018learning). We refer to (bengio2018machine) for a more complete survey. From a more abstract perspective, there has been research on integrating the fields of predictive, prescriptive and descriptive analytics (lombardi2017empirical; bertsimas2019predictive). More directly related to our work, statistical learning has also been applied to the DC Optimal Power Flow Problem (misra2018learning).

In this work, instead of completely replacing existing MIP methods by ML models, as done in previous research with artificial neural networks and genetic algorithms, we propose the usage of ML techniques to significantly enhance the warm-start capabilities of modern MIP solvers when solving SCUC. More specifically, we develop three ML models that can automatically identify different types of patterns in previously solved instances and subsequentially use these patterns to improve MIP performance. The first model is designed to predict, based on statistical data, which constraints are necessary in the formulation and which constraints can be safely omitted. The second model constructs, based on a large number of previously obtained optimal solutions, a (partial) solution which is likely to work well as warm start. The third model identifies, with very high confidence, a smaller-dimensional affine subspace where the optimal solution is likely to lie, leading to the elimination of a large number of decision variables and significantly reducing the complexity of the problem.

We start, in Section 2, by presenting the MIP formulation of SCUC and some brief introduction to the ML concepts we use. In Section 3, we formalize the setting in which the learning takes place and we describe the proposed learning strategies. In Section 4, we evaluate the practical impact of the proposed ML approach by performing extensive computational experiments on a diverse set of large-scale and realistic SCUC instances, containing up to 6515 buses and 1388 units, under realistic uncertainty scenarios. Using the proposed technique, we show that SCUC can be solved on average 4.3x faster than conventional methods with optimality guarantees and 10.2x faster without optimality guarantees, but no observed reduction in solution quality. We also conduct some experiments where the training and test samples are drawn from similar, but not exactly the same distribution, to provide some evidence that the method is somewhat robust against moderate dataset shifts. Finally, in Section LABEL:sec:limitations, we discuss some limitations of the method, as well as directions for future work. Although the techniques are presented in the context of one particular application, they can be easily adapted to other challenging optimization problems.

2 Background

2.1 Mixed-Integer Linear Programming Formulation of SCUC

Unit Commitment (UC) refers to a broad class of optimization problems dealing with the scheduling of power generating units (CohenSherkat87). For each hour in the planning horizon, the problem is to decide which units should be operational and how much power should they produce (OstrowskiAnjosVannelli2012). Each unit typically has minimum and maximum power outputs, a quadratic production cost curve, and fixed startup and shutdown costs. Other constraints typically include ramping rates, which restrict the maximum variation in power production from hour to hour, and minimum-up and minimum-down constraints, which prevents units from starting and shutting down too frequently. The Security-Constrained Unit Commitment (SCUC) problem is a subclass of UC that, in addition to all the constraints previously mentioned, also guarantees the deliverability of power, by enforcing that the power flow in each transmission line does not exceed its safe operational limits (Shaw95). To guarantee sufficient security margin in case of component failure, the problem considers not only transmission constraints under normal operating conditions, but also when there is one transmission line failure in the system (so-called N-1 transmission contingency) (BatutRenaud92). In this way, even if a transmission line unexpectedly fails and the power flow in the network changes, there will be no violations of the transmission line thermal limits, which is a reliability requirement enforced by North American Electric Reliability Corporation (NERCwebsite).

A number of MIP formulations and strong valid inequalities for SCUC have been proposed. In this work, we use the formulation proposed in (MoralesEspanaLatorreRamos2013), since it presents good computational performance even when transmission and N-1 security constraints are present. The full MIP formulation is shown in Appendix A. Here, we present a simplified summary, containing key decision variables and constraints that are most influential to the computational performance of the problem. We note, however, that the methods we present later sections make very few assumptions about the formulation, and can still be used with additional decision variables and constraints.

Consider a power system composed by a set $B$ of buses, a set $G$ of generators and a set $L$ of transmission lines. Let $T$ be the set of hours within the planning horizon. For each generator $g\in G$ and hour $t\in T$, let $x_{gt}$ be a binary decision variable indicating if the generator is operational, and $y_{gt}$ be a continuous decision variable indicating the amount of power being produced. The problem can then be formulated as

[TABLE]

In the objective function, $c_{g}$ is a piecewise-linear function which includes the startup, shutdown and production costs for generator $g$. The notation $(x_{g\bullet},y_{g\bullet})$ is a shorthand for the vector $(x_{gt_{1}},\ldots,x_{gt_{k}},y_{gt_{1}},\ldots,y_{gt_{k}})$, where $\{t_{1},\ldots,t_{k}\}=T$. Constraint (2) enforces a number of constraints which are local for each generator, including production limits, minimum-up and minimum-down running times, ramping rates, among others. Constraint (3) enforces that the total amount of power produced equals the sum of the demands $d_{bt}$ at each bus. Constraints (4) enforces the deliverability of power under normal system conditions. The set $G_{b}$ denotes the set of generators attached to bus $b$. We recall that power flows are governed by physical laws, i.e., Ohm’s laws and Kirchhoff’s Laws. Using the DC linearization of these laws, it is possible to express the thermal limits for each transmission line as (4), where $\delta^{0}_{bl}$ are real numbers known as Injection Shift Factors (ISFs). Similarly, constraints (5) enforce the deliverability of power under the scenario that line $c$ has suffered an outage. Note that the bounds and the ISFs may be different under each outage scenario.

To improve the strength of the linear relaxation, the full formulation contains additional auxiliary binary decision variables. However, once $x_{gt}$ is determined, all the other binary variables are immediately implied. Furthermore, once all the $x_{gt}$ variables are fixed, the problem reduces to a linear programming problem (known as Economic Dispatch) which can be solved efficiently.

Constraints (4) and (5) have a significant impact on the computational performance of SCUC. The total number of such constraints is quadratic on the number of transmission lines. Since large power systems can have more than 10,000 lines, the total number of constraints can easily exceed hundreds of millions. To make matters worse, these constraints are typically very dense, which makes adding them all to the formulation impractical for large-scale instances, not only to due to degraded MIP performance, but even due to memory requirements. Fortunately, it has been observed that enforcing only a very small subset of these constraints is already sufficient to guarantee that all the remaining ones are automatically satisfied (BouffardGalianaArroyo2005). Identifying this critical subset of constraints can be very challenging, and usually involves either solving a relaxation of SCUC multiple times (TejadaSanchezRamos2018; XavierFengWangThimmapuram2018) or solving auxiliary optimization problems (ArdakaniBouffard2015; ZhaiGuanChengWu2010). In practice, system operators still rely on experience and external processes to pre-select a subset of constraints to include in the optimization problem.

2.2 Machine Learning Classifiers

In this work, we employ two classical supervised learning methods for binary classification: $k$-Nearest Neighbors ($k$NN) and Support Vector Machines (SVM). Given a training dataset containing labeled examples $(p^{1},h_{1}),\ldots,(p^{s},h_{s})\in\mathbb{R}^{d}\times\{-1,1\}$, each method constructs a function $\phi:\mathbb{R}^{d}\to\{-1,1\}$ which can be used to label new samples $\tilde{p}\in\mathbb{R}^{d}$. In the following, we present a brief description of each method. For a more complete description, we refer to Alpaydin2014.

$k$NN works by identifying the most similar examples in the training dataset and conducting a simple majority vote. More precisely, given a point $\tilde{p}\in\mathbb{R}^{d}$, the method first sorts the labeled examples by increasing distance $\|p^{i}-\tilde{p}\|$, according to some norm, to obtain a sequence $(p^{i_{1}},h_{i_{1}}),\ldots,(p^{i_{s}},h_{i_{s}})$. Then, it assigns label $1$ to $\tilde{p}$ if $\frac{h_{i_{1}}+\ldots h_{i_{k}}}{k}>0$, and label $-1$ otherwise. More efficient implementations, which avoid sorting the training dataset for each evaluation, have also been described. In Subsections 3.1 and 3.2, we use variations of this method for prediction of redundant constraints and warm starts.

Support Vector Machines construct a hyperplane that best separates training examples into two convex regions, according to their labels, and use this hyperplane for further classification. More precisely, the method first finds a hyperplane $(\pi,\pi_{0})\in\mathbb{R}^{d}\times\mathbb{R}$ by solving the quadratic optimization problem

[TABLE]

where $C$ is a configurable penalty term. Then, it assigns label $1$ to $\tilde{p}$ if $\pi^{T}\tilde{p}\geq\pi_{0}$ and $-1$ otherwise. This description is commonly referred to SVMs with linear kernels. Non-linear extensions have also been described. In Subsection 3.3, we use SVMs to predict affine subspaces.

2.3 Model Evaluation and Cross-Validation

To evaluate the performance of a prediction method on a specific training dataset, we employ $k$-fold cross cross-validation. Given a training dataset $\mathcal{D}=\{(p^{1},h_{1}),\ldots,(p^{s},h_{s})\}$, first we partition $\mathcal{D}$ into $k$ disjoint subsets $\mathcal{D}_{1},\ldots,\mathcal{D}_{k}$ of equal size (we assume that $s$ is divisible by $k$) and train $k$ predictors $\phi_{1},\ldots,\phi_{k}$, using the set $\bigcup_{j\in\{1,\ldots,k\}\setminus\{i\}}\mathcal{D}_{j}$ as training dataset for predictor $\phi_{i}$. The performance of the method on $\mathcal{D}$ is then defined as

[TABLE]

where $f:\{-1,1\}^{d}\times\{-1,1\}^{d}\to[0,1]$ is a performance measure. The idea, for each sample $(p,h)\in\mathcal{D}_{i}$, is to compare the actual label $h$ against the predicted label $\phi_{i}(p)$, produced by a predictor $\phi_{i}$ that did not have access to the pair $(p,h)$ during training. In this work, we use two classical performance measures: precision and recall, given, respectively, by

[TABLE]

3 Setting and Learning Strategies

In this section, we formalize the setting in which the learning takes place, in addition to our proposed learning strategies. The theoretical setting below was designed to match the realistic conditions that energy market operators face on a daily basis.

Assume that every instance $I(p)$ of SCUC is completely described by a vector of real-valued parameters $p$. These parameters include only the data that is unknown to the market operators the day before the market clearing process takes places. For example, they may include the peak system demand and the production costs for each generator. In this work, we assume that network topology and the total number of generators is known ahead of time, so there is no need to encode them via parameters. We may, therefore, assume that $p\in\mathbb{R}^{n}$, where $n$ is fixed. We are given a set of training instances $I(p^{1}),\ldots,I(p^{s})$, which are similar to the instance we expect to solve during actual market clearing. These training instances can either come directly from historical data, if enough similar instances have been solved in the past, or they can be artificially generated, by sampling the probability distribution of $p$. Knowledge about this distribution is not an unreasonable assumption, since market operators have been solving SCUC daily for years and have large amounts of accumulated data which can be readily analyzed.

We also assume that we have at our disposal a customized MIP solver which can receive, in addition to the instance $I(p)$ to be solved, a vector of hints. These hints may affect the performance of the MIP solver and will, hopefully, improve its running time. The hints may also affect the quality or optimality of the solutions produced by the MIP solver. Examples of hints may include a list of warm starts, a list variables that should be fixed, or a list of initial constraints to enforce.

The day before market clearing, a training phase takes place. During this phase, our task is to construct a prediction function $\phi$ which maps each possible parameter $p\in\mathbb{R}^{n}$ into a vector of hints. In order to build this function, we solve the training instances $I(p^{1}),\ldots,I(p^{s})$, capturing any data we may consider useful. During actual market clearing, when we are under time pressure to solve new instances of the problem, the test phase begins. In this phase, one particular vector of parameters $\tilde{p}$ is given, the instance $I(\tilde{p})$ is constructed, and our goal is to solve $I(\tilde{p})$ very quickly. To do this, the vector of hints $\phi(\tilde{p})$ is computed and the pair $(I(\tilde{p}),\phi(\tilde{p}))$ is given to the MIP solver. The total running time is measured, including the time taken to construct the vector of hints from the parameters, as well as the time taken by the MIP solver to solve the instance.

In the following, we propose three different ML models, which target different challenges of solving SCUC. The first model is focused on predicting violated transmission and N-1 security constraints. The second model focuses on producing initial feasible solutions, which can be used as warm starts. The third model aims at predicting an affine subspace where the solution is likely to lie.

3.1 Learning Violated Transmission Constraints

As explained in Subection 2.1, one of the most complicating factors when solving SCUC is handling the large number of transmission and security constraints that need to be enforced. Our first model is designed to predict which transmission constraints should be initially added to the relaxation and which constraints can be safely omitted.

As baseline, we use the iterative contingency screening algorithm presented in (XavierFengWangThimmapuram2018), which was independently implemented and evaluated at MISO, where it presented better performance than internal methods. The method works by initially solving a relaxed version of SCUC where no transmission or N-1 security constraints are enforced. At each iteration, all violated constraints are identified, but only a small and carefully selected subset is added to the relaxation. The method stops when no further violations are detected, in which case the original problem is solved.

In the following, we modify this baseline method to work with hints from a machine learning predictor. Let $L$ be the set of transmission lines. We recall that each transmission constraint takes the form

[TABLE]

We assume that the customized MIP solver accepts a vector of hints

[TABLE]

indicating whether the thermal limits of transmission line $l$, under contingency scenario $c$, at time $t$, is likely to be binding in the optimal solution. Given these hints, the solver proceeds as described in Algorithm 1. In the first iteration, only the transmission constraints specified by the hints are enforced. For the remaining iterations, the solver follows exactly the same strategy as the baseline method, and keeps adding small subsets of violated constraints until no further violations are found. If the predictor is successful and correctly predicts a superset of the binding transmission constraints, then only one iteration is necessary for the algorithm to converge. If the predictor fails to identify some necessary constraints, the algorithm may require additional iterations, but still returns the correct solution. The predictor should still strive to keep a low false positive rate, since the inclusion of too many redundant constraints during the first iteration of the procedure may significantly slow down the MIP optimization time.

For this prediction task, we propose a simple learning strategy based on $k$-Nearest Neighbors. Let $p^{1},\ldots,p^{s}$ be the training parameters. During training, the instances $I(p^{1}),\ldots,I(p^{s})$ are solved using Algorithm 1, where the initial hints are set to RELAX for all constraints. During the solution of each training instance $I(p^{i})$, we record the set of constraints $\mathcal{L}^{i}\subseteq L\times(L\cup\{0\})\times T$ that was added to the relaxation by Algorithm 1. During test, given a vector of parameters $\tilde{p}$, the predictor finds a list of the $k$ vectors $p^{i_{1}},\ldots,p^{i_{k}}$ which are the closest to $\tilde{p}$ and sets $h_{l,c,t}=\texttt{ENFORCE}$ if $(l,c,t)$ appears in at least $p\%$ of the sets $\mathcal{L}^{i_{1}},\ldots,\mathcal{L}^{i_{k}}$, where $p$ is a hyperparameter. When $p=50$, this strategy corresponds to the traditional $k$NN algorithm with simple majority voting. Note that the cost of incorrectly identifying a constraint as necessary (false positive) is typically much smaller than the cost of incorrectly identifying a constraint as redundant (false negative): in the former case, the dimension of the problem simply increases a little, while in the latter, an additional iteration may be necessary. For this reason, in our experiments we try to avoid false negatives (at the cost of more false positives) by setting $p=10$. That is, a constraint is added if is necessary for at least 10% of the nearest training examples.

Note that this $k$NN description generalizes a number of other very simple strategies. For example, when $k=1$, the predictor looks for the most similar instance in the training data. When $k=s$, the strategy is equivalent to counting how many times each transmission constraint was necessary during the solution of the entire training dataset. When $k=1$ and $p=0$, the predictor simply memorizes all constraints that were ever necessary during training. In Subsection 4.3, we run experiments for several of these variations. As we discuss later, these simple $k$NN strategies already performed very well in our experiments, which is the reason we do not recommend more elaborated methods.

We note that this $k$NN transmission predictor is also very suitable for online learning, and can be used even when only a very small number of training samples is available. During the solution of the training instances, for example, after a small number of samples has been solved, subsequent samples can already use hints from this predictor to accelerate the solution process.

3.2 Learning Warm Starts

Another considerably challenging aspect of solving large-scale SCUC instances is obtaining high-quality initial feasible solutions to the problem. Over the last decades, researchers have proposed a number of strong mixed-integer linear programming formulations for SCUC, leading to very strong dual bounds being obtained fairly quickly during the optimization process. For challenging instances, however, we have observed that a significant portion of the running time is usually spent in finding a feasible solution that has an objective value close to this dual bound. These solutions are currently typically found through generic primal heuristics included in the MIP solver, such as feasibility pump (FischettiGloverLodi2005) or large neighborhood search (Shaw1998). Modern MIP solvers also allow users to specify user-constructed solutions, which can be used by the solver to accelerate various parts of the branch-and-cut process. Very importantly, the solutions provided by the user do not need to be complete, or even feasible. Modern solvers can solve restricted subproblems to find values of missing variables, or even repair infeasible solutions using various heuristics. The strategy we present in this subsection takes advantage of these capabilities, instead of trying to replace them. Our aim is to produce, based on a large number of previous optimal solution, a (partial) solution which works well as warm start.

For this prediction task we propose, once again, a strategy based on $k$-Nearest Neighbors. Let $p^{1},\ldots,p^{s}$ be the training parameters and let $(x^{1},y^{1}),\ldots,(x^{s},y^{s})$ be their optimal solutions. In principle, the entire list of solutions $(x^{1},y^{1}),\ldots,(x^{s},y^{s})$ could be provided as warm starts to the solver; however, processing each solution requires a non-negligible amount of time, as we show in Subsection LABEL:subsec:eval_ws. We propose instead the usage of $k$NN to construct a single (partial and possibly infeasible) solution, which should be repaired by the solver. The idea is to set the values only for the variables which we can predict with high confidence, and to let the MIP solver determine the values for all the remaining variables. To maximize the likelihood of our warm starts being useful, we focus on predicting values for the binary variables $x$; all continuous variables $y$ are determined by the solver. Given a test parameter $\tilde{p}$, first we find the $k$ solutions $(x^{i_{1}},y^{i_{1}}),\ldots,(x^{i_{k}},y^{i_{k}})$ corresponding to the $k$ nearest training instances $I(p^{i_{1}}),\ldots,I(p^{i_{k}})$. Then, for each variable $x_{gt}$, we verify the average $\tilde{x}_{gt}$ of the values $x^{i_{1}}_{gt},\ldots,x^{i_{k}}_{gt}$. If $\tilde{x}_{gt}>p$, where $p\in[0.5,1.0]$ is a hyperparameter, then the predictor sets $x_{g,t}$ to one in the warm start solution. If $\tilde{x}_{gt}\leq(1-p)$, then the value is set to zero. In any other case, the value is left undefined.

When $p=0.5$, this strategy corresponds to the traditional $k$NN method, where all variables are defined, even if the neighboring solutions significantly disagree. As $p$ increases, a higher threshold for a consensus is required, and variables for which consensus is not reached are left undefined. More complete warm starts can be processed significantly faster by the solver, but have higher likelihood of being sub-optimal or irreparably infeasible. In Subsection LABEL:subsec:eval_ws, we evaluate several choices of $p$ to quantify this trade-off. In our computational experiments, $p=0.9$ provided the best results.

Like the transmission predictor presented previously, this warm-start predictor is also very suitable for online learning, and can be used with a very small number of training samples.

3.3 Learning Affine Subspaces

Through their experience, market operators have naturally learned a number of patterns in the optimal solutions of their own power systems. For example, they know that certain power generating units (known as base units) will almost certainly be committed throughout the day, while other units (known as peak units) typically come online only during peak demand. These and many other intuitive patterns, however, are completely lost when SCUC is translated into a mathematical problem. Given the fixed characteristics of a particular power system and a realistic parameter distribution, new constraints could probably be safely added to the MIP formulation of SCUC to restrict its solution space without affecting its set of optimal solutions. In this subsection, we develop a ML model for finding such constraints automatically, based on statistical data.

More precisely, let $\tilde{p}$ be a vector of parameters. Our goal is to develop a model which predicts a list of hyperplanes $(h^{1},h^{1}_{0}),\ldots,(h^{k},h^{k}_{0})$ such that, with very high likelihood, the optimal solution $(x,y)$ of $I(\tilde{p})$ satisfies $\langle h^{i},x\rangle=h^{i}_{0}$, for $i=1,\ldots,k$. In this work, we restrict ourselves to hyperplanes that can be written as

[TABLE]

where $g\in G$ and $t\in T$. In other words, the predictor tries to determine whether each variable $x_{gt}$ should be fixed to zero, to one, or to the next variable $x_{g,t+1}$. Furthermore, to prevent conflicting assignments the predictor makes at most one of these three recommendations for each variable $x_{g,t}$.

Let $\mathcal{H}$ be the set of all hyperplanes considered for inclusion. During the training phase, our goal is to construct, for each hyperplane $(h,h_{0})\in\mathcal{H}$, a prediction function $\phi_{(h,h_{0})}:\mathbb{R}^{n}\to\{\texttt{ADD},\texttt{SKIP}\}$, which receives a vector of parameters $\tilde{p}\in\mathbb{R}^{n}$ and returns the label ADD to indicate that the equality constraint $\langle h,x\rangle=h_{0}$ is very likely to be satisfied at the optimal solution, and SKIP otherwise. Given these hints and a vector of parameteres $p$, our customized MIP solver simply adds to the relaxation all the equality constraints $\langle h,x\rangle=h_{0}$ such that $\phi_{(h,h_{0})}(p)=\texttt{ADD}$. These constraints are added at the very beginning of the optimization process and kept until the final solution is obtained; that is, they are never removed. While, in principle, adding these constraints could lead to sub-optimal solutions, in our computational experiments we have observed that, by using a reasonably large number of training samples and by carefully constructing high-precision predictors, as described in detail below, this strategy does not lead to any noticeable degradation in solution quality, while providing significant speedups.

Next we describe how can such high-precision predictors can be constructed in practice. Let $(h,h_{0})\in\mathcal{H}$. Furthermore, let $I(p^{1}),\ldots,I(p^{s})$ be the training instances and let $(x^{1},y^{1}),\ldots,(x^{s},y^{s})$ be their respective optimal solutions. For every $i\in\{1,\ldots,s\}$, let $z^{i}\in\{0,1\}$ be the label indicating whether solution $(x^{i},y^{i})$ lies in the hyperplane $(h,h_{0})$ or not. That is, $z^{i}=1$ if $\langle h,x^{i}\rangle=h_{0}$, and $z^{i}=0$ otherwise. We also denote by $\bar{z}$ the average value of the variables $z^{1},\ldots,z^{s}$.

The proposed algorithm for constructing the function $\phi_{(h,h_{0})}$ is the following. First, if $\bar{z}$ is very close to one, then the predictor always suggests adding this hyperplane to the formulation. More precisely, if $\bar{z}\geq z^{\text{FIX}}$, where $z^{\text{FIX}}$ is a fixed hyper-parameter, then the prediction function $\phi_{(h,h_{0})}$ returns ADD for every input. Next, the the algorithm verifies whether the labels are sufficiently balanced to reliably perform supervised training. More precisely, if $\bar{z}\notin[z^{\text{MIN}},z^{\text{MAX}}]$, where $z^{\text{MIN}}$ and $z^{\text{MAX}}$ are hyper-parameters, then the prediction function $\phi_{(h,h_{0})}$ always returns SKIP. Finally, if $\bar{z}\in[z^{\text{MIN}},z^{\text{MAX}}]$, then the algorithm constructs a binary classifier and evaluates its precision and recall on the training set using $k$-fold cross validation. If the binary classifier proves to have sufficiently high precision and recall, the prediction function $\phi_{(h,h_{0})}$ returns its predictions during the test phase. Otherwise, the algorithm discards the binary classifier, and the function $\phi_{(h,h_{0})}$ returns SKIP for every input. The minimum acceptable recall is given by a fixed hyperparameter $\alpha^{\text{R}}$. The minimum acceptable precision is computed by the expression

[TABLE]

where $\alpha^{\text{P}}$ is a hyper-parameter. Intuitively, the algorithm only accepts the trained binary classifier if it significantly outperforms a dummy classifier which always returns the same label for every input. When $\alpha^{\text{P}}=0$, the trained binary classifier is acceptable as long as it has same precision as the dummy classifier. When $\alpha^{\text{P}}=1$, the classifier is only acceptable if it has perfect precision.

Table 1 shows the precise hyper-parameters used in our experiments, for each type of hyperplane. Due to the high dimensionality of the parameters, we do not train the binary classifiers directly the original vector of parameters $p$. Instead, we propose using as training features only the following subset of parameters: the peak system load, the hourly system loads, the average production cost of generator $g$ and the average production costs of the remaining generators. We propose the usage of support vector machines (with linear kernels), since, in our experiments, these models performed better than logistic regression and random decision forests models, while remaining computationally friendly. Neural networks were not considered given the small number of training samples available.

Unlike the previous two predictors, the affine subspace predictor described in this subsection requires a substantial number of samples in order to give reliable predictions, making it less suitable for online learning.

4 Computational Experiments

In this section we evaluate the practical impact of the ML models introduced in Section 3 by performing extensive computational experiments on a diverse set of large-scale and realistic SCUC instances. In Subsections 4.1 and 4.2, we describe the our implementation of the proposed methods, the computational environment, and the instances used for benchmark. In Subsections 4.3, LABEL:subsec:eval_ws and LABEL:subsec:eval_aff we evaluate, respectively, the performance of the three predictors proposed in Section 3, in addition to several variations. In Subsection LABEL:subsec:ood, we evaluate the out-of-distribution performance of these predictors, to measure their robustness against moderate dataset shift.

4.1 Computational Environment and Instances

The three proposed machine-learning models described in Section 3 were implemented in Julia 1.2 (julialang) and scikit-learn (scikitlearn). Package ScikitLearn.jl (scikitlearnjl) was used to access scikit-learn predictors from Julia, while package DataFrames.jl (dataframesjl) was used for tabular data manipulation. IBM ILOG CPLEX 12.8.0 (cplex) was used as the MIP solver. The code responsible for loading the instances, constructing the MIP, querying the machine learning models and translating the given hints into instructions for CPLEX was also written in Julia, using JuMP (JuMP) as modeling language. Solver features not available through JuMP were accessed through CPLEX’s C API. Training instances were generated and solved at the Bebop cluster at Argonne National Laboratory (Intel Xeon E5-2695v4 2.10GHz, 18 cores, 36 threads, 128GB DDR4). During the training phase, multiple training instances were solved in parallel at a time, on multiple nodes. To accurately measure the running times during the test phase, a dedicated node at the same cluster was used. CPLEX was configured to use a relative MIP gap tolerance of 0.1% during test and 0.01% during training.

A total of nine realistic instances from MATPOWER (matpower) were selected to evaluate the effectiveness of the method. These instances correspond to realistic, large-scale European test systems. Table 2 presents their main characteristics, including the number of buses, units and transmission lines. Some generator data necessary for SCUC was missing in these instances, and was artificially generated based on publicly available data from PjmWebsite and MisoWebsite.

4.2 Training and Test Samples

During training, 300 variations for each of the nine original instances were generated and solved. We considered four sources of uncertainty: i) uncertain production and startup costs; ii) uncertain geographical load distribution; iii) uncertain peak system load; and iv) uncertain temporal load profile. The precise randomization scheme is described below. The specific parameters were chosen based on our analysis of publicly available bid and hourly demand data from PjmWebsite, corresponding to the month of January, 2017.

1. (i)

Production and startup costs. In the original instances, the production cost for each generator $g\in G$ is modelled as a convex piecewise-linear function described by the parameters $c^{0}_{g}$, the cost of producing the minimum amount of power, and $c^{1}_{g},\ldots,c^{k}_{g}$, the costs of producing each additional MW of power within the piecewise interval $k$. In addition, each generator has a startup cost $c^{s}_{g}$. In our data analysis, we observed that the daily changes in bid prices rarely exceeded $\pm 5\%$. Therefore, random numbers $\alpha_{g}$ were independently drawn from the uniform distribution in the interval $[0.95,1.05]$, for each generator $g\in G$, and the costs of $g$ were linearly scaled by this factor. That is, the costs for generator $g$ were set to $\alpha_{g}c^{0}_{g},\alpha_{g}c^{1}_{g},\ldots,\alpha_{g}c^{k}_{g},\alpha_{g}c^{s}_{g}$. 2. (ii)

Geographical load distribution. In the original instances, each bus $b\in B$ is responsible for a certain percentage $d_{b}$ of the total system load. To generate variations of these parameters, random numbers $\beta_{b}$ were independently drawn from the uniform distribution in the interval $[0.90,1.10]$. The percentages $d_{b}$ were then linearly scaled by the $\beta_{b}$ factors and later normalized. More precisely, the demand for each bus $b\in B$ was set to $\frac{\beta_{b}d_{b}}{\sum_{i\in b}\beta_{i}}$. 3. (iii)

Peak system load and temporal load profile. For simplicity, assume $T=\{1,\ldots,24\}$. Let $D_{1},\ldots,D_{24}$ denote the system-wide load during each hour of operation, and let $\gamma_{t}=\frac{D_{t+1}}{D_{t}}$ for $t=1,\ldots,23$, be the hourly variation in system load. In order to generate realistic distribution for these parameters, we analyzed hourly demand data from PjmWebsite. More specifically, we computed the mean $\mu_{t}$ and variance $\sigma^{2}_{t}$ of each parameter $\gamma_{t}$, for $t=1,\ldots,23$. To generate instance variations, random numbers $\gamma^{\prime}_{t}$ were then independently drawn from the Gaussian distribution $\mathcal{N}(\mu_{t},\sigma^{2}_{t})$. Note that the $\gamma^{\prime}$ parameters only specify the variation in system load from hour to hour, and are not sufficient to determine $D_{1},\ldots,D_{24}$. Therefore, in addition to these parameters, a random number $\rho$, corresponding to the peak system load $\max\{D_{1},\ldots,D_{24}\}$, was also generated. In the original instances, the peak system load is always 60% of the total capacity. Based on our data analysis, we observed that the actual peak load rarely deviates more than $\pm 7.5\%$ from the day-ahead forecast. Therefore, to generate instance variations, $\rho$ was sampled from the uniform distribution in the interval $[0.6\times 0.925C,0.6\times 1.075C]$, where $C$ is the total capacity. Note that the $\rho$ and $\gamma^{\prime}$ parameters are now sufficient to construct $D_{1},\ldots,D_{24}$. Figure 1 shows a sample of some artificially generated load profiles.

In this work, we do not consider changes to the network topology, since this information is assumed to be know by the operators the day before the market clearing takes places. If the network topology is uncertain, an additional parameter $p_{l}\in\{0,1\}$ can be added for each transmission line $l\in L$, indicating whether the line is available or not.

During the test phase, 50 additional variations of each original instance were generated. In Subsection 4.3, LABEL:subsec:eval_ws and LABEL:subsec:eval_aff, the test samples were generated using the same randomization scheme outline above, but with a different random seed. In Subsection LABEL:subsec:ood, as described in more detailed in that section, the test samples were generated based on a similar, but different distribution.

4.3 Evaluation of Transmission Predictor

We start by evaluating the performance of different strategies for predicting necessary transmission and N-1 security constraints. For this prediction task, we compared seven different strategies. Strategy zero corresponds to the case where no machine learning is used. In strategy tr:nearest, we add to the initial relaxation all constraints that were necessary during the solution of the most similar training variation. In tr:all, we add all constraints that were ever necessary during training. Strategies tr:knn:$k$ correspond to the method outlined in Subsection 3.1, where $k\in\{10,50,300\}$ is the hyperparameter describing the number of neighbors. We recall that a constraint is initially added to the relaxation it was necessary for at least 10% of the nearest neighbors. Since there are only 300 variations for each instance, strategy tr:knn:300 is equivalent to simply counting, over all variations, how many times each constraint was necessary. Finally, strategy tr:perf shows the performance of the theoretically perfect predictor, which can predict the exact set of constraints found by the iterative method, with perfect accuracy, under zero running time. To implement this predictor, we gave it access to the optimal solution produced by the zero strategy. For each strategy, Table 4.3 shows the average wallclock running time (in seconds), the average number of constraints added to the relaxation per time period, as well as the average number of iterations required by Algorithm 1 to converge. Table LABEL:table:trSpeed shows a more detailed breakdown of the average running times per instance.