# A mixed-integer approximation of robust optimization problems with mixed-
integer adjustments

Jan Kronqvist, Boda Li, Jan Rolfes  KTH - Royal Institute of Technology,
Stockholm, Sweden, Department of Mathematics, Lindtstedtsvägen 25, SE-100 44
Stockholm; Sweden<EMAIL_ADDRESS>ABB Corporate Research Center,
Ladenburg, Germany, Wallstadter Str. 59, 68526 Ladenburg
<EMAIL_ADDRESS>
(Date: August 27, 2024)

###### Abstract.

In the present article we propose a mixed-integer approximation of adjustable-
robust optimization (ARO) problems, that have both, continuous and discrete
variables on the lowest level. As these trilevel problems are notoriously hard
to solve, we restrict ourselves to weakly-connected instances. Our approach
allows us to approximate, and in some cases exactly represent, the trilevel
problem as a single-level mixed-integer problem. This allows us to leverage
the computational efficiency of state-of-the-art mixed-integer programming
solvers. We demonstrate the value of this approach by applying it to the
optimization of power systems, particularly to the control of smart
converters.

##  1\. Introduction

Optimization under uncertainty is a very active and growing subfield in
mathematical optimization as in many real world applications, the input data
is either not known in advance or may be subject to perturbations. In the
current literature two major approaches have been developed in order to
address this challenge, stochastic optimization (SO), see e.g. [BL11] for
further details, and robust optimization (RO), see e.g. [BBC11]. In stochastic
optimization it is generally assumed that the underlying probability
distribution of the uncertain parameters is known, whereas robust optimization
methods do not require any knowledge about the distribution but assume that
the uncertain data is contained in a predefined set of scenarios. In addition,
robustly optimal solutions tend to be rather conservative as one considers the
worst-case scenario within the predefined uncertainty set.

In order to address the conservativeness inherent in robust optimization,
adjustable robust optimization (ARO) was introduced in [Ben+04] and grew to a
very active subarea in robust optimization. On the one hand, it serves as a
natural extension of the concept of RO, where uncertainties are not only
addressed a priori by making first stage or here-and-now decisions, but also
allowing to react to the realization of the uncertainty by second stage or
wait-and-see decisions. Hence, ARO is often used to provide more competitive
solutions without sacrificing the uncertainty protection that is provided by
RO. This is particularly useful in applications such as gas networks, see e.g.
[Aßm+18] and the seminal unit commitment problem in power systems, see e.g.
[LS17].

Methodically, AROs are approached either by approximation schemes or by exact
reformulations. In [Ben+04] an approximation scheme was presented, where the
second stage decisions are restricted to affinely depend on the uncertain
parameters. This so-called affine decision rule significantly simplifies the
underlying trilevel problem and often leads to tractable optimization problems
such as linear or semidefinite programs. Additionally, more elaborate decision
rules such as piecewise-linear decision rules, see [BG15], can also address
mixed-integer AROs. For further details, we refer to the excellent survey
[YGH19].

Exact reformulations of AROs are usually computationally intractable,
particularly if non-convexities such as binary variables are involved.
Consequently the literature pivots to decomposition approaches, see [BÖ08],
[ZZ13] and [Ber+13]. Furthermore, in [AP19] a parametric programming approach
is given, that solves AROs with a small number of binary variables to global
optimality. For larger instances however, parametric programming as of now
seems to fall behind decomposition approaches in terms of computational
runtime. Hence, the present work aims to contribute to close this
computational gap. To this end, we introduce a single-level MIP, that
approximates a significant class of AROs and can be proven to be exact under
some further assumptions.

As mentioned above, adjustable robust optimization has been active research
topic in recent years and it has been found useful in a large variety of
areas. Applications dealing with adjustable robustness include, model
predictive control [Tej+19] for dealing with process-model mismatch, process
scheduling under uncertainty [LG16, LI08, LJF04, Gro+16], multi-task
scheduling with imperfect tasks [Lap+19], for analysing resiliency and
flexibility of chemical process [GCG14, ZGL16], to handle uncertainty in
supply chains [Ben+05, BH17], and portfolio optimization [TTT08]. Here, we
will not go into more details on these applications, but will elaborate on the
application of AROs to find the optimal operating points of smart converters
in power systems.

Renewable energy, particularly wind and photovoltaics, has seen rapid growth
in recent years [Bha+15]. These energy sources are widely connected to the
power grid, which can greatly reduce the use of fossil fuels and electricity
costs. However, this integration also introduces power uncertainty that can
threaten the stability and safety of the grid [Wei+14].

Smart inverters [Zha+18] and energy storage technology [Wei+15] provide a
solution to this problem. The former controls the power output of renewables,
while the latter could smooth power fluctuations through charging and
discharging. One of the goals of power system optimization is to balance the
use of renewable energy and system safety by setting the operating state and
set points of inverters and storages.

Historically, power system operation-related optimization involves
deterministic optimization, see e.g. [AC17] for a broader overview on the
seminal unit commitment problem and [Wan+95] for an application to security-
constrained unit commitment or [SM87] for security-constrained economic
dispatch, assuming accurate renewable energy forecasting. However, this stream
of research does not account for potential uncertainties in the system, which
may lead to infeasible operating states. Addressing these challenges has
gained increased attention due to the incorporation of renewables in modern
power systems.

Robust optimization has gained attention for its potential to address this
issue and has been extensively studied in power system operation. We refer to
Conejo et. al. [CW22] for a brief survey. Moreover, multi-level robust
optimization has also be applied to a various challenges that arise in power
system operation, such as robust unit commitment [ZG13a, LS16] or robust
optimal power flow [ZG13] to name a few. It is also applied to the energy
management of power systems with smart inverters and energy storage [Yu+20,
Lek+21, YS21]. While these studies provide excellent use cases of multi-level
optimization models in power systems, they tend to oversimplify the discrete
variables in inverter and storage control models by converting them into
continuous variables for increased tractability. However, this may result in
infeasible solutions in practical applications since we overestimate the
possibilities, the grid operator has at its disposal to readjust the grid.
Conversely, the proposed method in this work preserves discrete variables and
thus guarantees the feasibility of the computed operating points. Moreover, as
our only approximation comes from relaxing the adversarial problem posed by
the uncertainties in the system, we provide (potentially overly conservative)
operating points, whose feasibility can be guaranteed.

This work is structured as follows. In Section 2, the adjustable robust
problem formulation is introduced. Moreover, we discuss the key assumption,
that these problems are weakly connected and present our main results. Section
3 introduces the application to smart converters in networks and demonstrates
the consequences of the approximation results from Section 2 in this setting.
Subsequently, Section 4 illustrates the results and presents numerical
evidence, that the presented method outperforms existing algorithms on weakly-
connected trilevel problems.

##  2\. A MIP approach to robust optimization with MIP adjustments

Adjustable robust programming consist of three separate types or _levels_ of
variables. The first-level variables, denoted in the present paper by
$x\in\mathcal{X}$, describe an initial planning approach, which has to be
decided first. Then, an uncertainty affects the outcome of this initial
planning. Here, we model this uncertainty by a random vector
$h\in\Omega\subseteq\mathbb{R}^{I}$, which is distributed by a probability
distribution $\mathbb{P}\in\mathcal{P}(\Omega)$ on a compact domain
$\Omega\subseteq\mathbb{R}^{I}$. The set $\mathcal{P}(\Omega)$ is called the
_ambiguity set_ of probability measures and the vector $h$ is referred to as
the second-level variable. Lastly, the initial planning $x$ can be adjusted to
the uncertainty $h$ by choosing the third-level variables
$y\in\mathcal{Y}(x,h)$, where in the present article, we suppose that
$\mathcal{Y}(x,h)\subseteq\mathbb{R}^{n}\times\\{0,1\\}^{l}$ is assumed to be
defined through linear constraints. Thus, in total the basic DRO setting can
be defined as follows:

$\min_{x\in\mathcal{X}}G(x)+\max_{\mathbb{P}\in\mathcal{P}(\Omega)}\left(\mathbb{E}_{\mathbb{P}}\left(\min_{y\in\mathcal{Y}(x,h)}c^{\top}y\right)\right),$
(1)

where $\mathcal{X}$ and $\Omega$ are assumed to be a polytopes and
$\mathcal{Y}(x,h)$ a polytope intersected with an integer lattice. However, in
the present article, we restrict ourselves to a standard robust ambiguity set,
i.e., for the polytope $\Omega$ the ambiguity set is defined as a set of Dirac
measures $\mathcal{P}(\Omega)\coloneqq\\{\delta_{\\{h\\}}:h\in\Omega\\}$ and
aim to solve instances of

$\min_{x\in\mathcal{X}}G(x)+\max_{h\in\Omega}\min_{y\in\mathcal{Y}(x,h)}c^{\top}y.$
(2)

Note, that (2) is still considered to be a very challenging problem as it
contains the $\mathrm{NP}$-complete MIP $\min_{y\in\mathcal{Y}(x,h)}c^{\top}y$
as a subproblem, see Section 6 in [YGH19] for further details. In the present
article, we focus on turning the above trilevel-problem into a single-level
MIP. To this end, we fix the integral part of $y$ and dualize the resulting LP
in order to achieve a bilevel problem with a bilinear objective. This
objective is then relaxed with McCormick envelopes and dualized again in order
to achieve a single-level LP for every fixed integral assignment of $y$.
However, instead of evaluating the resulting LP relaxations for every integer
combination, we can simply reincorporate the integrality of $y$ as an
additional constraint and, thus, obtain a single-level MIP. This approach may
lead to overly pessimistic outcomes as using the McCormick relaxations
strengthens the adverserial in the multilevel problem, but can be shown to be
exact under some assumptions on $h$ and the dual variables of
$\min_{y\in\mathcal{Y}(x,h)}c^{\top}y$.

However, this method may not work on strongly connected multilevel problems
since broadly speaking the hardness in multilevel optimization problems comes
from strong connections between the levels, we consider the following subclass
of adjustable robust problems:

We define a problem of type (2) _weakly connected_ , if the only relation
between the first-,second- and third level variables is a linear relation in
$\mathcal{Y}(x,h)$:

$\mathcal{Y}(x,h)=\\{y\in\mathbb{R}^{n}\times\\{0,1\\}^{l}:A^{\prime}y\geq
b^{\prime},By\geq B_{x}x+B_{h}h+b_{0}\\},$

i.e., the first-level variables $x$ and the second-level random vector only
affect the right-hand side of the constraints defining $\mathcal{Y}(x,h)$.
Weakly connected adjustable robust problems include among others, problems
with affine decision rules
($B^{\top}=(I_{n},-I_{n}),B_{x}=0,B_{h}^{\top}=(Q^{\top},-Q^{\top}),b_{0}=(q^{\top},-q^{\top})$,
where $Q,q$ can be chosen arbitrarily.

Moreover, in order to work with the levels separately, we denoted by
$A^{\prime}y\geq b^{\prime}$ the constraints, that solely deal with third-
level variables, i.e., that are neither affected by the first-level variables
$x$ nor by the second-level variables $h$ and by $By\geq B_{x}x+B_{h}h+b_{0}$
the constraints, that are affected by the upper levels. Furthermore, instances
that are particularly well suited for our approach should have relatively few
non-zero entries in the matrix $B_{h}$. Few non-zero entries in $B_{h}$ is not
a strict requirement, but both the approximation quality and computational
complexity can scale with the number of non-zero elements.

In addition, let us consider an LP inner-approximation of the third-level
program $\min_{y\in\mathcal{Y}(x,h)}c^{\top}y$, i.e. we fix the integer
variables $y_{n+1},\ldots,y_{n+l}\in\\{0,1\\}$ by a set of linear constraints
denoted by $A_{f}y\geq b_{f}$. Consequently, we consider the following LP:

$\displaystyle\min\ $ $\displaystyle c^{\top}y$ (3a) s.t. $\displaystyle
A^{\prime}y\geq b^{\prime},$ (3b) $\displaystyle A_{f}y\geq b_{f}$ (3c)
$\displaystyle By\geq B_{x}x+B_{h}h+b_{0}$ (3d)

Let
$A^{\top}=((A^{\prime})^{\top},A_{f}^{\top}),b^{\top}=((b^{\prime})^{\top},b_{f}^{\top})$
and $\alpha$ denote the dual variables that correspond to Constraints (3b) and
(3c), i.e. to $Ay\geq b$. Additionally, we denote by $\beta$ the dual
variables corresponding to (3d) and obtain as the dual program of (3):

$\displaystyle\max\ $ $\displaystyle
b^{\top}\alpha+(B_{x}x+B_{h}h+b_{0})^{\top}\beta$ (4a) s.t.
$\displaystyle\begin{pmatrix}A^{\top}&B^{\top}\end{pmatrix}\begin{pmatrix}\alpha\\\
\beta\end{pmatrix}=c,$ (4b) $\displaystyle\alpha,\beta\geq 0,$ (4c)

We observe, that on the second level the objective decomposes into a bilinear
part $(B_{h}h)^{\top}\beta=h^{\top}B_{h}\beta$ and a linear one
$b^{\top}\alpha+(B_{x}x+b_{0})^{\top}\beta$. Moreover, the amount of nonzero
entries of $B_{h}$ crucially determines the number of bilinear terms and
thereby the difficulty of computing (4).

Let $I$ and $J$ denote the index sets of $h$ and $\beta$ respectively. Then,
in order to relax the bilinear term $h^{\top}B_{h}\beta$, we make the
following key assumption:

######  Assumption 1.

Both, the dual variables $\beta\in\mathbb{R}^{J}$ as well as the second-level
variables $h\in\mathbb{R}^{I}$ are bounded, i.e., there are
$\beta^{-},\beta^{+}$ and $h^{-},h^{+}$ such that adding

$h^{-}\leq h\leq h^{+}\text{ and }\beta^{-}\leq\beta\leq\beta^{+}$

to $\Omega$ or (4) does not affect the outcome of (1).

Note that a compact $\Omega$ already implies $h^{-}\leq h\leq h^{+}$. Hence,
in those cases it suffices to check Assumption 1 for $\beta$. Later, we show
that finite, and meaningful, bounds $\beta^{-}$ and $\beta^{+}$ can be easily
determined for the considered application of finding optimal operating points
in power systems. In addition, we mention, that every sharpening of the bounds
in Assumption 1 improves the quality of the upcoming results. This is due to
the fact, that the boundedness of both, $\beta$ and $h$ enables us to relax
the second-level problem
$\max_{h\in\Omega}\min_{y\in\mathcal{Y}(x,h)}c^{\top}y$ further using
McCormick envelopes for the bilinear terms. In particular, this idea gives
rise to the following theorem:

######  Theorem 1.

Let the ambiguity set $\Omega$ be a polytope defined by
$\Omega=\\{h\in\mathbb{R}^{I}:A_{\Omega}^{\top}h+B_{\Omega}^{\top}\eta=b_{\Omega},\
\eta\geq 0\\}$ with $b_{\Omega}\in\mathbb{R}^{k}$, i.e. $\Omega$ is compact
and $\eta\geq 0$ denote potential nonnegative slack variables in the rows. Let
further $\beta^{-},\beta^{+}$ be a lower/upper bound for $\beta$. Then, the
following linear program provides an upper bound to (2):

$\displaystyle\min\ $ $\displaystyle
G(x)+(\beta^{+})^{\top}u_{\beta}^{+}+(\beta^{-})^{\top}u_{\beta}^{-}+b_{\Omega}^{\top}u_{\Omega}+c^{\top}y$
(5a) $\displaystyle-\sum_{i\in I,j\in
J}\left(h_{i}^{-}\beta_{j}^{-}(u_{\text{env}})_{ij}^{1}+h_{i}^{+}\beta_{j}^{+}(u_{\text{env}})_{ij}^{2}+h_{i}^{+}\beta_{j}^{-}(u_{\text{env}})_{ij}^{3}+h_{i}^{-}\beta_{j}^{+}(u_{\text{env}})_{ij}^{4}\right)$
s.t. $\displaystyle Ay\geq b,$ (5b) $\displaystyle
By+u_{\beta}^{+}+u_{\beta}^{-}-h_{i}^{-}(u_{\text{env}})_{i}^{1}-h_{i}^{+}(u_{\text{env}})_{i}^{2}-h_{i}^{+}(u_{\text{env}})_{i}^{3}-h_{i}^{-}(u_{\text{env}})_{i}^{4}\geq
B_{x}x+b_{0}$ $\displaystyle\text{for every }i\in I$ (5c) $\displaystyle
u_{\beta}^{+}\geq 0,\ -u_{\beta}^{-}\geq 0,$ (5d)
$\displaystyle(A_{\Omega}u_{\Omega})_{i}-\sum_{j\in
J}\left(\beta_{j}^{-}(u_{\text{env}})_{ij}^{1}+\beta_{j}^{+}(u_{\text{env}})_{ij}^{2}+\beta_{j}^{-}(u_{\text{env}})_{ij}^{3}+\beta_{j}^{+}(u_{\text{env}})_{ij}^{4}\right)\geq
0$ $\displaystyle\text{for every }i\in I,$ (5e) $\displaystyle
B_{\Omega}u_{\Omega}\geq 0,$ (5f)
$\displaystyle(u_{\text{env}})_{ij}^{1}+(u_{\text{env}})_{ij}^{2}+(u_{\text{env}})_{ij}^{3}+(u_{\text{env}})_{ij}^{4}\geq(B_{h})_{ij}$
$\displaystyle\text{for every }i\in I,j\in J,$ (5g)
$\displaystyle-(u_{\text{env}})_{ij}^{1},-(u_{\text{env}})_{ij}^{2}\geq 0$
$\displaystyle\text{for every }i\in I,j\in J,$ (5h)
$\displaystyle(u_{\text{env}})_{ij}^{3},(u_{\text{env}})_{ij}^{4}\geq 0$
$\displaystyle\text{for every }i\in I,j\in J,$ (5i) $\displaystyle
x\in\mathcal{X}.$ (5j)

###### Proof.

We observe that with Assumption 1 the second-level
$\max_{h\in\Omega}\min\\{c^{\top}y:Ay\geq b,By\geq B_{x}x+B_{h}h+b_{0}\\}$ can
be written as

$\displaystyle\max\ $ $\displaystyle
b^{\top}\alpha+(B_{x}x+b_{0})^{\top}\beta+\sum_{i\in I,j\in
J}(B_{h})_{ij}h_{i}\beta_{j}$ (6a) s.t. $\displaystyle
A^{\top}\alpha+B^{\top}\beta=c,$ (6b) $\displaystyle\beta+\gamma=\beta^{+},$
(6c) $\displaystyle\beta-\delta=\beta^{-},$ (6d) $\displaystyle
A_{\Omega}^{\top}h+B_{\Omega}^{\top}\eta=b_{\Omega},$ (6e)
$\displaystyle\alpha,\beta,\gamma,\delta,h,\eta\geq 0,$ (6f)

where we assumed w.l.o.g. that $h^{-}=0$ since otherwise, we could substitute
$h$ by $h-h^{-}$ and adjust $b_{0}$ and $b_{\Omega}$ accordingly. Next, we
substitute $\kappa_{ij}\coloneqq h_{i}\beta_{j}$ in the objective and relax
the resulting constraint $\kappa_{ij}=h_{i}\beta_{j}$ by a McCormick envelope.
Note, that w.l.o.g. $\kappa\geq 0$. If we further introduce suitable
nonnegative slack variables $\rho\geq 0$, we obtain the following LP:

$\displaystyle\max\ $ $\displaystyle
b^{\top}\alpha+(B_{x}x+b_{0})^{\top}\beta+\sum_{i\in I,j\in
J}(B_{h})_{ij}\kappa_{ij}$ (7a) s.t. $\displaystyle
A^{\top}\alpha+B^{\top}\beta=c,$ (7b) $\displaystyle\beta+\gamma=\beta^{+},$
(7c) $\displaystyle\beta-\delta=\beta^{-},$ (7d) $\displaystyle
A_{\Omega}^{\top}h+B_{\Omega}^{\top}\eta=b_{\Omega},$ (7e)
$\displaystyle\kappa_{ij}=h_{i}^{-}\beta_{j}+h_{i}\beta_{j}^{-}-h_{i}^{-}\beta_{j}^{-}+\rho_{ij}^{1}$
$\displaystyle\text{for every }i\in I,j\in J$ (7f)
$\displaystyle\kappa_{ij}=h_{i}^{+}\beta_{j}+h_{i}\beta_{j}^{+}-h_{i}^{+}\beta_{j}^{+}+\rho_{ij}^{2}$
$\displaystyle\text{for every }i\in I,j\in J,$ (7g)
$\displaystyle\kappa_{ij}=h_{i}^{+}\beta_{j}+h_{i}\beta_{j}^{-}-h_{i}^{+}\beta_{j}^{-}-\rho_{ij}^{3}$
$\displaystyle\text{for every }i\in I,j\in J,$ (7h)
$\displaystyle\kappa_{ij}=h_{i}\beta_{j}^{+}+h_{i}^{-}\beta_{j}-h_{i}^{-}\beta_{j}^{+}-\rho_{ij}^{4}$
$\displaystyle\text{for every }i\in I,j\in J,$ (7i)
$\displaystyle\alpha,\beta,\gamma,\delta,h,\eta,\kappa,\rho\geq 0.$ (7j)

If we denote the dual variables of (7) by
$y,u_{\beta}^{+},u_{\beta}^{-},u_{\Omega},(u_{\text{env}})_{ij}^{1},(u_{\text{env}})_{ij}^{2},(u_{\text{env}})_{ij}^{3},(u_{\text{env}})_{ij}^{4}$
respectively, then the result follows by strong duality and including the
first-level variables and objectives. ∎

We observe, that Theorem 1 provides an LP inner approximation of (2) in the
sense that the adversarial is overestimated, whereas the space of decision
variables is underestimated, potentially leading to a more conservative
solution. In particular, we simply fixed the discrete or w.l.o.g. binary
decisions in $\mathcal{Y}(x,h)$ denoted by $y_{m+1},\ldots
y_{m+l}\in\\{0,1\\}$. To this end, we briefly illustrate the key idea behind
the following Theorem 2, that addresses integralities:

Observe, that solving (2) is equivalent to solving exponentially many trilevel
LPs – one LP for every fixed choice of $y_{m+1},\ldots y_{m+l}$. This is due
to the fact that the proof of Theorem 1 does not depend on the fixing of
$(y_{m+1},\ldots,y_{m+l})^{\top}$ to a vector in $y^{\prime}\in\\{0,1\\}^{l}$.
Consequently, we obtain $2^{l}$ single-level LPs, where each of these LPs
yields an inner-approximation of the trilevel LPs given by the fixing of
$y_{m+1},\ldots,y_{m+l}$. Finally, instead of solving each of those
exponentially many LPs separately, the following MIP gives the same result.

######  Theorem 2.

Let the ambiguity set $\Omega$ be a polytope defined by
$\Omega=\\{h\in\mathbb{R}^{I}:A_{\Omega}^{\top}h+B_{\Omega}^{\top}\eta=b_{\Omega},\
\eta\geq 0\\}$ with $b_{\Omega}\in\mathbb{R}^{k}$, i.e.
$\Omega\subseteq\mathbb{R}^{I}$ is compact and $\eta\geq 0$ denote potential
nonnegative slack variables in the rows. Let further $\beta^{-},\beta^{+}$ be
a lower/upper bound for $\beta$. Then, the following MIP provides an upper
bound to the tri-level MIP (2):

$\displaystyle\min\ $ $\displaystyle
G(x)+(\beta^{+})^{\top}u_{\beta}^{+}+(\beta^{-})^{\top}u_{\beta}^{-}+b_{\Omega}^{\top}u_{\Omega}+c^{\top}y$
(8a)
$\displaystyle-\sum_{i,j\in[m]}\left(h_{i}^{-}\beta_{j}^{-}(u_{\text{env}})_{ij}^{1}+h_{i}^{+}\beta_{j}^{+}(u_{\text{env}})_{ij}^{2}+h_{i}^{+}\beta_{j}^{-}(u_{\text{env}})_{ij}^{3}+h_{i}^{-}\beta_{j}^{+}(u_{\text{env}})_{ij}^{4}\right)$
s.t. $\displaystyle Ay\geq b,$ (8b) $\displaystyle
By+u_{\beta}^{+}+u_{\beta}^{-}-h_{i}^{-}(u_{\text{env}})_{i}^{1}-h_{i}^{+}(u_{\text{env}})_{i}^{2}-h_{i}^{+}(u_{\text{env}})_{i}^{3}-h_{i}^{-}(u_{\text{env}})_{i}^{4}\geq
B_{x}x+b_{0}$ $\displaystyle\text{for every }i\in[m]$ (8c) $\displaystyle
u_{\beta}^{+}\geq 0,\ -u_{\beta}^{-}\geq 0,$ (8d)
$\displaystyle(A_{\Omega}u_{\Omega})_{i}-\sum_{j\in[m]}\left(\beta_{j}^{-}(u_{\text{env}})_{ij}^{1}+\beta_{j}^{+}(u_{\text{env}})_{ij}^{2}+\beta_{j}^{-}(u_{\text{env}})_{ij}^{3}+\beta_{j}^{+}(u_{\text{env}})_{ij}^{4}\right)\geq
0$ $\displaystyle\text{for every }i\in[m],$ (8e) $\displaystyle
B_{\Omega}u_{\Omega}\geq 0,$ (8f)
$\displaystyle(u_{\text{env}})_{ij}^{1}+(u_{\text{env}})_{ij}^{2}+(u_{\text{env}})_{ij}^{3}+(u_{\text{env}})_{ij}^{4}\geq(B_{h})_{ij}$
$\displaystyle\text{for every }i,j\in[m],$ (8g)
$\displaystyle-(u_{\text{env}})_{ij}^{1},-(u_{\text{env}})_{ij}^{2}\geq 0$
$\displaystyle\text{for every }i,j\in[m],$ (8h)
$\displaystyle(u_{\text{env}})_{ij}^{3},(u_{\text{env}})_{ij}^{4}\geq 0$
$\displaystyle\text{for every }i,j\in[m],$ (8i) $\displaystyle
x\in\mathcal{X},u_{\beta}^{+},u_{\beta}^{-}\in\mathbb{R}^{J},u_{\Omega}\in\mathbb{R}^{k},(u_{\text{env}})^{1},(u_{\text{env}})^{2},(u_{\text{env}})^{3},(u_{\text{env}})^{4}\in\mathbb{R}^{I\times
J},y\in\mathbb{R}^{n+l},$ (8j) $\displaystyle
y_{m+1},\ldots,y_{m+l}\in\\{0,1\\}.$ (8k)

###### Proof.

First, we fix $(y_{m+1},\ldots,y_{m+l})^{\top}=y^{\prime}$ with a set of
linear (in-)equalities $A_{f}y\geq b_{f}$. Here, $y^{\prime}\in\\{0,1\\}^{l}$
denotes an arbitrary integer assignment. Suppose we replace (8k) by this
system of inequalities and denote the resulting feasible set by
$\mathcal{X}(y^{\prime})$, i.e.

$\mathcal{X}(y^{\prime})=\left\\{x,u_{\beta}^{+},u_{\beta}^{-},u_{\Omega},(u_{\text{env}})^{1},(u_{\text{env}})^{2},(u_{\text{env}})^{3},(u_{\text{env}})^{4},y:\eqref{Constr:
single_level_MIP_general_C1}-\eqref{Constr:
single_level_MIP_general_env2_lb},(y_{m+1},\ldots,y_{m+l})^{\top}=y^{\prime}\right\\}.$

We further denote the corresponding objective formed by (8a) by $c(x)$, i.e.
the LP $\min_{x\in\mathcal{X}(y^{\prime})}c(x)$ describes (8) with fixed
values of $(y_{m+1},\ldots,y_{m+l})^{\top}$. Hence, Theorem 1 implies, that
$\min_{x\in\mathcal{X}(y^{\prime})}c(x)$ inner-approximates the following
trilevel LP

$\min_{x\in\mathcal{X}}G(x)+\max_{h\in\Omega}\min_{y\in\mathcal{Y}(x,h,y^{\prime})}c^{\top}y,$

where

$\mathcal{Y}(x,h,y^{\prime})=\left\\{y\in\mathcal{Y}(x,h):(y_{m+1},\ldots,y_{m+l})^{\top}=y^{\prime}\right\\}.$

In particular, we obtain by Theorem 1 that
$\min_{x\in\mathcal{X}(y^{\prime})}c(x)\geq\min_{x\in\mathcal{X}}G(x)+\max_{h\in\Omega}\min_{y\in\mathcal{Y}(x,h,y^{\prime})}c^{\top}y$
for every $y^{\prime}\in\\{0,1\\}^{l}$. Subsequently, we conclude

$\displaystyle\eqref{Prob: single_level_MIP_general}$
$\displaystyle=\min_{y^{\prime}\in\\{0,1\\}^{l}}\min_{x\in\mathcal{X}(y^{\prime})}c(x)$
$\displaystyle\geq\min_{y^{\prime}\in\\{0,1\\}^{l}}\min_{x\in\mathcal{X}}G(x)+\max_{h\in\Omega}\min_{y\in\mathcal{Y}(x,h,y^{\prime})}c^{\top}y$
$\displaystyle=\min_{x\in\mathcal{X}}G(x)+\min_{y^{\prime}\in\\{0,1\\}^{l}}\max_{h\in\Omega}\min_{y\in\mathcal{Y}(x,h,y^{\prime})}c^{\top}y$
$\displaystyle\overset{(*)}{\geq}\min_{x\in\mathcal{X}}G(x)+\max_{h\in\Omega}\min_{y^{\prime}\in\\{0,1\\}^{l},y\in\mathcal{Y}(x,h,y^{\prime})}c^{\top}y=\eqref{Prob:
entire_adjustable_DRO}.$

As stated above, the first inequality is based on Theorem 1 whereas the second
one is a consequence of the max-min-inequality. In the remainder of this
proof, we will argue, that $(*)$ is even sharp and thus, any potential
differences between the relaxation (8) and the original problem (1) solely
originate from the McCormick envelopes. To this end, let us consider the
Lagrangian relaxation with penalty terms instead of hard constraints:

$f_{x}(h,y^{\prime})\coloneqq\min_{y\in\mathcal{Y}(x)}c^{\top}y+u_{0}^{\top}(A_{\Omega}^{\top}h+B_{\Omega}^{\top}\eta-
b_{\Omega})+u_{1}^{\top}((y_{m+1},\ldots,y_{m+l})^{\top}-y^{\prime}),$

where $\mathcal{Y}(x)\supseteq\mathcal{Y}(x,h)$ is defined as the relaxation
of $\mathcal{Y}(x,h)$ occuring if Constraint
$A_{\Omega}^{\top}h+B_{\Omega}^{\top}\eta=b_{\Omega}$ is dropped. Observe that
since $\mathcal{Y}(x)$ is bounded, for a given $h\in\Omega$ and
$y^{\prime}\in\\{0,1\\}^{l}$, according to Theorem 21 in [Lem01], there exist
$u_{0}\in\mathbb{R}^{k},u_{1}\in\mathbb{R}^{l}$ that satisfy

$\min_{y\in\mathcal{Y}(x,h,y^{\prime})}c^{\top}y=f_{x}(h,y^{\prime}).$

Then, since $f_{x}(h,y^{\prime})$ is concave in $h$, i.e. on the nonempty set
$\Omega$ and further is convex in $y^{\prime}$ on the compact nonempty set
$\\{0,1\\}^{l}$, we apply the Ky-Fan theorem [Fan53] and conclude equality in
$(*)$. ∎

We note, that the approximation quality given by Theorem 2 relies solely on
the strenght of the relaxation due to the McCormick envelopes. Moreover, as
the number of nonzero entries in $B_{h}$ determines the number of bilinear
terms in the second-level, it crucially affects the approximation quality
given by the McCormick relaxation. In particular, a large number of bilinear
terms may lead to an overly conservative solution. Thus, exploiting problem
specific information that improves the bounds
$\beta^{-},\beta^{+},h^{-},h^{+}$ directly increases the solution quality.
Moreover, more elaborate approximations of $h^{\top}B_{h}\beta$ may lead to a
direct improvement of Theorem 2 and are subject to future research.

##  3\. Application to smart converters in power systems networks

In the present section, we apply our results to the questions of how to
operate a power grid integrated with smart inverters and storage.

Usually, the optimal operation problem of the power system can be divided into
two phases.

  * •

The first phase is about day-ahead scheduling. In this phase, the power demand
of the power grid needs to be determined and reported to the day-ahead
electricity market. In addition, internal conventional generators need to make
a power generation plan for the next day to reasonably allocate diesel
consumption.

  * •

The second phase is about intra-day system operation. Given the day-ahead
decisions, operators need to further set operating points and working status
of smart inverters and energy storage systems to realize a real-time power
balance and reduce the waste of renewables.

As can be seen from our analysis, the day-ahead decisions influence the intra-
day decisions. The renewable-led uncertainties will further affect the quality
of decision-making. In order to realize the safe and economical operation of
the power grid, we need to overcome the following three difficulties:

  1. (a)

Sequential: The formulation of the day-ahead strategy is subject to
uncertainty realizations and intra-day operations. A multi-layer model is
needed to describe the decision sequence as well as uncertainty realizations
in the real world.

  2. (b)

Uncertain: The uncertainties threaten the safe operation of power systems.
Power fluctuations induced by renewables may deteriorate power quality and
increase electricity costs. Hence, the operation strategy should be robust
enough to handle different situations.

  3. (c)

Discrete: The intra-day operation in the second phase often involves many
state-switching operations. The states can be modeled as integer variables in
the optimization problem. However, these variables destroy the convexity of
the model. A new solution method should be developed in such a case.

The power system operation problem can be solved using the robust optimization
approach with MIP adjustments proposed in Section 2. A multi-layer model is
built to describe the sequential decision-making process; the adaptive robust
programming is employed to account for worst-case scenarios; the MIP
adjustment is then applied to consider the state switching in the system.

At first, we formulate the mathematical model of the power system operation.
We mainly follow the notation by [YW19] but use the simpler DC approximation
of Kirchhoff laws in order to model the power flow when operating the grid.

Let $\mathcal{B}$ denote a set of buses and $\mathcal{L}$ denote a set of
lines/branches in a power grid. Additionally, power is generated within the
grid either 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}$. Moreover, the transmission system operator (TSO) may
decide to store or release power via a set $\mathcal{N}_{S}$ of storages, e.g.
batteries. The last potential sources and sinks of power is a trading node
with other connected regional transmission grids. Here, the TSO may purchase
power on the day-ahead market or intra-day. As the day-ahead market is often
called first-level market, we denote the amount of power traded day-ahead by
$P_{fl}$ and the corresponding market price by $p_{fl}$. Similarly, the amount
of power traded intra-day is denoted by $P_{sl}$ and its market price by
$p_{sl}$. Note, that a positive value for $P_{fl},P_{sl}$ is interpreted as a
purchase and negative values for $P_{fl},P_{sl}$ correspond to a sell of
energy.

For the SO’s initial planning, one considers the first-level variables
$x=(P_{G}^{\top},P_{fl})^{\top}\in\mathbb{R}^{\mathcal{N}_{G}}\times\mathbb{R}$
combined with the estimated renewable energy production
$P_{DG,forecast}\in\mathbb{R}^{\mathcal{N}_{DG}}$ and ensures that a given
total demand $\sum_{i\in\mathcal{B}}P_{d_{i}}$ is met. Hence,

$\mathcal{X}=\\{x\in\mathbb{R}^{n}:\ \eqref{Constr: fl1_define_Pfl}\ \&\
\eqref{Constr: fl2_bounds_P_G}\\},$

where

$\displaystyle
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}^{t}$
$\displaystyle\text{for every }t\in T,$ (9a) $\displaystyle
P_{G_{i},\text{min}}^{t}\leq P_{G_{i}}^{t}\leq P_{G_{i},\text{max}}^{t}$
$\displaystyle\text{for every }i\in\mathcal{N}_{G},t\in T.$ (9b)

with given parameters
$P_{G,\min},P_{G,\max}\in\mathbb{R}^{\mathcal{N}_{G}\times T}$ and the
subsequent constraints (9a) and (9b). We like emphasize that the _market-
clearing condition_ (9a) ensures the active power balance in the whole system.
Moreover, the objective of the first-level is given by

$G(x)=\sum_{i\in\mathcal{N}_{G},t\in
T}c_{G_{i}^{t},2}(P_{G_{i}}^{t})^{2}+c_{G_{i}^{t},1}P_{G_{i}}^{t}+c_{G_{i}^{t},0}+\sum_{t\in
T}p_{fl}^{t}P_{fl}^{t},$

where $c_{G_{i}^{t},2},c_{G_{i}^{t},1},c_{G_{i}^{t},0}\in\mathbb{R}$ are given
generator cost parameters.

Since the uncertainties will impact the initial planning, we will consider the
uncertain capacity of the renewable generators, i.e. we denote the second-
level variable by $h=P_{DG,\text{max}}$. 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 high proportion of renewable
energies [PHK14, INO20]. As the grid stability is crucial, we will address
this uncertainty in a robust manner and set the domain of the second-level
variables $P_{DG,\max}$ to

$\Omega=\left\\{P_{DG,\text{max}}\in\mathbb{R}^{\mathcal{N}_{DG}\times T}:\
0\leq P_{DG_{i},\text{max}}^{t}\leq P_{i}^{+}\forall\
i\in\mathcal{N}_{DG},t\in
T,\sum_{i\in\mathcal{N}_{DG}}P_{DG_{i},\text{max}}^{t}\geq
R\sum_{i\in\mathcal{N}_{DG}}P_{DG_{i},forecast}^{t}\right\\},$ (10)

where $P_{i}^{+}$ denotes the technical limit of the renewable generator, i.e.
its capacity under optimal conditions and $R\in[0,1]$ denotes a maximal
forecast error.

On the third level, the TSO is able to react to this uncertainty and adjust
the initial planning accordingly. In particular, instead of producing $P_{G}$,
the TSO might regulate the energy output to $P_{G,\text{reg}}$ by either
increasing the production by adding $P_{G}^{+}\geq 0$ or decreasing the
production by adding $P_{G}^{-}\leq 0$. However, this can only be done at a
cost $r^{+}$ or $r^{-}$ respectively. Similarly, $P_{DG}$ is the adjusted
energy production that deviates from its forecast by $P_{DG}^{+}$ or
$P_{DG}^{-}$, where deviations are penalized by $f^{+},f^{-}$ respectively.
Despite of these regulations, the TSO might trade power intra-day ($P_{sl}$)
or decide whether ($\mu_{ch},\mu_{dch}\in\\{0,1\\}^{\mathcal{N}_{S}}$) and by
how much $P_{ch},P_{dch}$ to charge or discharge the storages for balancing
potential power deficiency and surplus. The state of charge of a storage is
denoted by soc. Lastly, the power on a line $(k,l)\in\mathcal{L}$ is denoted
by $p_{kl}$ and the phase angles of the system by $\vartheta$. Hence, the TSO
can adjust 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}:\ \eqref{Constr:
adjust_fossil_fuel_generators}-\eqref{Constr: mu_bound}\right\\},$

where the constraints (11a)-(15g) are given below:

  1. (a)

First, we consider the following _Generator and DG output constraints_ :

$\displaystyle
P_{G_{i},\text{reg}}^{t}=P_{G_{i}}^{t}+P_{G_{i}}^{t,+}+P_{G_{i}}^{t,-}$
$\displaystyle\text{for every }i\in\mathcal{N}_{G},t\in T,$ (11a)
$\displaystyle P_{G_{i},\text{min}}^{t}\leq P_{G_{i},\text{reg}}^{t}\leq
P_{G_{i},\text{max}}^{t}$ $\displaystyle\text{for every
}i\in\mathcal{N}_{G},t\in T,$ (11b) $\displaystyle
P_{DG_{i}}^{t}=P_{DG_{i},\text{forecast}}^{t}+P_{DG_{i}}^{t,+}+P_{DG_{i}}^{t,-}$
$\displaystyle\text{for every }i\in\mathcal{N}_{DG},t\in T,$ (11c)
$\displaystyle P_{DG_{i},\text{min}}^{t}\leq P_{DG_{i}}^{t}\leq
P_{DG_{i},\text{max}}^{t}$ $\displaystyle\text{for every
}i\in\mathcal{N}_{DG},t\in T,$ (11d) $\displaystyle
P_{sl}^{t}-P_{fl}^{t}=-\mathbbm{1}^{\top}(P_{G}^{t,+}+P_{G}^{t,-})-\mathbbm{1}^{\top}(P_{DG}^{t,+}+P_{DG}^{t,-})-\mathbbm{1}^{\top}(P_{dch}^{t}-P_{ch}^{t})$
$\displaystyle\text{for every }t\in T.$ (11e)

Constraints (11a)-(11d) describe the output range of the conventional and
renewable generators and the potential impact of uncertainties. In particular,
(11d) shows that the third-level variables $P_{DG_{i}}^{t}$ are restricted by
the uncertainties realized in the second level ($P_{DG_{i},\text{max}}^{t}$).
Constraint (11e) illustrates the market clearing on the intra-day market,
which reflect the actual power demand-supply relations.

  2. (b)

Second, we consider the _operation constraints_ :

$\displaystyle 0\leq P_{G_{i}}^{t,+}\leq P_{G_{i},\text{max}}^{t,+}$
$\displaystyle\text{for every }i\in\mathcal{N}_{G},t\in T,$ (12a)
$\displaystyle P_{G_{i},\text{min}}^{t,-}\leq P_{G_{i}}^{t,-}\leq 0$
$\displaystyle\text{for every }i\in\mathcal{N}_{G},t\in T,$ (12b)
$\displaystyle P_{DG_{i}}^{t,+}\geq 0$ $\displaystyle\text{for every
}i\in\mathcal{N}_{DG},t\in T,$ (12c) $\displaystyle P_{DG_{i}}^{t,-}\leq 0$
$\displaystyle\text{for every }i\in\mathcal{N}_{DG},t\in T.$ (12d)

Constraints (12a) -(12d) limit the real output derivations of conventional and
renewable generators.

  3. (c)

Third, we consider the _power flow constraints_ with a DC approximation for
the given, constant line reactance ($x_{ij}>0$) and demand in active power
($P_{k}^{d,t}$) at every time step $t$ and bus $k$. The nodal power flow
balance is established in (13a) and separately in (13b) for the root node. The
branch power flow is established in (13c). The constraints are listed as
follows:

$\displaystyle\sum_{i\in\mathcal{N}_{G}:i\sim
k}P_{G_{i},\text{reg}}^{t}+\sum_{i\in\mathcal{N}_{DG}:i\sim
k}P_{DG_{i}}^{t}+\sum_{i\in\mathcal{N}_{S}:i\sim
k}(P_{dch_{i}}^{t}-P_{ch_{i}}^{t})-P_{k}^{d,t}=\sum_{l\in\delta(k)}p_{kl}^{t}$
$\displaystyle\forall k\in\mathcal{B}\setminus\\{0\\},t\in T,$ (13a)
$\displaystyle\sum_{i\in\mathcal{N}_{G}:i\sim
0}P_{G_{i},\text{reg}}^{t}+\sum_{i\in\mathcal{N}_{DG}:i\sim
0}P_{DG_{i}}^{t}+\sum_{i\in\mathcal{N}_{S}:i\sim
0}(P_{dch_{i}}^{t}-P_{ch_{i}}^{t})+P_{sl}^{t}-P_{0}^{d,t}=\sum_{l\in\delta(k)}p_{0l}^{t}$
$\displaystyle\forall t\in T,$ (13b) $\displaystyle
p_{ij}^{t}=\frac{1}{x_{ij}}(\theta_{i}^{t}-\theta_{j}^{t})$
$\displaystyle\forall\\{i,j\\}\in\mathcal{L},t\in T$ (13c)

  4. (d)

Fourth, to guarantee the safe operation of the branch, the power flow should
not exceed the branch’s capacities. Hence, we consider the _branch thermal
constraints_ :

$\displaystyle-s_{ij,\max}\leq p_{ij}^{t}\leq s_{ij,\max}$
$\displaystyle\text{for every }\\{i,j\\}\in\mathcal{L},t\in T.$ (14a)

  5. (e)

Fifth, we consider the _storage constraints_. Noticing that the storage
operation involves two actions, the action shifts need to be considered. Two
binary variables $\mu_{ch},\mu_{dch}$ are defined to represent the storage
action, where $\mu_{ch},\mu_{dch}\in\\{0,1\\}^{\mathcal{N}_{S}\times T}$.
Moreover, only one of those variables can be 1 at time $t$. Thus, we have the
following constraints:

$\displaystyle\text{soc}_{i,\text{min}}^{t}\leq\text{soc}_{i}^{t}\leq\text{soc}_{i,\text{max}}^{t}$
$\displaystyle\text{for every }i\in\mathcal{N}_{S},t\in T,$ (15a)
$\displaystyle\text{soc}_{i}^{t}=\text{soc}_{i}^{t-1}+\frac{(P_{ch_{i}}^{t}-P_{dch_{i}}^{t})}{E_{i}}\Delta
T$ $\displaystyle\text{for every }i\in\mathcal{N}_{S},t\in T,$ (15b)
$\displaystyle P_{ch_{i}}^{t},P_{dch_{i}}^{t}\geq 0$ $\displaystyle\text{for
every }i\in\mathcal{N}_{S},t\in T,$ (15c)
$\displaystyle\mu_{ch_{i}}^{t},\mu_{dch_{i}}^{t}\in\\{0,1\\}$
$\displaystyle\text{for every }i\in\mathcal{N}_{S},t\in T,$ (15d)
$\displaystyle\mu_{ch_{i}}^{t}P_{ch_{i},\text{min}}^{t}\leq
P_{ch_{i}}^{t}\leq\mu_{ch_{i}}^{t}P_{ch_{i},\text{max}}^{t}$
$\displaystyle\text{for every }i\in\mathcal{N}_{S},t\in T,$ (15e)
$\displaystyle\mu_{dch_{i}}^{t}P_{dch_{i},\text{min}}^{t}\leq
P_{dch_{i}}^{t}\leq\mu_{dch_{i}}^{t}P_{dch_{i},\text{max}}^{t}$
$\displaystyle\text{for every }i\in\mathcal{N}_{S},t\in T,$ (15f)
$\displaystyle\mu_{ch_{i}}^{t}+\mu_{dch_{i}}^{t}\leq 1$
$\displaystyle\text{for every }i\in\mathcal{N}_{S},t\in T.$ (15g)

Constraints (15a) and (15b) set the upper/lower bounds of soc and give the
relationships between soc and charging/discharging actions. Constraints (15c)
– (15g) depict the connection between storage actions $\mu_{ch}$ and
$\mu_{dch}$ and their corresponding real power output
$P_{ch_{i}}^{t},P_{dch_{i}}^{t}$.

Lastly, the SO’s cost function is given by

$c^{\top}y\coloneqq\sum_{t\in
T}\sum_{i\in\mathcal{N}_{G}}(r_{i}^{+,t}P_{G_{i}}^{+,t}+r_{i}^{-,t}P_{G_{i}}^{-,t})+p_{sl}^{t}(P_{sl}^{t}-P_{fl}^{t})+\sum_{i\in\mathcal{N}_{DG}}(f_{i}^{+}P_{DG_{i}}^{t,+}+f_{i}^{-}P_{DG_{i}}^{t,-}).$

The cost function aims to minimize the electricity cost and 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,h)}c^{\top}y.$

Following the structure from Section 2, we denote by $Ay\geq b$ the
constraints, that solely deal with third-level variables, i.e., every
constraint of (11a)-(15g) despite of (11d). We observe, that the only
remaining constraint is the upper bound on the second-level (adverserial)
variables $h=P_{DG,\max}$ given by (11d) implying
$B=\begin{pmatrix}0&-I_{DG}\end{pmatrix},B_{x}=0,B_{h}=-I_{DG,\max},b_{0}=0$.
Consequently, the third-level program in this particular case reads

$\displaystyle\min\ $ $\displaystyle c^{\top}y$ (16a) s.t. $\displaystyle
Ay\geq b,$ (16b) $\displaystyle-P_{DG_{i}}^{t}\geq-P_{DG_{i},\max}^{t}$
$\displaystyle\text{ for every }i\in\mathcal{N}_{DG},t\in T.$ (16c)

In addition, we denote by $a_{DG}$ the columns of $A$ corresponding to
$P_{DG}$ and the remaining columns by $A^{\prime}$, i.e.,
$A=\begin{bmatrix}A^{\prime}&a_{DG}\end{bmatrix}$. Let further $\alpha$ denote
the dual variables that correspond to Constraints (11a) – (15c), i.e. to
$Ay\geq b$ and $\beta\in\mathbb{R}^{\mathcal{N}_{DG}\times T}$ denote the dual
variables corresponding to (16c). Then the dual program of (16) is

$\displaystyle\max\ $ $\displaystyle
b^{\top}\alpha-\sum_{i\in\mathcal{N}_{DG},t\in
T}P_{DG_{i},\max}^{t}\beta_{DG_{i},t}$ (17a) s.t.
$\displaystyle(A^{\prime})^{\top}\alpha=c,$ (17b) $\displaystyle
a_{DG_{i}}^{\top}\alpha-\beta_{DG_{i},t}=0$ $\displaystyle\text{ for every
}i\in\mathcal{N}_{DG},t\in T,$ (17c) $\displaystyle\alpha,\beta\geq 0,$ (17d)

We observe, that, as in Section 2, on the second level the objective
decomposes into a bilinear part and a linear one. However, one can argue that
the dual variables $\beta_{DG}$ as well as the maximal capacity $P_{DG,\max}$
of the distributed generators are bounded, i.e., Assumption 1 holds. In
particular, we will show that there are $P_{i,t}^{+},\beta_{i,t}^{+}$ such
that

$P_{DG_{i},\min}^{t}\leq P_{DG_{i},\max}^{t}\leq P_{i,t}^{+}\text{ and
}\beta_{i,t}^{-}\leq\beta_{DG_{i},t}\leq\beta_{i,t}^{+}.$

On the one hand, this is because distributed generators have technical limits.
For instance, the power outputs of wind turbines are restricted by the cut-out
wind speed. The outputs will not exceed the power corresponding to this wind
speed. As for the solar panels, their outputs are also restricted by the rated
power of the devices themselves. Hence, $P_{DG_{i},\max}^{t}$ is always
bounded.

On the other hand, we may prove the existence of an upper bound for
$\beta_{DG_{i},t}$ and thereby verify Assumption 1 in our application as
follows:

######  Lemma 1.

For every optimal solution $(\alpha^{*},\beta^{*})$ to (17), we have that

$p_{sl}^{t}-f_{i}^{+}\leq\beta_{DG_{i},t}^{*}\leq\max\\{r_{i}^{+,t},r_{i}^{-,t},p_{sl}^{t}\\}-\min\\{f_{i}^{+},f_{i}^{-}\\}\text{
for every }i\in\mathcal{N}_{DG},t\in T.$

###### Proof.

Consider the following relaxed version of (16), where we penalized violations
in (16c) instead of incorporating (16c) as a hard constraint:

$\displaystyle\min\ $ $\displaystyle c^{\top}y+\sum_{i\in\mathcal{N}_{DG},t\in
T}\beta_{i,t}^{+}\lambda_{i,t}$ (18a) s.t. $\displaystyle Ay\geq b,$ (18b)
$\displaystyle-P_{DG_{i}}^{t}+\lambda_{i,t}\geq-P_{DG_{i},\max}^{t}$
$\displaystyle\text{for every }i\in\mathcal{N}_{DG},t\in T$ (18c)
$\displaystyle\lambda_{i,t}\geq 0$ $\displaystyle\text{for every
}i\in\mathcal{N}_{DG},t\in T,$ (18d)

where
$\beta_{i,t}^{+}\coloneqq\max\\{r_{i}^{+,t},r_{i}^{-,t},p_{sl}^{t}\\}-\min\\{f_{i}^{+},f_{i}^{-}\\}$.
Suppose $\lambda_{i,t}>0$, then we can increase the value of $P_{DG_{i}}^{t}$
by at most $\lambda_{i,t}$, thereby at least decreasing the objective value by
$\min\\{f_{i}^{+},f_{i}^{-}\\}\lambda_{i,t}$. Due to the power balance
equations (13a) and (13b), we have to either sell the energy on the (second-
level) market, i.e., decrease $P_{sl}^{t}$, which results in a benefit of
$p_{sl}^{t}\lambda_{i,t}$, (in-)decrease
$P_{ch_{i}},P_{dch_{i}}^{t},P_{G_{i},\text{reg}}$ resulting either in a
benefit of $0$ or at most $\max\\{r_{i}^{+,t},r_{i}^{-,t}\\}\lambda_{i,t}$
respectively. Given this $\beta_{i,t}^{+}$, we obtain by strong duality:

$\displaystyle c^{\top}y^{*}=\max\ $ $\displaystyle
b^{\top}\alpha-\sum_{i\in\mathcal{N}_{DG},t\in
T}P_{DG_{i},\max}^{t}\beta_{DG_{i},t}$ (19a) s.t.
$\displaystyle(A^{\prime})^{\top}\alpha=c,$ (19b) $\displaystyle
a_{DG_{i}}^{\top}\alpha-\beta_{DG_{i},t}=0$ $\displaystyle\text{ for every
}i\in\mathcal{N}_{DG},t\in T,$ (19c)
$\displaystyle\beta_{DG_{i},t}+\gamma_{i,t}=\beta_{i,t}^{+}$
$\displaystyle\text{ for every }i\in\mathcal{N}_{DG},t\in T,$ (19d)
$\displaystyle\alpha,\beta,\gamma\geq 0.$ (19e)

Here, the last two constraints imply

$\beta_{DG_{i},t}\leq\beta_{i,t}^{+}\text{ for every
}i\in\mathcal{N}_{DG},t\in T,$

i.e. adding a sufficiently large upper bound on $\beta_{DG}$ does not change
the outcome of the dual program and hence, we can safely assume
$0\leq\beta_{DG}\leq\beta^{+}$.

For the inequality $p_{sl}^{t}-f_{i}^{+}\leq\beta_{DG_{i},t}^{*}$, we observe
that (18) is unbounded whenever $\beta_{i,t}^{+}<p_{sl}^{t}-f_{i}^{+}$:
Consider an optimal solution $(y^{*},P_{DG}^{*})$ of (16), which is feasible
for (18) with $\lambda_{i,t}^{*}=0$ for every $i\in\mathcal{N}_{DG},t\in T$.
If we fix $k\in\mathcal{N}_{DG},t^{\prime}\in T$, then we observe that
$P_{DG_{k}}^{t^{\prime}}(\mu)\coloneqq(P_{DG_{k}}^{t^{\prime}})^{*}+\mu$,
$P_{DG_{k}}^{+,t^{\prime}}(\mu)\coloneqq(P_{DG_{k}}^{+,t^{\prime}})^{*}+\mu$,
$\lambda_{k,t^{\prime}}(\mu)\coloneqq\lambda_{k,t^{\prime}}^{*}+\mu$,
$P_{sl}^{t^{\prime}}(\mu)\coloneqq(P_{sl}^{t^{\prime}})^{*}-\mu$ is also
feasible for (18) for every $\mu>0$. Its objective value is

$c^{\top}y^{*}-p_{sl}^{t^{\prime}}\mu+f_{k}^{+}\mu+\sum_{i\in\mathcal{N}_{DG},t\in
T}\beta_{i,t}^{+}\lambda_{i,t}^{*}+\beta_{k,t^{\prime}}^{+}\mu,$

which tends to $-\infty$ for $\mu\rightarrow\infty$, if
$\beta_{k,t^{\prime}}^{+}<p_{sl}^{t^{\prime}}-f_{k}^{+}$. By strong duality,
we obtain, that for $\beta_{k,t^{\prime}}^{+}<p_{sl}^{t^{\prime}}-f_{k}^{+}$
(19) does not have a feasible solution. Since $k,t^{\prime}$ have been chosen
arbitrarily, the following problem is equivalent to (19):

$\displaystyle c^{\top}y^{*}=\max\ $ $\displaystyle
b^{\top}\alpha-\sum_{i\in\mathcal{N}_{DG},t\in
T}P_{DG_{i},\max}^{t}\beta_{DG_{i},t}$ (20a) s.t.
$\displaystyle(A^{\prime})^{\top}\alpha=c,$ (20b) $\displaystyle
a_{DG_{i}}^{\top}\alpha-\beta_{DG_{i},t}=0$ $\displaystyle\text{ for every
}i\in\mathcal{N}_{DG},t\in T,$ (20c)
$\displaystyle\beta_{DG_{i},t}+\gamma_{i,t}=\max\\{r_{i}^{+,t},r_{i}^{-,t},p_{sl}^{t}\\}-\min\\{f_{i}^{+},f_{i}^{-}\\}$
$\displaystyle\text{ for every }i\in\mathcal{N}_{DG},t\in T,$ (20d)
$\displaystyle\beta_{DG_{i},t}-\delta_{i,t}=p_{sl}^{t}-f_{i}^{+}$
$\displaystyle\text{ for every }i\in\mathcal{N}_{DG},t\in T,$ (20e)
$\displaystyle\alpha,\beta,\gamma,\delta\geq 0.$ (20f)

∎

Moreover, this bound might even be sharp due to the following observation:

######  Remark 1.

Depending on the market operations, that is considered, we may have that both,
$r_{i}^{+,t},r_{i}^{-,t}\leq p_{sl}^{t}$ and $f_{i}^{+}\leq f_{i}^{-}$ holds.
In this case the two bounds in Lemma 1 coincide and

$\beta_{DG_{i},t}^{*}=p_{sl}^{t}-f_{i}^{+}.$

Thus, in those cases the McCormick envelopes are exact and the upcoming
corollaries provide an exact reformulation of (2). However, even if the
McCormick envelope is not exact, the boundedness of both, $\beta_{DG}$ and
$P_{DG,\max}$ enables us to relax the second level problem
$\max_{h\in\Omega}D(x,h)$ with the McCormick envelope. This gives rise to the
following corollary of Theorem 1:

######  Corollary 1.

Let $\Omega$ be a robust ambiguity set as in (10), where $\Omega$ denotes a
polytope with potential nonnegative slack variables $\eta$ in the rows $i\in
I$, i.e. $\Omega=\left\\{P_{DG,\max}\in\mathbb{R}^{\mathcal{N}_{DG}\times
T},\eta\in\mathbb{R}^{I}_{\geq
0}:A_{\Omega}^{\top}P_{DG,\max}+\eta=b_{\Omega}\right\\}$. Let further
$\beta_{i,t}^{-},\beta_{i,t}^{+}$ be a lower/upper bound for
$\beta_{DG_{i},t}$. Then, the following linear program provides an upper bound
to (2):

$\displaystyle\min\ $ $\displaystyle\sum_{i\in\mathcal{N}_{G},t\in
T}c_{G_{i}^{t},2}(P_{G_{i}}^{t})^{2}+c_{G_{i}^{t},1}P_{G_{i}}^{t}+c_{G_{i}^{t},0}+\sum_{t\in
T}p_{fl}^{t}P_{fl}^{t}+c^{\top}y$
$\displaystyle+\sum_{i\in\mathcal{N}_{DG},t\in
T}\beta_{i,t}^{+}u_{i,t}^{\beta,+}+\sum_{i\in\mathcal{N}_{DG},t\in
T}\beta_{i,t}^{-}u_{i,t}^{\beta,-}+b_{\Omega}^{\top}u_{\Omega}$
$\displaystyle-\sum_{i\in\mathcal{N}_{DG},t\in
T}P_{DG_{i},\min}^{t}\beta_{i,t}^{-}(u_{\text{env}})_{i,t}^{1}-\sum_{i\in\mathcal{N}_{DG},t\in
T}P_{i,t}^{+}\beta_{i,t}^{+}(u_{\text{env}})_{i,t}^{2}$
$\displaystyle-\sum_{i\in\mathcal{N}_{DG},t\in
T}P_{i,t}^{+}\beta_{i,t}^{-}(u_{\text{env}})_{i,t}^{3}-\sum_{i\in\mathcal{N}_{DG},t\in
T}P_{DG_{i},\min}^{t}\beta_{i,t}^{+}(u_{\text{env}})_{i,t}^{4}$ (21a) s.t.
$\displaystyle Ay\geq b,$ (21b) $\displaystyle-
P_{DG_{i}}^{t}+u_{i,t}^{\beta,+}+u_{i,t}^{\beta,-}-P_{DG_{i},\min}^{t}(u_{\text{env}})_{i,t}^{1}-P_{i,t}^{+}(u_{\text{env}})_{i,t}^{2}$
$\displaystyle\qquad-
P_{i,t}^{+}(u_{\text{env}})_{i,t}^{3}-P_{DG_{i},\min}^{t}(u_{\text{env}})_{i,t}^{4}\geq
0$ $\displaystyle\forall i\in\mathcal{N}_{DG},t\in T$ (21c) $\displaystyle
u_{i,t}^{\beta,+}\geq 0,\ -u_{i,t}^{\beta,-}\geq 0$ $\displaystyle\forall
i\in\mathcal{N}_{DG},t\in T$ (21d)
$\displaystyle(A_{\Omega}u_{\Omega})_{i}-\beta_{i,t}^{-}((u_{\text{env}})_{i,t}^{1}+(u_{\text{env}})_{i,t}^{3})-\beta_{i,t}^{+}((u_{\text{env}})_{i,t}^{2}+(u_{\text{env}})_{i,t}^{4})\geq
0$ $\displaystyle\forall i\in\mathcal{N}_{DG},t\in T$ (21e) $\displaystyle
u_{\Omega}\geq 0$ (21f)
$\displaystyle(u_{\text{env}})_{i,t}^{1}+(u_{\text{env}})_{i,t}^{2}+(u_{\text{env}})_{i,t}^{3}+(u_{\text{env}})_{i,t}^{4}\geq-1$
$\displaystyle\forall i\in\mathcal{N}_{DG},t\in T$ (21g)
$\displaystyle-(u_{\text{env}})_{i,t}^{1},-(u_{\text{env}})_{i,t}^{2}\geq 0$
$\displaystyle\forall i\in\mathcal{N}_{DG},t\in T$ (21h)
$\displaystyle(u_{\text{env}})_{i,t}^{3},(u_{\text{env}})_{i,t}^{4}\geq 0$
$\displaystyle\forall i\in\mathcal{N}_{DG},t\in T$ (21i)
$\displaystyle\eqref{Constr: fl1_define_Pfl},\ \&\ \eqref{Constr:
fl2_bounds_P_G}$ (21j)

We note, that Corollary 1 is a direct consequence of Theorem 1 and thus its
proof follows the same lines. However, we include the full proof here, as the
notation varies a bit and the proof illustrates the impact of Lemma 1.

###### Proof.

We observe that with (17) and Lemma 1 the second-level
$\max_{P_{DG,\max}\in\Omega}\min_{y\in\mathcal{Y}(x,h)}c^{\top}y$ can be
written as

$\displaystyle\max\ $ $\displaystyle
b^{\top}\alpha-\sum_{i\in\mathcal{N}_{DG},t\in
T}P_{DG_{i},\max}^{t}\beta_{DG_{i},t}$ (22a) s.t.
$\displaystyle(A^{\prime})^{\top}\alpha=c,$ (22b) $\displaystyle
a_{DG_{i}}^{\top}\alpha-\beta_{DG_{i},t}=0$ $\displaystyle\text{ for every
}i\in\mathcal{N}_{DG},t\in T,$ (22c)
$\displaystyle\beta_{DG_{i},t}+\gamma_{i,t}=\beta_{i,t}^{+}$
$\displaystyle\text{ for every }i\in\mathcal{N}_{DG},t\in T,$ (22d)
$\displaystyle\beta_{DG_{i},t}-\delta_{i,t}=\beta_{i,t}^{-}$
$\displaystyle\text{ for every }i\in\mathcal{N}_{DG},t\in T,$ (22e)
$\displaystyle A_{\Omega}^{\top}P_{DG,\max}+\eta=b_{\Omega},$ (22f)
$\displaystyle\alpha,\beta,\gamma,P_{DG,\max},\delta,\eta\geq 0,$ (22g)

where $\beta^{-},\beta^{+}$ are chosen as in Lemma 1. Next, we substitute
$\kappa_{i,t}\coloneqq P_{DG_{i},\max}^{t}\beta_{DG_{i},t}$ in the objective
and relax the resulting constraint $\kappa_{i,t}\coloneqq
P_{DG_{i},\max}^{t}\beta_{DG_{i},t}$ by a McCormick envelope. Note, that since
$P_{DG_{i},\max}^{t},\beta_{DG_{i},t}\geq 0$, we can immediately conclude
$\kappa\geq 0$, which simplifies our notation a bit. If we further introduce
suitable nonnegative slack variables, we obtain the following dual LP:

$\displaystyle\max\ $ $\displaystyle
b^{\top}\alpha-\sum_{i\in\mathcal{N}_{DG},t\in
T}P_{DG_{i},\max}^{t}\beta_{DG_{i},t}$ (23a) s.t.
$\displaystyle(A^{\prime})^{\top}\alpha=c,$ (23b) $\displaystyle
a_{DG_{i}}^{\top}\alpha-\beta_{DG_{i},t}=0$ $\displaystyle\text{ for every
}i\in\mathcal{N}_{DG},t\in T,$ (23c)
$\displaystyle\beta_{DG_{i},t}+\gamma_{i,t}=\beta_{i,t}^{+}$
$\displaystyle\text{ for every }i\in\mathcal{N}_{DG},t\in T,$ (23d)
$\displaystyle\beta_{DG_{i},t}-\delta_{i,t}=\beta_{i,t}^{-}$
$\displaystyle\text{ for every }i\in\mathcal{N}_{DG},t\in T,$ (23e)
$\displaystyle A_{\Omega}^{\top}P_{DG,\max}+\eta=b_{\Omega},$ (23f)
$\displaystyle\kappa_{i,t}=P_{DG_{i},\min}^{t}\beta_{DG_{i},t}+P_{DG_{i},\max}^{t}\beta_{i,t}^{-}-P_{DG_{i},\min}^{t}\beta_{i,t}^{-}+\eta_{i,t}^{1}$
$\displaystyle\text{for every }i\in\mathcal{N}_{DG},t\in T,$ (23g)
$\displaystyle\kappa_{i,t}=P_{i,t}^{+}\beta_{DG_{i},t}+P_{DG_{i},\max}^{t}\beta_{i,t}^{+}-P_{i,t}^{+}\beta_{i,t}^{+}+\eta_{i,t}^{2}=0$
$\displaystyle\text{for every }i\in\mathcal{N}_{DG},t\in T,$ (23h)
$\displaystyle\kappa_{i,t}=P_{i,t}^{+}\beta_{DG_{i},t}+P_{DG_{i},\max}^{t}\beta_{i,t}^{-}-P_{i,t}^{+}\beta_{i,t}^{-}-\eta_{i,t}^{3}$
$\displaystyle\text{for every }i\in\mathcal{N}_{DG},t\in T,$ (23i)
$\displaystyle\kappa_{i,t}=P_{DG_{i},\max}^{t}\beta_{i,t}^{+}+P_{DG_{i},\min}\beta_{DG_{i},t}-P_{DG_{i},\min}^{t}\beta_{i,t}^{+}-\eta_{i,t}^{4}$
$\displaystyle\text{for every }i\in\mathcal{N}_{DG},t\in T,$ (23j)
$\displaystyle\alpha,\beta,\gamma,\delta,P_{DG,\max},\rho,\eta,\kappa\geq 0.$
(23k)

If we denote the dual variables of (23a) – (23h) by
$y,u_{i,t}^{\beta,+},u_{i,t}^{\beta,-},u_{\Omega},(u_{\text{env}})_{i,t}^{1},(u_{\text{env}})_{i,t}^{2},(u_{\text{env}})_{i,t}^{3},(u_{\text{env}})_{i,t}^{4}$
respectively, then the result follows by strong duality and including the
first-level variables and objectives. ∎

Again, we observe that Corollary 1 provides an LP inner approximation of the
linear relaxation of (1). Now, the following direct corollary of Theorem 2
incorporates the discrete (binary) decisions $\mu_{ch}$ and $\mu_{dch}$.

######  Corollary 2.

Let $\Omega$ be a robust ambiguity set as in (10), where $\Omega$ denotes a
polytope with potential nonnegative slack variables $\eta$ in the rows $i\in
I$, i.e. $\Omega=\left\\{P_{DG,\max}\in\mathbb{R}^{\mathcal{N}_{DG}\times
T},\eta\in\mathbb{R}^{I}_{\geq
0}:A_{\Omega}^{\top}P_{DG,\max}+\eta=b_{\Omega}\right\\}$. Let further
$\beta_{i,t}^{-},\beta_{i,t}^{+}$ be a lower/upper bound for
$\beta_{DG_{i},t}$. Then, the following MIP provides an upper bound to the
tri-level MIP (1) with the given parameters from Section 3:

$\displaystyle\min\ $ $\displaystyle\sum_{i\in\mathcal{N}_{G},t\in
T}c_{G_{i}^{t},2}(P_{G_{i}}^{t})^{2}+c_{G_{i}^{t},1}P_{G_{i}}^{t}+c_{G_{i}^{t},0}+\sum_{t\in
T}p_{fl}^{t}P_{fl}^{t}+c^{\top}y$
$\displaystyle+\sum_{i\in\mathcal{N}_{DG},t\in
T}\beta_{i,t}^{+}u_{i,t}^{\beta,+}+\sum_{i\in\mathcal{N}_{DG},t\in
T}\beta_{i,t}^{-}u_{i,t}^{\beta,-}+b_{\Omega}^{\top}u_{\Omega}$
$\displaystyle-\sum_{i\in\mathcal{N}_{DG},t\in
T}P_{DG_{i},\min}^{t}\beta_{i,t}^{-}(u_{\text{env}})_{i,t}^{1}-\sum_{i\in\mathcal{N}_{DG},t\in
T}P_{i,t}^{+}\beta_{i,t}^{+}(u_{\text{env}})_{i,t}^{2}$
$\displaystyle-\sum_{i\in\mathcal{N}_{DG},t\in
T}P_{i,t}^{+}\beta_{i,t}^{-}(u_{\text{env}})_{i,t}^{3}-\sum_{i\in\mathcal{N}_{DG},t\in
T}P_{DG_{i},\min}^{t}\beta_{i,t}^{+}(u_{\text{env}})_{i,t}^{4}$ (24a) s.t.
$\displaystyle Ay\geq b,$ (24b) $\displaystyle-
P_{DG_{i}}^{t}+u_{i,t}^{\beta,+}+u_{i,t}^{\beta,-}-P_{DG_{i},\min}^{t}(u_{\text{env}})_{i,t}^{1}-P_{i,t}^{+}(u_{\text{env}})_{i,t}^{2}$
$\displaystyle\qquad-
P_{i,t}^{+}(u_{\text{env}})_{i,t}^{3}-P_{DG_{i},\min}^{t}(u_{\text{env}})_{i,t}^{4}\geq
0$ $\displaystyle\forall i\in\mathcal{N}_{DG},t\in T$ (24c) $\displaystyle
u_{i,t}^{\beta,+}\geq 0,\ -u_{i,t}^{\beta,-}\geq 0$ $\displaystyle\forall
i\in\mathcal{N}_{DG},t\in T$ (24d)
$\displaystyle(A_{\Omega}u_{\Omega})_{i}-\beta_{i,t}^{-}((u_{\text{env}})_{i,t}^{1}+(u_{\text{env}})_{i,t}^{3})-\beta_{i,t}^{+}((u_{\text{env}})_{i,t}^{2}+(u_{\text{env}})_{i,t}^{4})\geq
0$ $\displaystyle\forall i\in\mathcal{N}_{DG},t\in T$ (24e) $\displaystyle
u_{\Omega}\geq 0$ (24f)
$\displaystyle(u_{\text{env}})_{i,t}^{1}+(u_{\text{env}})_{i,t}^{2}+(u_{\text{env}})_{i,t}^{3}+(u_{\text{env}})_{i,t}^{4}\geq-1$
$\displaystyle\forall i\in\mathcal{N}_{DG},t\in T$ (24g)
$\displaystyle-(u_{\text{env}})_{i,t}^{1},-(u_{\text{env}})_{i,t}^{2}\geq 0$
$\displaystyle\forall i\in\mathcal{N}_{DG},t\in T$ (24h)
$\displaystyle(u_{\text{env}})_{i,t}^{3},(u_{\text{env}})_{i,t}^{4}\geq 0$
$\displaystyle\forall i\in\mathcal{N}_{DG},t\in T$ (24i)
$\displaystyle\eqref{Constr: fl1_define_Pfl},\ \&\ \eqref{Constr:
fl2_bounds_P_G}$ (24j)
$\displaystyle\mu_{ch},\mu_{dch}\in\\{0,1\\}^{\mathcal{N}_{S}\times T}$ (24k)
$\displaystyle
P_{G_{i}}^{t},P_{fl}^{t},y^{\prime},P_{DG_{i}}^{t},u^{\beta,+}_{i,t},u^{\beta,-}_{i,t},u^{P}_{i,t},(u_{\text{env}})^{1}_{i,t},(u_{\text{env}})^{2}_{i,t},(u_{\text{env}})^{3}_{i,t},(u_{\text{env}})^{4}_{i,t}\in\mathbb{R}.$
(24l)

Thus, our approximation technique is applicable to the optimization of smart
converters in power grids. Moreover, the only strict relaxation comes from
approximating the uncertainty set $\Omega$ by McCormick envelopes of the
bilinear terms and may, depending on the second-level (intra-day) market price
$p_{sl}$ and the penalizations for adjustments $f^{+},f^{-},r^{+},r^{-}$, even
be sharp. Hence, it seems natural to test our approximations numerically.

##  4\. Computational results

We present the results of three case studies in this section. All of these
instances are considered on a daily basis divided into hourly (24 period) or
15min (96 period) time intervals. The first benchmark is a 5-bus instance,
based on the “case5.m” instance from the matpower library [ZMT11]. Second, we
consider a 30-bus instance, based on the “case_ieee30.m” instance of the
matpower library and finally a modified version of the IEEE 118-bus system
similar to the one in [Cob+18]. 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.

###  4.1. 5-bus example

The topology of the test system is shown in Figure 1, where we consider bus
$1$ to be the root node. We observe, that the conventional 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\\}$.

Figure 1. The case5.m network with its corresponding generators and flows at
3am. $P^{d}=0$ at buses $1$ and $5$

Whether a generator in “case5.m” is a conventional/renewable one or a storage
was decided by the authors. Both, the day-ahead and intra-day market prices
$p_{fl},p_{sl}$ were taken as averages from the Pecan street data base’s
[pec22] “miso” data set for October 3rd, 2022. The daily deviations in
$P^{d}$, denoted by $\Delta_{d}^{t}$, or daily deviations in
$P_{DG,\max}^{t}$, denoted by $\Delta_{DG}^{t}$ were similarly taken from the
Pecan street data base’s “california_iso” dataset for October 3rd, 2022.

Then, the demand varying over the day was modeled as
$P^{d}=P^{d}\cdot\Delta_{d}^{t}$ and the varying potential renewable energy
production was modeled by
$P_{DG,\text{forecast}}^{t}\coloneqq\min\\{\frac{P_{DG,\min}+P_{DG,\max}}{2}\cdot\Delta_{DG}^{t},P_{i,t}^{+}\\}$.

Solving Problem (24) takes less than 1s. However, its objective value highly
depends on the uncertainty imposed on the system. We illustrate in Figure 2
how sensitive the optimal solution reacts to changes in the maximal forecast
error $R$, that crucially determines $\Omega$ through (10). In particular,
since in “case5.m” we have neglectable upward and downward regulation costs
and $f^{-}\geq f^{+}$, Lemma 1 implies that the McCormick relaxation is sharp
if and only if $r^{-},r^{+}\leq p_{sl}$. As realistic penalties $r^{+}=r^{-}$,
we assume $r_{1}^{+}=r_{1}^{-}=14$$ p.u., which are the costs of operating the
first generator at bus $1$. For the sake of a better analysis, we replace the
natural choice $r_{4}^{+}=r_{4}^{-}=40$$ p.u., which are the costs of
operating the generator at bus $4$ by $r_{4}^{+}=r_{4}^{-}=20$$ p.u. as then
$r^{-},r^{+}\leq p_{sl}$ and we can compare variations in the penalties to an
optimal robust solution.

In particular, for this instance, $R=1$ yields the nominal optimal solution
with an objective value of $727082\$$.

Figure 2. Objective value of (24) on “case5.m” under varying $R$

Moreover, the blue line in Figure 2 illustrates the perfect linear relation
between the lower bound of the uncertainty set $\Omega$ and the objective
value of (24). We conclude, that the DSO may cut its worst-case costs by
almost $50\%$ with a perfectly accurate weather forecast, i.e.
$\Omega=\\{P_{DG,forecast}\\}$. As this is unrealistic with present
forecasting methods, we would like to highlight that one may gain already
significant cost reductions in the worst-case by incorporating more
information on $\Omega$.

Moreover, the red line in Figure 2 shows an upper bound to (2) given by
Corollary 2 in case $r^{+,t}_{4}>p_{sl}^{t}$ for some $t$, i.e. in case Lemma
1 is violated. Since the objective value with respect to this penalty is
contained between the red and the blue line, Figure 2 thereby shows a rather
strong approximation quality for this particular instance. However, we want to
stress, that this only holds for this particular example and may not be a
general pattern.

###  4.2. 30-bus example

Similarly, as in the 5-bus example, the topology of the 30-bus test system is
taken from “case_ieee30.m”, a system with 41 transmission lines and after
modification 4 dispatchable generators as well as 4 energy storages. The only
renewable generator is placed at bus $2$, i.e. $\mathcal{N}_{DG}=\\{2\\}$. In
addition, we choose $\mathcal{N}_{G}=\\{5,8,11,13\\}$ and the four energy
storage units to be connected to buses $1,2,8,13$, i.e.
$\mathcal{N}_{S}=\\{1,2,8,13\\}$. The estimation procedure of market prices
and demands are kept from the “case5.m” example. Since also “case_ieee30.m”
does not include upward or downward regulation costs for the generators, the
McCormick envelope is sharp and for $R=1$, i.e.
$\Omega=\\{P_{DG,forecast}\\}$, the nominal value of $104,088\$$ is attained.

Moreover, we would like to illustrate the dependency of the worst-case revenue
with respect to choices of $R$ in Figure 3. We want to highlight, that also in
this more elaborate example, the runtime was $<1$s.

Figure 3. Objective values of (24) with $|\mathcal{N}_{S}|\in\\{0,2,4\\}$ on
“case_ieee30.m”

Note, that different slopes may occur due to the different capabilities of the
storages. In particular, the improved performance from
$\mathcal{N}_{S}=\emptyset$ to $\mathcal{N}_{S}\neq\emptyset$ indicates that
the capability of storing all renewable energy produced within the
transmission system is more valuable than storing energy from prior purchases,
i.e. externally produced energy.

###  4.3. Analyzing the runtime on larger instances

After illustrating the behavior of the objective value under different
uncertainties and storages, we focus on the main advantage of the proposed MIP
approach, namely its speed. As the scaling of the runtime is crucial in
industrial applications, we demonstrate the applicability of our algorithm to
larger power systems, particularly a 118-bus (“case118.m”), a 200-bus
(“case200.m”) and 300-bus (“case300.m”) test system.

To this end, we aim to keep the considered instances as comparable as we can.
In particular, we again keep the estimation procedure of market prices and
demands from the “case5.m” example. Additionally, none of the considered
instances contains upward or downward regulation costs for the generators and
thereby due to Lemma 1 we always compute the exact solutions to (2). Lastly,
we note that the number of storages $|\mathcal{N}_{S}|$ determines the amount
of binary variables in (24) and consequently should crucially impacts the
runtime. Thus, we equipped our test systems with $|\mathcal{N}_{S}|=6$,
$|\mathcal{N}_{S}|=10$ and $|\mathcal{N}_{S}|=15$ storages respectively in
order to achieve an approximately constant ratio of buses to storages
$|V|/|\mathcal{N}_{S}|$, i.e. $|V|/|\mathcal{N}_{S}|\approx 20$. To further
improve comparability, we also recomputed the $30$-bus system with $2$
storages instead of $4$. The following figure illustrates the achieved
runtime:

Figure 4. runtime comparison “case5.m” with $|\mathcal{N}_{S}|=1$,
“case_ieee30.m” with $|\mathcal{N}_{S}|=2$, “case118.m” with
$|\mathcal{N}_{S}|=6$, “case200.m” with $|\mathcal{N}_{S}|=10$, “case300.m”
with $|\mathcal{N}_{S}|=15$

As (2) is a notoriously challenging problem, see Question in [CW22] and
Section 6 in [YGH19] benchmarks on the exact problem setting are, to the best
of our knowledge, rare. However, Cobos et. al. [Cob+18] applied a nested
column generation approach in order to solve a strongly related problem on
“case118.m” and achieved runtimes between $200s$ and $800s$. The considered
instances in [Cob+18] contain $|\mathcal{N}_{G}||T|+3|\mathcal{N}_{S}||T|$
first-level, $2|\mathcal{N}_{DG}||T|$ second-level and $|\mathcal{N}_{S}||T|$
third-level binary variables – a significantly larger amount of binary
variables compared to our instances since in [Cob+18] we have
$|\mathcal{N}_{G}|=54,|\mathcal{N}_{DG}|=10,|\mathcal{N}_{S}|=6$ and $|T|=24$.
In summary, the nested column generation in [Cob+18] addresses $54\cdot
24+3\cdot 6\cdot 24+2\cdot 10\cdot 24+6\cdot 24=98\cdot 24=2352$ binary
decisions on a comparable instance. As Figure 4 illustrates, the proposed
algorithm solves the instance in $\approx 3$s, but with significantly fewer,
namely $2\cdot 6\cdot 24=12\cdot 24=288$ binary decisions. Thus, a direct
comparison with [Cob+18] seems rather inappropriate.

Furthermore, the parametric programming approach presented in [AP19] can be
used to solve the ARO (2), even if (2) is not weakly-connected. The same
authors demonstrate, that an instance with $60$ binary variables on various
levels can be solved within $15$s, see Table 7, Problem P5 in [AP20]. As the
number of binary variables of this instance is still significantly smaller
than “case_ieee30.m”, which the presented MIP framework can solve within
$<1$s, it seems natural to conjecture, that the MIP approach outperforms the
parametric programming approach on weakly-connected AROs in terms of runtime.
However, we would like to highlight, that the work in [AP19] rather aims at
wide applicability as the authors present a significantly more general
approach.

##  5\. Conclusion

This article presents a new MIP framework to approximate adjustable robust
programs with integer variables in the innermost (adjustment) stage. It is
based on a weak connection between the separate stages and uses a McCormick
envelope to strengthen the adversarial, thereby relaxing the ARO. We have
proven that the resulting MIP provides feasible solutions for the first stage,
that can be adjusted to a solution satisfying an objective at least as good as
the ARO objective regardless of the realization of the uncertainty. In
addition, we have provided a sufficient criterion, for the exactness of our
approximation.

Moreover, we applied our results to model discrete adjustments of smart
converters in a power system and provided numerical evidence, that our
approach is competitive to previous methods such as the nested column
generation or parametric programming in terms of runtime.

## Acknowledgments

We are grateful to Robert Burlacu and Qianwen Xu for stimulating discussions.
The authors gratefully acknowledge Digital Futures for financially supporting
the project.

## References

  * [AC17] Miguel Anjos and Antonio J Conejo  “Unit Commitment in Electric Energy Systems”  Now FoundationsTrends, 2017  DOI: 10.1561/3100000014
  * [AP19] Styliani Avraamidou and Efstratios N. Pistikopoulos  “Multi-parametric global optimization approach for tri-level mixed-integer linear optimization problems”  In _J Glob Optim_ 74, 2019, pp. 443–465  DOI: https://doi.org/10.1007/s10898-018-0668-4
  * [AP20] Styliani Avraamidou and Efstratios N. Pistikopoulos  “Adjustable robust optimization through multi-parametric programming”  In _Optimization Letters_ 14, 2020, pp. 873–887  DOI: https://doi.org/10.1007/s11590-019-01438-5
  * [Aßm+18] Denis Aßmann, Frauke Liers, Michael Stingl and Juan C. Vera  “Deciding Robust Feasibility and Infeasibility Using a Set Containment Approach: An Application to Stationary Passive Gas Network Operations”  In _SIAM Journal on Optimization_ 28.3, 2018, pp. 2489–2517  DOI: 10.1137/17M112470X
  * [BBC11] Dimitris Bertsimas, David B Brown and Constantine Caramanis  “Theory and applications of robust optimization”  In _SIAM review_ 53.3 SIAM, 2011, pp. 464–501
  * [Ben+04] Aharon Ben-Tal, Alexander Goryashko, Elana Guslitzer and Arkadi Nemirovski  “Adjustable robust solutions of uncertain linear programs”  In _Mathematical programming_ 99.2 Springer, 2004, pp. 351–376
  * [Ben+05] Aharon Ben-Tal, Boaz Golany, Arkadi Nemirovski and Jean-Philippe Vial  “Retailer-supplier flexible commitments contracts: A robust optimization approach”  In _Manufacturing & Service Operations Management_ 7.3 INFORMS, 2005, pp. 248–271
  * [Ber+13] Dimitris Bertsimas et al.  “Adaptive Robust Optimization for the Security Constrained Unit Commitment Problem”  In _IEEE Transactions on Power Systems_ 28.1, 2013, pp. 52–63  DOI: 10.1109/TPWRS.2012.2205021
  * [BG15] Dimitris Bertsimas and Angelos Georghiou  “Design of near optimal decision rules in multistage adaptive mixed-integer optimization”  In _Operations Research_ 63.3 INFORMS, 2015, pp. 610–627
  * [BH17] Viktoryia Buhayenko and Dick Hertog  “Adjustable robust optimisation approach to optimise discounts for multi-period supply chain coordination under demand uncertainty”  In _International Journal of Production Research_ 55.22 Taylor & Francis, 2017, pp. 6801–6823
  * [Bha+15] Binayak Bhandari et al.  “Optimization of hybrid renewable energy power systems: A review”  In _International journal of precision engineering and manufacturing-green technology_ 2 Springer, 2015, pp. 99–112
  * [BL11] J.R. Birge and F. Louveaux  “Introduction to Stochastic Programming”, Springer Series in Operations Research and Financial Engineering  Springer New York, 2011  URL: https://books.google.se/books?id=Vp0Bp8kjPxUC
  * [BÖ08] Daniel Bienstock and Nuri Özbay  “Computing robust basestock levels” In Memory of George B. Dantzig  In _Discrete Optimization_ 5.2, 2008, pp. 389–414  DOI: https://doi.org/10.1016/j.disopt.2006.12.002
  * [Cob+18] Noemi G. Cobos, José M. Arroyo, Natalia Alguacil and Jianhui Wang  “Robust Energy and Reserve Scheduling Considering Bulk Energy Storage Units and Wind Uncertainty”  In _IEEE Transactions on Power Systems_ 33.5, 2018, pp. 5206–5216  DOI: 10.1109/TPWRS.2018.2792140
  * [CW22] A.J. Conejo and X. Wu  “Robust optimization in power systems: a tutorial overview”  In _Optim Eng_ 23, 2022, pp. 2051–2073  DOI: https://doi.org/10.1007/s11081-021-09667-3
  * [Fan53] Ky Fan  “Minimax Theorems”  In _Proceedings of the National Academy of Sciences of the United States of America_ 39, 1953, pp. 42–47
  * [GCG14] Ignacio E Grossmann, Bruno A Calfa and Pablo Garcia-Herreros  “Evolution of concepts and models for quantifying resiliency and flexibility of chemical processes”  In _Computers & Chemical Engineering_ 70 Elsevier, 2014, pp. 22–34
  * [Gro+16] Ignacio E Grossmann et al.  “Recent advances in mathematical programming techniques for the optimization of process systems under uncertainty”  In _Computers & Chemical Engineering_ 91 Elsevier, 2016, pp. 3–14
  * [INO20] Semich Impram, Secil Varbak Nese and Bülent Oral  “Challenges of renewable energy penetration on power system flexibility: A survey”  In _Energy Strategy Reviews_ 31 Elsevier, 2020, pp. 100539
  * [Lap+19] Nikos H Lappas, Luis A Ricardez-Sandoval, Ricardo Fukasawa and Chrysanthos E Gounaris  “Adjustable Robust Optimization for multi-tasking scheduling with reprocessing due to imperfect tasks”  In _Optimization and Engineering_ 20.4 Springer, 2019, pp. 1117–1159
  * [Lek+21] Amir Aris Lekvan et al.  “Robust optimization of renewable-based multi-energy micro-grid integrated with flexible energy conversion and storage devices”  In _Sustainable Cities and Society_ 64 Elsevier, 2021, pp. 102532
  * [Lem01] Claude Lemaréchal  “Lagrangian Relaxation”  In _Computational Combinatorial Optimization_ , 2001
  * [LG16] Nikolaos H Lappas and Chrysanthos E Gounaris  “Multi-stage adjustable robust optimization for process scheduling under uncertainty”  In _AIChE Journal_ 62.5 Wiley Online Library, 2016, pp. 1646–1667
  * [LI08] Zukui Li and Marianthi G Ierapetritou  “Robust optimization for process scheduling under uncertainty”  In _Industrial & Engineering Chemistry Research_ 47.12 ACS Publications, 2008, pp. 4148–4157
  * [LJF04] Xiaoxia Lin, Stacy L Janak and Christodoulos A Floudas  “A new robust optimization approach for scheduling under uncertainty:: I. Bounded uncertainty”  In _Computers & chemical engineering_ 28.6-7 Elsevier, 2004, pp. 1069–1085
  * [LS16] Alvaro Lorca and Xu Andy Sun  “Multistage robust unit commitment with dynamic uncertainty sets and energy storage”  In _IEEE Transactions on Power Systems_ 32.3 IEEE, 2016, pp. 1678–1688
  * [LS17] Álvaro Lorca and Xu Andy Sun  “Multistage Robust Unit Commitment With Dynamic Uncertainty Sets and Energy Storage”  In _IEEE Transactions on Power Systems_ 32.3, 2017, pp. 1678–1688  DOI: 10.1109/TPWRS.2016.2593422
  * [pec22] "pecanstreet.org"  In _[Online]_ Pecan street dataport, 2022  URL: https://www.pecanstreet.org/
  * [PHK14] Stefan Pfenninger, Adam Hawkes and James Keirstead  “Energy systems modeling for twenty-first century energy challenges”  In _Renewable and Sustainable Energy Reviews_ 33 Elsevier, 2014, pp. 74–86
  * [SM87] CW Sanders and CA Monroe  “An algorithm for real-time security constrained economic dispatch”  In _IEEE transactions on power systems_ 2.4 IEEE, 1987, pp. 1068–1074
  * [Tej+19] Manuel Tejeda-Iglesias, Nikolaos H Lappas, Chrysanthos E Gounaris and Luis Ricardez-Sandoval  “Explicit model predictive controller under uncertainty: An adjustable robust optimization approach”  In _Journal of Process Control_ 84 Elsevier, 2019, pp. 115–132
  * [TTT08] Akiko Takeda, Shunsuke Taguchi and RH Tütüncü  “Adjustable robust optimization models for a nonlinear two-period system”  In _Journal of Optimization Theory and Applications_ 136.2 Springer, 2008, pp. 275–295
  * [Wan+95] SJ Wang et al.  “Short-term generation scheduling with transmission and environmental constraints using an augmented Lagrangian relaxation”  In _IEEE Transactions on Power Systems_ 10.3 IEEE, 1995, pp. 1294–1301
  * [Wei+14] Wei Wei, Feng Liu, Shengwei Mei and Yunhe Hou  “Robust energy and reserve dispatch under variable renewable generation”  In _IEEE Transactions on Smart Grid_ 6.1 IEEE, 2014, pp. 369–380
  * [Wei+15] Stefan Weitemeyer, David Kleinhans, Thomas Vogt and Carsten Agert  “Integration of Renewable Energy Sources in future power systems: The role of storage”  In _Renewable Energy_ 75 Elsevier, 2015, pp. 14–20
  * [YGH19] İhsan Yanıkoğlu, Bram L. Gorissen and Dick Hertog  “A survey of adjustable robust optimization”  In _European Journal of Operational Research_ 277.3, 2019, pp. 799–813  DOI: https://doi.org/10.1016/j.ejor.2018.08.031
  * [YS21] Jun Yang and Changqi Su  “Robust optimization of microgrid based on renewable distributed power generation and load demand uncertainty”  In _Energy_ 223 Elsevier, 2021, pp. 120043
  * [Yu+20] Dongmin Yu et al.  “Energy management of wind-PV-storage-grid based large electricity consumer using robust optimization technique”  In _Journal of Energy Storage_ 27 Elsevier, 2020, pp. 101054
  * [YW19] Yue Yang and Wenchuan Wu  “A Distributionally Robust Optimization Model for Real-Time Power Dispatch in Distribution Networks”  In _IEEE Transactions on Smart Grid_ 10.4, 2019, pp. 3743–3752  DOI: 10.1109/TSG.2018.2834564
  * [ZG13] Yu Zhang and Georgios B Giannakis  “Robust optimal power flow with wind integration using conditional value-at-risk”  In _2013 IEEE International Conference on Smart Grid Communications (SmartGridComm)_ , 2013, pp. 654–659  IEEE
  * [ZG13a] Chaoyue Zhao and Yongpei Guan  “Unified stochastic and robust unit commitment”  In _IEEE Transactions on Power Systems_ 28.3 IEEE, 2013, pp. 3353–3361
  * [ZGL16] Qi Zhang, Ignacio E Grossmann and Ricardo M Lima  “On the relation between flexibility analysis and robust optimization for linear systems”  In _AIChE Journal_ 62.9 Wiley Online Library, 2016, pp. 3109–3123
  * [Zha+18] Xin Zhao, Liuchen Chang, Riming Shao and Katelin Spence  “Power system support functions provided by smart inverters—A review”  In _CPSS Transactions on Power Electronics and Applications_ 3.1 CPSS, 2018, pp. 25–35
  * [ZMT11] Ray Daniel Zimmerman, Carlos Edmundo Murillo-Sánchez and Robert John Thomas  “MATPOWER: Steady-State Operations, Planning, and Analysis Tools for Power Systems Research and Education”  In _IEEE Transactions on Power Systems_ 26.1, 2011, pp. 12–19  DOI: 10.1109/TPWRS.2010.2051168
  * [ZZ13] Bo Zeng and Long Zhao  “Solving two-stage robust optimization problems using a column-and-constraint generation method”  In _Operations Research Letters_ 41.5 Elsevier, 2013, pp. 457–461