Alternating mixed-integer programming and neural network training for
  approximating stochastic two-stage problems
The presented work addresses two-stage stochastic programs (2SPs), a broadly
applicable model to capture optimization problems subject to uncertain
parameters with adjustable decision variables. In case the adjustable or
second-stage variables contain discrete decisions, the corresponding 2SPs are
known to be NP-complete. The standard approach of forming a single-stage
deterministic equivalent problem can be computationally challenging even for
small instances, as the number of variables and constraints scales with the
number of scenarios. To avoid forming a potentially huge MILP problem, we build
upon an approach of approximating the expected value of the second-stage
problem by a neural network (NN) and encoding the resulting NN into the
first-stage problem. The proposed algorithm alternates between optimizing the
first-stage variables and retraining the NN. We demonstrate the value of our
approach with the example of computing operating points in power systems by
showing that the alternating approach provides improved first-stage decisions
and a tighter approximation between the expected objective and its neural
network approximation.
Jan Kronqvist, Boda Li, Jan Rolfes, Shudian Zhao
Alternating mixed-integer programming and neural network training for approximating stochastic two-stage problems+

###### Abstract

The presented work addresses two-stage stochastic programs (2SPs), a broadly applicable model to capture optimization problems subject to uncertain parameters with adjustable decision variables. In case the adjustable or second-stage variables contain discrete decisions, the corresponding 2SPs are known to be NP-complete. The standard approach of forming a single-stage deterministic equivalent problem can be computationally challenging even for small instances, as the number of variables and constraints scales with the number of scenarios. To avoid forming a potentially huge MILP problem, we build upon an approach of approximating the expected value of the second-stage problem by a neural network (NN) and encoding the resulting NN into the first-stage problem. The proposed algorithm alternates between optimizing the first-stage variables and retraining the NN. We demonstrate the value of our approach with the example of computing operating points in power systems by showing that the alternating approach provides improved first-stage decisions and a tighter approximation between the expected objective and its neural network approximation.

Keywords:Stochastic Optimization Neural Network Power Systems.

## 1 Introduction

The area of mathematical optimization often assumes perfect information on the parameters defining the optimization problem. However, in many real-world applications, the input data may be uncertain. The uncertain parameters when optimizing power systems can, e.g., be the power output of wind and solar power units or hourly electricity demands. To address this challenge, two major approaches have been developed in the current literature, namely, robust optimization (RO), see, e.g., [1] for further details and stochastic optimization
(SO), see, e.g., the seminal surveys [3] and [21]. In RO, only limited information on the underlying distribution of uncertainties is available, e.g., that the uncertain parameter vector belongs to a polyhedron, and the goal is to protect against the worst possible outcome of the uncertainty. In SO, it is typically assumed that the underlying probability distribution of the uncertain parameters is known, or at least that representative samples are available. Here, often the goal is to optimize over the expectation of the uncertain parameters, which results in less conservative solutions compared to RO. The ability to handle uncertainties makes both methods valuable in power-system-related optimization. A comparison between RO and SO for unit commitment with uncertain renewable production is presented in [13]. However, SO tends to be more prevalent in power system operation optimization, as operating scenarios are readily available. Applications of SO include the optimal operation of high-proportion renewable energy systems [23, 32, 38], pricing of energy and reserves in electricity markets [35, 7, 39] or day-ahead power scheduling [36, 30, 24] - to name a few.

Moreover, in some real-world applications, not all the decisions have to be taken before the uncertainty realizes. Here, we distinguish between _first-stage_ or here-and-now decisions and _second-stage_ or wait-and-see decisions. These _stochastic two-stage_ problems or _2SP_s have drawn significant attention over the last years, both due to the practical applicability and computational challenges. Methods for dealing with 2SPs include L-shaped methods [20], Lagrangian relaxation-based methods [6], and the sample average approach of Kleywegt et al. [15]. For more information on 2SP and algorithmic approaches, we refer to [33, 19, 8] and the references therein.

Here, we focus on 2SPs, where the uncertainties either belong to a bounded discrete distribution or can be approximated sufficiently close by such a distribution. Consequently, we can restrict ourselves to a finite set of scenarios. Furthermore, we are assuming some of the second-stage variables are restricted to integer values and that all constraints are linear. For such problems, it is possible to represent the 2SP as a classical single-stage mixed-integer linear program (MILP) through the so-called deterministic equivalent [3]. The deterministic equivalent formulation creates one copy of the second-stage variables and constraints for each scenario, thus resulting in a potentially huge MILP problem. Especially if the second stage involves many integer variables, it may be computationally intractable to consider a large number of scenarios. To overcome this issue, Dumouchelle et al. [8] proposed the so-called Neur2SP approach,where the value function, i.e., the optimal objective value of the second stage, is approximated by a neural network (NN). By selecting a MILP representable NN architecture, e.g., ReLU activation functions [22, 9, 34], the NN approximation of the second-stage problem along with the first-stage problem can then be integrated into a single stage MILP. The advantage of this approach is that the MILP presentation of the NN can require far fewer variables and constraints than forming the deterministic equivalent. The potential drawback is that we only utilize an approximation of the expected value of the second-stage problem.
We build upon the approach presented by Dumouchelle et al. [8] and investigate the proposed algorithm for an application in optimizing power systems. We show that proper training of the NN is crucial, and we propose a dynamic sampling as a training method. Here, we directly train a NN to approximate the expected value function. Our approach can be viewed as a surrogate-based optimization [2, 10] approach for 2SPs, as we are iteratively forming a surrogate of the expected value function which is globally optimized to find a new solution and to improve the approximated expected value function. The proposed algorithm is not limited to the specific application in Section 3, but applicable to any 2SP of suitable structure, e.g., see the applications used in [8].

For the application of finding optimal operating points of a power grid considered in the present article, the existing approaches in the power system literature usually approximate the integer variables in the second stage, e.g. by ignoring the action variable of energy storage, see [5] or only search for a local optimum of the model, e.g. by applying an augmented Benders decomposition algorithm, see [4]. Thus, the Neur2SP approach keeps the underlying integral structure of the second-level, while at the same time aims for a global optimum.

The outline of the present article is as follows. In Section 2, we present the stochastic programming formulation [8] and discuss our methodology. Then, we connect this modeling to achieve data-driven operating states of smart converters in power systems in Section 3. Subsequently, we present a numerical comparison between these data-driven results with a baseline algorithm based on the Neur2SP framework in Section 4.

## 2 Replacing the second-stage in stochastic programming by a neural network

The two-stage stochastic programming consists of two types or _levels_ of decision variables. The _first-level_ variables, denoted by \(x\in\mathcal{X}\), describe an initial planning approach. Then, a random vector \(\xi\in\Omega\) supported on a compact domain \(\Omega\subseteq\mathbb{R}^{I}\) and distributed by a probability distribution \(\mathbb{P}\in\mathcal{P}(\Omega)\) affects the outcome of this initial planning. However, one can adjust the initial planning \(x\) after the random vector \(\xi\) realizes by choosing the _second-level or recourse_ variables \(y\in\mathcal{Y}(x,\xi)\), wherein the present article, we suppose that \(\mathcal{Y}(x,\xi)\subseteq\mathbb{R}^{n}\times\{0,1\}^{l}\) is assumed to be defined through linear constraints. Thus, we can summarize a _stochastic two-stage problem_ as follows:

\[\min_{x\in\mathcal{X}}\;G(x)+\mathbb{E}_{\mathbb{P}}\left(\min_{y\in\mathcal{Y }(x,\xi)}c^{\top}y\right), \tag{1}\]

where in this article \(\mathcal{X}\) and \(\Omega\) are assumed to be polytopes, \(c\) denotes the objective coefficients, and \(\mathcal{Y}(x,\xi)\) a polytope intersected with an integer lattice. Note, that (1) is considered to be a challenging problem as it contains the NP-complete MIP \(\min_{y\in\mathcal{Y}(x,\xi)}c^{\top}y\) as a subproblem. Moreover, stochastic two-stage problems are in general notoriously challenging, see e.g., Section 3 in [21] for an overview.
In the present article, we restrict ourselves to the estimated expected value

\[\min_{x\in\mathcal{X}}\ G(x)+\frac{1}{m}\sum_{j=1}^{m}\left(\min_{y\in\mathcal{Y}( x,\xi_{j})}c^{\top}y\right), \tag{2}\]

where \(\xi_{j}\) denote scenarios sampled from \(\mathbb{P}\). For the numerical experiments in Section 4, we limit the number of scenarios to \(m=400\). We then aim to learn the mapping

\[Q(x)\coloneqq\frac{1}{m}\sum_{j=1}^{m}\left(\min_{y\in\mathcal{Y}(x,\xi_{j})}c ^{\top}y\right)\]

by training a series of neural networks inspired by the Neur2SP framework developed in [8]. However, in contrast to the approach in [8], we aim to learn the mapping \(x\mapsto Q(x)\), by training a series of neural networks with iteratively generated training data instead of predefined data, see Algorithm 1.

### Embedding the neural network into 2SP

We approximate the map \(x\to Q(x)\) by training a fully-connected neural network with 2-hidden layers, where both layers have 40 neurons, denoted as NN \(2\times 40\). Hence, the relation between \(l\)-th layer input \(x^{l}\in\mathbb{R}^{n_{l}}\) and output \(x^{l+1}\in\mathbb{R}^{n_{l+1}}\) in such neural network with the ReLu activation function is

\[x^{l+1}=\max\{0,W^{l}x^{l}+b^{l}\}, \tag{3}\]

where \(W^{l}\in\mathbb{R}^{n_{l+1}\times n_{l}}\) is the weight matrix and \(b^{l}\in\mathbb{R}^{n_{l+1}}\) is the bias vector. The ReLu activation function is piece-wise linear and the big-M formulation was the first presented to encode it and is still common [9, 22]. In this way, with binary variable \(\sigma^{l}\in\mathbb{R}^{n_{l+1}}\), (3) is equivalent to

\[\begin{split}(w_{i}^{l})^{\top}x^{l}+b_{i}^{l}& \leq x_{i}^{l+1},\\ (w_{i}^{l})^{\top}x^{l}+b_{i}^{l}-(1-\sigma_{i}^{l})LB_{i}^{l+1}& \geq x_{i}^{l+1},\\ x_{i}^{l+1}&\leq\sigma_{i}^{l}UB_{i}^{l+1},\\ \sigma_{i}^{l}&\in\{0,1\},\ x_{i}^{l+1}\geq 0, \forall i\in[n_{l+1}],\end{split} \tag{4}\]

where \(w_{i}^{l}\) is the \(i\)-th row vector of \(W^{l}\), \(b_{i}^{l}\) is the \(i\)-th entry of \(b^{l}\), \(LB_{i}^{l+1}\) and \(UB_{i}^{l+1}\) are upper and lower bounds on the pre-activation function over \(x_{i}^{l+1}\), such that \(LB_{i}^{l+1}\leq(w_{i}^{l})^{\top}x^{l}+b_{i}^{l}\leq UB_{i}^{l+1}\). Encoding ReLu functions has been an active research topic in recent years, as it enables a wide range of applications of mixed-integer programming to analyze and enhance NNs, such as generating adversarial examples [9], selecting optimal inputs of the training data [18, 40], lossless compression [29], and robust training [28].

By encoding the trained NN as (4), we can approximate (1) by the following MIP:
\[\min G(x^{1})+x^{L}\] (5a) s.t. \[W^{l}x^{l}+b^{l}\leq x^{l+1}, \forall\ l\in[L-1], \tag{5b}\] \[W^{l}x^{l}+b^{l}-\mathrm{diag}(LB^{l+1})(\mathbf{1}-\sigma^{l+1}) \geq x^{l+1}, \forall\ l\in[L-1],\] (5c) \[x^{l}\leq\mathrm{diag}(UB^{l})\sigma^{l}, \forall\ l\in\{2,\ldots,L\},\] (5d) \[x^{L}=W^{L-1}x^{L-1}+b^{L-1},\] (5e) \[\sigma^{l}\in\{0,1\}^{n_{l}}, \forall\ l\in\{2,\ldots,L\},\] (5f) \[x^{l}\in\mathbb{R}_{+}^{n_{l}}, \forall\ l\in[L-1],\] (5g) \[x^{1}\in\mathcal{X},x^{L}\in\mathbb{R}, \tag{5h}\]

where \(x^{1}\coloneqq x\) is the input-layer variable that belongs to the polytope \(\mathcal{X}\), and \(x^{L}\) is the output of the neural network which approximates the expected value \(Q(x^{1})\).

### Alternating MIP and NN training

Since the approximation quality of (5) with regards to (2) highly depends on the accuracy of the neural network, the approach has significant drawbacks for problems with high dimensional first-level decisions, i.e., high-dimensional \(\mathcal{X}\). This drawback becomes even more significant as optimal solutions for (2) are often attained at the boundary or even vertices of \(\mathcal{X}\). Here, the approximation by the neural network tends to be worse due to a lack of training data around the vertices. Note that, sampling around each vertex is typically computationally intractable as there are potentially exponentially many such vertices.

The approach we propose in the present article uses intermediate solutions of (5) to inform the retraining of neural networks. To this end, after training the NN with an initial uniformly distributed sample of \(\mathcal{X}\), we solve (5) and add sample points around the computed optimal solution \(x_{k}^{*}\) to our training data. Algorithm 1 illustrates our approach.

We would like to stress, that the number of additional datapoints in every iteration \(k\) as well as the parameter \(\alpha\) towards \(x_{k}^{*}\) are chosen arbitrarily. Thus, more careful choices of these parameters may give room for significant improvement of the performance of Algorithm 1.

Here, we would like to briefly comment on the relative computational expenses of the different parts of Algorithm 1. We specify the considered application as well as our computational ressources in Section 4. The most time-consuming computation of Algorithm 1 includes solving optimization problems (5) and \(\min_{y\in\mathcal{Y}(x,\xi)}c^{\top}y\). With the default settings, Gurobi [11] can solve (5) encoding an NN \(2\times 40\) within 5 seconds on a standard notebook. Moreover, the generation of new training data (see line 6-10), i.e., the optimization for each scenario \(\xi\), can be solved in parallel, where solving the problem \(\min_{y\in\mathcal{Y}(x,\xi)}c^{\top}y\) with a pair \(x\) and \(\xi\) only takes 0.03 seconds.
```
1:Input Initial training data \(\tilde{\mathcal{X}}\), \(\alpha=0.99\)
2:Output Approximately optimal solution \(x_{100}^{*}\) for (1)
3: Train a NN \(2\times 40\) with training data set \(\tilde{\mathcal{X}}\)
4:for\(k=1,\ldots,100\)do
5:\(x_{k}^{*}\leftarrow\operatorname*{argmin}\eqref{eq:x_k}\)
6:for\(i=1,\ldots,50\)do
7: Sample \(x_{i}\in X\) uniformly
8:\(x_{i}\leftarrow\alpha x_{k}^{*}+(1-\alpha)x_{i}\)
9:\(\tilde{\mathcal{X}}\leftarrow\tilde{\mathcal{X}}\cup\{x_{i}\}\)
10:endfor
11: Retrain NN with \(\tilde{\mathcal{X}}\)
12:endfor
```

**Algorithm 1** Alternating MIP-NN Algorithm for 2SP (MIP-NN 2SP)

## 3 Application to smart converters in power system networks

In order to assess the practical value of the above approach, we apply it to a power system with smart inverters and energy storage. Its effects include reducing power generation costs and improving operational performance are then observed. Since the proposed heuristic given by Algorithm 1 may compute first-stage decisions \(x\) that lead to an empty second-stage, i.e., \(\mathcal{Y}(x,\xi)=\emptyset\), we assume in the remainder of this article that every feasible \(x\in\mathcal{X}\) leads to non-empty \(\mathcal{Y}(x,\xi)\) for every scenario \(\xi\). In particular, this assumption is valid whenever, we discuss a value function \(Q:\mathcal{X}\times\Omega\to\mathbb{R},Q(x,\xi)\coloneqq\min_{y\in\mathcal{Y} (x,\xi)}c^{\top}y\) as is often the case in the existing literature, see e.g. [8], and also justified as the discussion in Section 2.1 in [31] illustrates.

Nevertheless, it is crucial to verify whether this assumption actually holds. For the power system application studied below however, numerical experiments indicate, that a feasible, but potentially costly \(y\in\mathcal{Y}(x,\xi)\), exists and thus the use of neural networks will not have an adverse impact on the feasibility of the problem. Moreover, from an engineering perspective, the power system has a slack node, see \(b_{1}\) in Figure 1, that ensures stable operation when the system deviates from the set points, although again the operating cost may increase significantly. On the contrary, it can effectively improve the computational efficiency of the solution. The reasons are as follows:

1. This paper establishes a two-stage optimization model to solve the day-ahead scheduling strategy of the power system (i.e., the first-stage strategy). We only use neural networks to replace the optimization model in the second stage. The first-stage decisions still satisfy the corresponding hard constraints (e.g., power balance).
2. The second-stage optimization model corresponds to the intra-day operation of the power system. In the cases studied by this paper, the system is connected to the external grid and has energy storage inside as a buffer to
balance the energy. From the perspective of practical engineering experience, when the decisions of the first stage are within a reasonable range, the operation strategy of the second stage is usually feasible. Therefore, it is reasonable to use the first-stage decision variables to estimate the operating cost of the second stage through a neural network.
3. In practical operations, the day-ahead operation (i.e., the first stage) needs to obtain results within tens of minutes. When there are a large number of integer variables in the second stage of the 2SP (energy storage requires frequent changes in operational states, which introduces a large number of binary variables into the second stage), this runtime requirement is difficult to guarantee. However, since replacing the second stage of (2) by a neural network can significantly reduce the runtime, Algorithm 1 has the potential to be applicable in practice.

A simplified DC power flow model, following Kirchhoff laws, is adopted to describe the energy balance and power flow in the power system. We follow closely the notations by [37] and the presentation in [16] in order to specify the parameters present in (1).

Consider a power grid equipped with a set of buses \(\mathcal{B}\) and a set of lines/branches \(\mathcal{L}\). Within the grid, power is generated by a set of conventional (fossil fuel) generators denoted by \(\mathcal{N}_{G}\) or by a set of distributed (renewable) generators denoted by \(\mathcal{N}_{DG}\). The _system operator_ (SysO) may decide, whether to

* store or release power to/from a set of batteries \(\mathcal{N}_{S}\)
* purchase or sell power on the day-ahead market or intra-day, this power enters/leaves the power system through a trading node with the main grid

The amount of power traded day-ahead, i.e., on the day-ahead or _first-level market_ is denoted by \(P_{fl}\) and the corresponding market price by \(p_{fl}\). Similarly, the amount of power and the market price traded intra-day is denoted by \(P_{sl}\) and \(p_{sl}\), respectively. A positive value for \(P_{fl},P_{sl}\) is interpreted as a purchase, negative values as a sell \(P_{fl},P_{sl}\) of energy.

The SO's initial planning, refers to deciding the first-level variables \(x=(P_{G}^{\top},P_{fl})^{\top}\in\mathbb{R}^{\mathcal{N}_{G}}\times\mathbb{R}\) based on the estimated renewable energy production \(P_{DG,forecast}\in\mathbb{R}^{\mathcal{N}_{DG}}\) in order to ensure that a given total demand \(\sum_{i\in\mathcal{B}}P_{d_{i}}\) is met. Hence,

\[\mathcal{X}=\{x\in\mathbb{R}^{n}:\ \eqref{eq:pfl}\ \&\ \eqref{eq:pfl}\},\]

where

\[P_{fl}^{t}=\sum_{i\in\mathcal{B}}P_{d_{i}}^{t}-\sum_{i\in \mathcal{N}_{G}}P_{G_{i}}^{t}-\sum_{i\in\mathcal{N}_{DG}}P_{DG_{i},forecast} \forall t\in T, \tag{6a}\] \[P_{G_{i},\min}^{t}\leq P_{G_{i}}^{t}\leq P_{G_{i},\max} \forall i\in\mathcal{N}_{G},t\in T, \tag{6b}\]

with given parameters \(P_{G,\min},P_{G,\max}\in\mathbb{R}^{\mathcal{N}_{G}\times T}\). Note that the _market-clearing condition_ (6a) ensures the active power balance in the whole system. The objective
of the first level in (1) is then given by

\[G(x)=x^{\top}\text{Diag}(c_{2})x+c_{1}^{\top}x+c_{0},\]

where \(c_{2},c_{1}\in\mathbb{R}^{\mathcal{N}_{G}},c_{0}\in\mathbb{R}\) are given generator cost parameters. Since the uncertainties will impact the initial planning and may jeopardize the power balance, we consider the capacity of the renewable generators as uncertain, i.e., we denote the second-level variable by \(\xi=P_{DG,\max}\in\mathbb{R}^{\mathcal{N}_{DG}\times T}\). As these generators are dependent on weather conditions, which are highly uncertain, this is one of the most common uncertainties faced by modern power grids with a high proportion of renewable energies [27, 12]. To this end, we draw samples from the following domain

\[\Omega=\left\{P_{DG,\max}\in\mathbb{R}^{\mathcal{N}_{DG}\times T}:\ 0\leq P_{DG_{i},\max}^{t}\leq P_{i}^{+}\ \forall\ i\in\mathcal{N}_{DG},t\in T\right\}, \tag{7}\]

where \(P_{i}^{+}\) denotes the technical limit of the renewable generator, i.e. its capacity under optimal conditions. In particular, since the capacity of the renewable generators is further limited by the weather conditions.

For our computational results in Section 4, we assume that without these limitations \(P_{DG_{i},\max}\) are independently distributed according to a normal distribution, i.e. \(P_{DG_{i},\max\text{-}\text{vol}}\sim\mathcal{N}(P_{DG_{i},forecast},0.1\cdot P _{DG_{i},forecast})\). Consequently, every realization of \(P_{DG_{i},\max\text{-}\text{vol}}\) leads naturally to a sample point

\[\xi_{DG_{i}}=P_{DG_{i},\max}=\begin{cases}0&\text{if }P_{DG_{i},\max\text{-} \text{vol}}\leq 0\\ P_{DG_{i},\max\text{-}\text{vol}}&\text{if }P_{DG_{i},\max\text{-}\text{vol}} \in(0,P_{i}^{+})\\ P_{i}^{+}&\text{if }P_{DG_{i},\max\text{-}\text{vol}}\geq P_{i}^{+}.\end{cases}\]

On the third level, the SysO adjusts the initial planning according to this realization. In particular, the SysO might regulate the energy output of the conventional generators \(P_{G}\) to \(P_{G,\text{reg}}\) by either increasing the production by adding \(P_{G}^{+}\geq 0\) at a cost \(r^{+}\) or decreasing the production by adding \(P_{G}^{-}\leq 0\) at a cost \(r^{-}\). Similarly, \(P_{DG}\in[0,\xi_{DG}]\) is the adjusted energy production that deviates from its forecast by \(P_{DG}^{+}\) or \(P_{DG}^{-}\) with deviations penalized by \(f^{+},f^{-}\) respectively. Moreover, the SysO might also trade power intra-day (\(P_{sl}\)) or decide to use the batteries by setting \((\mu_{ch},\mu_{dch}\in\{0,1\}^{\mathcal{N}_{S}})\) and change the charging/discharging quantity \(P_{ch},P_{dch}\) in order to balance a potential power deficiency or surplus. Based on these decisions, the following variables vary accordingly: The _state of charge_ of a storage, denoted by soc, the power on a line \((k,l)\in\mathcal{L}\), denoted by \(p_{kl}\) and the phase angles of the system, denoted by \(\theta\). Thus, the SysO adjusts the vector \(y=(P_{G_{i},\text{reg}},P_{G_{i}}^{+},P_{G_{i}}^{-},P_{DG_{i}},P_{DG_{i}}^{+}, P_{DG_{i}}^{-},P_{sl},P_{ch_{i}},P_{dch_{i}},p_{kl},\theta,\text{soc},\mu_{ ch},\mu_{dch})^{\top}\) in order to satisfy

\[y\in\mathcal{Y}(x,h)\coloneqq\left\{y\in\mathbb{R}^{m}:\ \text{(\ref{eq:p_1})}- \text{(\ref{eq:p_2})}\right\},\]

where the constraints (8a)-(12g) are given below:
1. First, we consider the _Generator and DG output constraints_: \[P^{t}_{G_{i},\mathrm{reg}}=P^{t}_{G_{i}}+P^{t,+}_{G_{i}}+P^{t,-}_{G_ {i}} \forall i\in\mathcal{N}_{G},t\in T,\] (8a) \[P^{t}_{G_{i},\min}\leq P^{t}_{G_{i},\mathrm{reg}}\leq P^{t}_{G_{i },\max} \forall i\in\mathcal{N}_{G},t\in T,\] (8b) \[P^{t}_{DG_{i}}=P^{t}_{DG_{i},\mathrm{forecast}}+P^{t,+}_{DG_{i}}+ P^{t,-}_{DG_{i}} \forall i\in\mathcal{N}_{DG},t\in T,\] (8c) \[P^{t}_{DG_{i},\min}\leq P^{t}_{DG_{i}}\leq P^{t}_{DG_{i},\max} \forall i\in\mathcal{N}_{DG},t\in T,\] (8d) \[P^{t}_{sl}-P^{t}_{fl}=-\mathbb{1}^{\top}(P^{t,+}_{G}+P^{t,-}_{G}) -\mathbb{1}^{\top}(P^{t,+}_{DG}+P^{t,-}_{DG})\] \[-\mathbb{1}^{\top}(P^{t}_{dch}-P^{t}_{ch}) \forall t\in T.\] (8e) Here, Constraints (8a)-(8d) describe the output range of the conventional and renewable generators. In particular, (8d) shows that the third-level variables \(P^{t}_{DG_{i}}\) are restricted by the uncertainties realized in the second level (\(P^{t}_{DG_{i},\max}\)). The market clearing on the intra-day market, i.e., the actual power demand-supply relations, is reflected by Constraint (8e).
2. Second, we consider the _operation constraints_: \[0\leq P^{t,+}_{G_{i}}\leq P^{t,+}_{G_{i},\max} \forall i\in\mathcal{N}_{G},t\in T,\] (9a) \[P^{t,-}_{G_{i},\min}\leq P^{t,-}_{G_{i}}\leq 0 \forall i\in\mathcal{N}_{G},t\in T,\] (9b) \[P^{t,+}_{DG_{i}}\geq 0 \forall i\in\mathcal{N}_{DG},t\in T,\] (9c) \[P^{t,-}_{DG_{i}}\leq 0 \forall i\in\mathcal{N}_{DG},t\in T,\] (9d) where Constraints (9a) -(9d) limit the real output derivations of conventional and renewable generators.
3. Third, we consider the _power flow constraints_. To this end, we apply a DC approximation for a given line reactance (\(x_{ij}>0\)) and demand in active power (\(P^{d,t}_{k}\)) at every time step \(t\) and bus \(k\): \[\sum_{l\in\delta(k)}p^{t}_{kl}= \sum_{i\in\mathcal{N}_{G}:i\sim k}P^{t}_{G_{i},\mathrm{reg}}+\sum _{i\in\mathcal{N}_{DG}:i\sim k}P^{t}_{DG_{i}}\] \[+\sum_{i\in\mathcal{N}_{S}:i\sim k}(P^{t}_{dch_{i}}-P^{t}_{ch_{i} })-P^{d,t}_{k} \forall k\in\mathcal{B}\setminus\{0\},t\in T,\] (10a) \[\sum_{l\in\delta(k)}p^{t}_{0l}= \sum_{i\in\mathcal{N}_{G}:i\sim 0}P^{t}_{G_{i},\mathrm{reg}}+\sum _{i\in\mathcal{N}_{DG}:i\sim 0}P^{t}_{DG_{i}}\] \[+\sum_{i\in\mathcal{N}_{S}:i\sim 0}(P^{t}_{dch_{i}}-P^{t}_{ch_{i} })+P^{t}_{sl}-P^{d,t}_{0} \forall t\in T,\] (10b) \[p^{t}_{ij}= \frac{1}{x_{ij}}(\theta^{t}_{i}-\theta^{t}_{j}) \forall\{i,j\}\in\mathcal{L},t\in T,\] (10c) where (10a) establishes the nodal power flow balance for every node except the root node, which is addressed separately in (10b). The branch power flow is modeled through (10c).
4. Fourth, we consider the _branch thermal constraints_ that guarantee, that the power flow does not exceed the branch's capacities: \[p^{t}_{ij}\leq s_{ij,\max} \forall\{i,j\}\in\mathcal{L},t\in T.\] (11a)
5. Fifth, we consider the _storage constraints_. To this end, note that the storage operation involves two actions, the storage action is represented by two binary variables \(\mu_{ch},\mu_{dch}\in\{0,1\}^{\mathcal{N}_{S}\times T}\), whereas the quantity of the charging/discharging is represented by continuous variables \(P^{t}_{ch_{i}},P^{t}_{dch_{i}}\). Thus, we have the following constraints: \[\operatorname{soc}^{t}_{i,\min}\leq\operatorname{soc}^{t}_{i} \leq\operatorname{soc}^{t}_{i,\max} \forall i\in\mathcal{N}_{S},t\in T,\] (12a) \[\operatorname{soc}^{t}_{i}=\operatorname{soc}^{t-1}_{i}+\frac{( P^{t}_{ch_{i}}-P^{t}_{dch_{i}})}{E_{i}}\Delta T \forall i\in\mathcal{N}_{S},t\in T,\] (12b) \[P^{t}_{ch_{i}},P^{t}_{dch_{i}}\geq 0 \forall i\in\mathcal{N}_{S},t\in T,\] (12c) \[\mu^{t}_{ch_{i}},\mu^{t}_{dch_{i}}\in\{0,1\} \forall i\in\mathcal{N}_{S},t\in T,\] (12d) \[\mu^{t}_{ch_{i}}P^{t}_{ch_{i},\min}\leq P^{t}_{ch_{i}}\leq\mu^{t} _{ch_{i}}P^{t}_{ch_{i},\max} \forall i\in\mathcal{N}_{S},t\in T,\] (12e) \[\mu^{t}_{dch_{i}}P^{t}_{dch_{i},\min}\leq P^{t}_{dch_{i}}\leq\mu^{t }_{dch_{i}}P^{t}_{dch_{i},\max} \forall i\in\mathcal{N}_{S},t\in T,\] (12f) \[\mu^{t}_{ch_{i}}+\mu^{t}_{dch_{i}}\leq 1 \forall i\in\mathcal{N}_{S},t\in T.\] (12g)

Here, the upper/lower bounds of \(\operatorname{soc}\), as well as the relationships between \(\operatorname{soc}\) and charging/discharging actions, are given by Constraints (12a) and (12b). The connection between the storage actions and the respective quantity is reflected by Constraints (12c) - (12g).

Lastly, the third-level cost function aims to both, minimize the electricity cost as well as reduce the deviation between the intra-day system operation strategy and the day-ahead planning. Thus, the whole adjustment can be summarized as solving

\[\min_{y\in\mathcal{Y}(x,\xi)}c^{\top}y, \tag{13}\]

where

\[c^{\top}y\coloneqq\sum_{t\in T}\sum_{i\in\mathcal{N}_{G}}(r^{+,t}_{i}P^{+,t}_ {G_{i}}+r^{-,t}_{i}P^{-,t}_{G_{i}})+p^{t}_{sl}(P^{t}_{sl}-P^{t}_{fl})+\sum_{i \in\mathcal{N}_{DG}}(f^{+}_{i}P^{t,+}_{DG_{i}}+f^{-}_{i}P^{t,-}_{DG_{i}}).\]

After having established the model parameters for (1), we continue with a case study in order to demonstrate the practical value of Algorithm 1.

## 4 Computational results

In this section, we present a case study, where we consider the daily power distribution of a 5-bus instance, based on the "case5.m" instance from the matpower library [41]. Here, the day is divided into hourly (24 period) time intervals. The MIP problems in Algorithm 1 are solved by Gurobi [11], which is
the state-of-art MIP solver. The computations were executed via Gurobi 10.0.0 under Python 3.7 on a Macbook Pro (2019) notebook with an Intel Core i7 2,8 GHz Quad-core and 16 GB of RAM. The neural network is trained with PyTorch [25].

The network topology of "case5.m" is illustrated in Figure 1, where we assume bus 1 to be the slack node, i.e., the trading node connected to other grids. Since the type of generators is not specified in [41] and this study solely serves as an academic example, we chose whether a generator in "case5.m" is a conventional/renewable one or a storage in the following convenient way: The conventional (fossil fuel) generators are connected to the buses 1 and 4, i.e. \(\mathcal{N}_{G}=\{1,4\}\), two distributed generators are connected to buses 1 and 5, i.e., \(\mathcal{N}_{DG}=\{1,5\}\) and an energy storage unit is connected to bus 3, i.e. \(\mathcal{N}_{S}=\{3\}\). In order to aid reproducibility and encourage further research on this topic, we provide the underlying system data for public use, see [17].

Incorporating the instance data closely follows the methodology used in [16]. Hence, both, the day-ahead and intra-day market prices \(p_{fl},p_{sl}\) were taken as averages from the Pecan street database's [26] "miso" data set for March 3rd, June 3rd, October 3rd and December 2nd, 2022. Daily deviations in \(P^{d}\), denoted by \(\Delta^{t}_{d}\), or daily deviations in \(P^{t}_{DG,\max}\), denoted by \(\Delta^{t}_{DG}\) were taken from the Pecan street database's "california_iso" dataset for the same days, 2022 in the same way. In particular, we also provide the underlying weather data for public use, see [17].

Figure 1: The case5.m network with its corresponding generators and an exemplary power flow. \(P^{d}=0\) at buses 1 and 5
Consequently, we model the daily varying demand by \(P^{d}=P^{d}\cdot\Delta_{d}^{t}\) and the varying potential renewable energy production by

\[P_{DG,\text{forecast}}^{t}\coloneqq\min\left\{\frac{P_{DG,\min}+P_{DG,\max}}{2} \cdot\Delta_{DG}^{t},P_{i,t}^{+}\right\}.\]

After incorporating this data, we compare the solutions given by Algorithm 1 to a baseline experiment, where we replace line 6-10 by adding 50 uniformly distributed sample points drawn from \(\mathcal{X}\) in each iteration \(k\). This algorithm is inspired by [8] and simply creates the same amount of sample points, but does not incorporate targeted information around \(x_{k}^{*}\). We chose 3000 uniformly distributed datapoints from \(\mathcal{X}\) as a starting dataset for both algorithms. We use trained neural networks with the same architecture for both experiments, i.e., NN \(2\times 40\).

In Figure 2, MIP-NN_2SP_NN (resp. Baseline_NN) denotes the objective function value of (5) solved by Algorithm 1 (resp. the baseline experiment) at

Figure 2: Comparison of Algorithm 1 (in green and red) to successively adding batches of 50 uniformly distributed data points (in blue and orange) with respect to operating costs in $ for weather data of different seasons.
each iteration. In particular, we include \(G(x_{k}^{*})+\frac{1}{400}\sum_{j=1}^{400}Q(x_{k}^{*},\xi_{j})\), denoted by MIP-NN_2SP_E (resp. Baseline_E), with optimal solutions \(x_{k}^{*}\) of line 5 from Algorithm 1 (resp. the baseline experiment). In this way, the difference between MIP-NN_2SP_NN (resp. Baseline_NN) and MIP-NN_2SP_E (resp. Baseline_E) measures the quality of the approximation by the neural networks at each iteration.

Observe, that after \(100\) iterations, the baseline experiment cannot close the gap between the estimated expected value of the second stage (i.e., Baseline_E) and its neural network approximation (i.e., Baseline_NN) since a gap of \(80,000\$\) remains. In other words, the difference between the NN approximation and the expected value is \(80,000\$\). While including targeted datapoints in the proposed sampling approach significantly reduces the gap between MIP-NN_2SP_NN and MIP-NN_2SP_E. Thus, showing that the dynamic sampling in Algorithm 1 clearly outperforms sampling uniformly distributed and results in a much more accurate approximation of the expected value. Moreover, the results indicate, that the operating state of the conventional generators computed from the baseline experiment is not optimal, and approximately \(20,000\$\) may be saved by using the operating states computed by Algorithm 1. The results, thus, show the importance of obtaining an accurate approximation of the expected value function.

## 5 Conclusion and Outlook

We have presented an algorithm for solving 2SP problems with integer second-stage variables by iteratively constructing a NN that approximates the expected value of the second-stage problem. This approach is inspired by the work of [8]. In the algorithm, we propose a dynamic sampling and retraining of the NN to improve the approximation in regions of interest. We numerically evaluate the algorithm in a case study of optimizing a power system. The numerical results highlight the importance of obtaining an accurate approximation of the expected value function and show that a poor approximation can lead to a suboptimal solution. The results strongly support the proposed sampling and retraining approach, and by the smart selection of data points we are able to obtain a much better approximation in the regions of interest compared to a simple uniform sampling. For the baseline, using uniform sampling, the accuracy does not seem to improve even after doubling the number of data points. Whereas, the proposed algorithm quickly obtains an NN that seems to closely approximate the expected value function, at least in the neighborhood of the minimizer of the NN approximated 2SP (5). The results, thus, support the idea of approximating the expected value function by an NN and using a MILP encoding of the NN to form a single-stage problem, but clearly show the importance of efficiently training a NN to high accuracy in the regions of interest.

As Algorithm 1 delivers promising numerical results on the considered power system, it is natural to ask, whether these results extend to the benchmark problems given in [8], i.e. Capacitated Facility location, Investment, Stochastic server location and pooling problems. Moreover, also extensions to ACOPF
models similarly to the ones in [14] may pose interesting challenges for future research.

#### Acknowledgements

We would like to thank the anonymous referees for their very valuable comments, that helped to significantly improve the quality of this article.
