###### Proof.

###### Data-Driven Two-Stage Conic Optimization with Zero-One Uncertainties

###### Abstract

We address high-dimensional zero-one random parameters in two-stage convex
conic optimization problems. Such parameters typically represent failures of
network elements and constitute rare, high-impact random events in several
applications. Given a sparse training dataset of the parameters, we motivate
and study a distributionally robust formulation of the problem using a
Wasserstein ambiguity set centered at the empirical distribution. We present a
simple, tractable, and conservative approximation of this problem that can be
efficiently computed and iteratively improved. Our method relies on a
reformulation that optimizes over the convex hull of a mixed-integer conic
programming representable set, followed by an approximation of this convex
hull using lift-and-project techniques. We illustrate the practical viability
and strong out-of-sample performance of our method on nonlinear optimal power
flow and multi-commodity network design problems that are affected by random
contingencies, and we report improvements of up to 20% over existing sample
average approximation and two-stage robust optimization methods.

##  1 Motivation

This work is motivated by optimization problems arising in applications that
are affected by an extremely large, yet finite, number of rare, high-impact
random events. In particular, we are motivated by applications in which the
decision-relevant random events consist of high-dimensional binary outcomes.
Such applications are ubiquitous in network optimization, where the uncertain
events amount to failures of the nodes or edges of the underlying network. For
example, in electric power networks, random node and edge failures have been
used to model losses of physical components such as substations, transmission
lines, generators, and transformers [13, 51]. Similarly, in natural gas [55],
wireless communication [27], and transportation networks [10], they can be
used to model failures of compressors and gas pipelines, antennas, and road
links, respectively.

The challenges in modeling and solving such uncertainty-affected optimization
problems are threefold. First, the number of possible failure states grows
_exponentially_ as the size of the input increases. For example, the number of
failure states in a network with $n$ failure-prone nodes is $2^{n}$. Second,
network failures are extremely _rare_ but _critical_ , and historical records
are often not rich enough to include observations for every possible failure
state. Third, failures of individual network elements are unlikely to be
independent of each other. For example, transmission line failures in electric
power systems often have a cascading effect that triggers the failure of other
transmission lines. Because of this _high dimensionality_ , the true
distribution governing the random parameters is often unknown and difficult to
estimate with a small number of historical observations.

In this high-dimensional context, we focus on _two-stage optimization_
problems under uncertainty. Suppose that there are $M$ uncertain parameters
$\tilde{\xi}_{1},\ldots,\tilde{\xi}_{M}$, and that $\Xi\subseteq\\{0,1\\}^{M}$
is the support of the underlying distribution of these parameters. If the
distribution $\mathbb{P}$ is known, then the two-stage problem takes the form

$\mathop{\text{minimize}}_{\bm{x}\in\mathcal{X}}\;c(\bm{x})+\mathbb{E}_{\mathbb{P}}\left[\mathcal{Q}(\bm{x},\tilde{\bm{\xi}})\right],$

where $\bm{x}$ represents the _first-stage_ decisions that must be made
agnostically to the realization of the random parameters,
$\mathcal{X}\subseteq\mathbb{R}^{N_{1}}$ is a convex compact set of feasible
first-stage decisions, $c:\mathcal{X}\mapsto\mathbb{R}$ is a convex function
representing the deterministic cost associated with these decisions, and
$\mathcal{Q}(\bm{x},\bm{\xi})$ is the _random loss_ (or _second-stage cost_)
corresponding to decisions $\bm{x}$ and a fixed realization $\bm{\xi}\in\Xi$
of the random parameters. We assume that the loss can be computed by solving a
convex conic optimization problem of the form

$\mathcal{Q}(\bm{x},\bm{\xi})=\inf_{\bm{y}\in\mathcal{Y}}\left\\{\bm{q}(\bm{\xi})^{\top}\bm{y}:\bm{W}(\bm{\xi})\bm{y}\geq\bm{T}(\bm{x})\bm{\xi}+\bm{h}(\bm{x})\right\\},$
(1)

where $\bm{y}$ denotes the second-stage decisions that can be made after the
realization $\bm{\xi}$ is known; $\mathcal{Y}\subseteq\mathbb{R}^{N_{2}}$ is a
proper (closed, convex, pointed, and full-dimensional) cone;
$\bm{q}:\Xi\mapsto\mathbb{R}^{N_{2}}$ and
$\bm{W}:\Xi\mapsto\mathbb{R}^{L\times N_{2}}$ are vector- and matrix-valued
affine functions respectively, while $\bm{h}:\mathcal{X}\mapsto\mathbb{R}^{L}$
and $\bm{T}:\mathcal{X}\mapsto\mathbb{R}^{L\times M}$ are componentwise
closed, proper, convex vector- and matrix-valued functions, respectively. We
allow uncertainty to affect only the affine constraints of the problem and,
similarly, the first- and second-stage decisions to interact only via the
affine part.

Since the true underlying distribution $\mathbb{P}$ is unknown, this two-stage
optimization formulation is ill-posed. Nevertheless, $\mathbb{P}$ is typically
observable through a finite amount of historical data. We assume that we have
access to $N$ such independent and identically distributed observations, which
we denote by $\\{\hat{\bm{\xi}}^{(1)},\ldots,\hat{\bm{\xi}}^{(N)}\\}$. We also
assume that generating additional data (e.g., via Monte Carlo computer
simulations) is either costly or impossible. Thus, it is imperative to use the
given data most efficiently.

A popular approach for solving the two-stage problem in such cases is to use
_sample average approximation_ [59]. This approach replaces the true
distribution with the empirical distribution
$\hat{\mathbb{P}}_{N}=\frac{1}{N}\sum_{i=1}^{N}\delta_{\hat{\bm{\xi}}^{(i)}}$,
where $\delta_{\hat{\bm{\xi}}^{(i)}}$ denotes the Dirac distribution at
$\hat{\bm{\xi}}^{(i)}$. In the context of rare events, however, obtaining
accurate estimates of the true distribution and, hence, the optimal solution
of the true two-stage problem may require unrealistically large amounts of
data. The following example illustrates this phenomenon even if there are only
$M=3$ uncertain parameters.

123$0.25\%$4$1\%$5$0.25\%$67100100100100$1$$1-\epsilon$$1-\epsilon$$1$$1$$1-\epsilon$$1-\epsilon$$1$
Figure 1: Illustrative network flow instance.

######  Example 1.

Consider the network shown in Figure 1. Let $A$ denote the set of arcs. The
goal is to decide arc capacities $\bm{x}\in\mathbb{R}^{\lvert A\rvert}_{+}$ so
as to route flow originating from the left layer of nodes to satisfy demand in
the right layer of nodes. The per-unit cost $\bm{c}\in\mathbb{R}^{\lvert
A\rvert}$ of installing capacity on each arc is denoted above the arc, and the
first-stage deterministic cost is given by $\bm{c}^{\top}\bm{x}$. The middle
layer of nodes can fail independently of each other, and the numbers above the
nodes denote their failure probabilities. If node $i$ fails (indicated by
${\xi}_{i}=1$), then all arcs incident to that node become unusable, and any
resulting supply shortfall is penalized at a cost of $1,000$ per unit. For a
given realization $\bm{\xi}\in\\{0,1\\}^{3}$ of the node failures, the loss
function $\mathcal{Q}(\bm{x},\bm{\xi})$ is simply the total penalty cost under
that realization and can be modeled as the optimal objective value of a linear
program.

We can show that the optimal solution assigns a capacity of 100 units to every
arc (see Appendix A). In practice, however, the true failure probabilities are
unknown. We therefore estimate the distribution using a sample average
approximation. Figure 2 shows the performance of the solutions obtained by
replacing the true distribution $\mathbb{P}$ with the empirical distribution
$\hat{\mathbb{P}}_{N}$ resulting from different sample sizes $N$. In
particular, the figure shows the difference between the total expected costs
(computed under the true distribution) of the sample average solution and the
true optimal solution. We observe that despite the small dimensionality
$(M=3)$, $N\geq 1000$ samples are required for the sample average solution to
estimate the true solution to an accuracy of 5%.

1234567100100100100$100$$100$$100$$100$ (a) The optimal arc capacities
determined by the sample average approximation when the empirical distribution
$\hat{\mathbb{P}}_{N}$ puts all weight on the realization $\bm{\xi}=\bm{0}$,
where no nodes fail (see Appendix A).

(b) The optimality gap of the sample average approximation with respect to
true optimal objective for increasing sample size. The mean (solid line) and
standard deviation (shaded) are estimated by using 1,000 statistically
independent sets of samples.

Figure 2: Performance of the sample average approximation on the network flow
instance from Figure 1.

##  2 Distributionally Robust Approach for Discrete Rare Events

The high dimensionality and rare occurrence of failure states render accurate
estimation of the underlying distribution difficult. To remedy this situation,
we adopt a _distributionally robust_ approach, and construct an _ambiguity
set_ $\mathcal{P}$ of possible distributions that are consistent with the
observed data. We then minimize the _worst-case_ expected costs over all
distributions in the ambiguity set. Specifically, we consider distributionally
robust two-stage conic optimization problems, of the form

$\mathop{\text{minimize}}_{\bm{x}\in\mathcal{X}}\;c(\bm{x})+\sup_{\mathbb{P}\in\mathcal{P}}\mathbb{E}_{\mathbb{P}}\left[\mathcal{Q}(\bm{x},\tilde{\bm{\xi}})\right].$
(2)

The ambiguity set $\mathcal{P}$ must be chosen such that it contains the true
distribution with high confidence or, at the very least, distributions that
assign nonzero probability to the rare events. We focus on the Wasserstein
ambiguity set, defined as the set of distributions that are close to the
empirical distribution $\hat{\mathbb{P}}_{N}$ with respect to the Wasserstein
distance $d_{W}$:

$\mathcal{P}=\left\\{\mathbb{Q}\in\mathcal{M}(\Xi):d_{W}(\mathbb{Q},\hat{\mathbb{P}}_{N})\leq\varepsilon\right\\}.$
(3)

Here, $\mathcal{M}(\Xi)$ denotes the set of all distributions supported on
$\Xi$, and $\varepsilon\geq 0$ is the radius of the Wasserstein ball. Given
any underlying metric $d(\cdot,\cdot)$ on the support set $\Xi$, the
Wasserstein distance $d_{W}(\mathbb{P},\mathbb{P}^{\prime})$ between two
distributions $\mathbb{P},\mathbb{P}^{\prime}$ can be defined as follows:

$d_{W}(\mathbb{P},\mathbb{P}^{\prime})=\min_{\Pi\in\mathcal{M}(\Xi\times\Xi)}\left\\{\sum_{\bm{\xi}\in\Xi}\sum_{\bm{\xi}^{\prime}\in\Xi}d(\bm{\xi},\bm{\xi}^{\prime})\Pi(\bm{\xi},\bm{\xi}^{\prime}):\begin{array}[]{l}\Pi\text{
is a coupling of $\mathbb{P}$ and $\mathbb{P}^{\prime}$}\end{array}\right\\},$

We motivate the choice of Wasserstein ambiguity sets in Section 2.1. The
crucial role of the metric $d(\cdot,\cdot)$ and the radius $\varepsilon$ of
the Wasserstein ball $\mathcal{P}$ is discussed in Section 2.2

###  2.1 Advantages of Wasserstein ambiguity sets for discrete rare events

One can construct several ambiguity sets of distributions that are consistent
with observed data. Broadly, these are either sets of distributions that
satisfy constraints on their moments [25, 66] or those defined as balls
centered on a reference distribution with respect to a metric such as the
$\phi$-divergence [9, 7] or the Wasserstein distance [67, 50, 30, 72]. Note
that the high dimensionality and sparsity of the training data in the case of
rare events prevent us from obtaining reliable estimates of moments, ruling
out moment-based sets. In contrast, metric-based sets have the attractive
feature of tying directly with available data; indeed, the empirical
distribution corresponding to the training data is a natural choice for the
reference distribution.

Ambiguity sets based on $\phi$-divergence, especially Kullback-Liebler
divergence, have certain shortcomings that are not shared by their Wasserstein
counterparts. First, the former can exclude the true distribution while
including pathological distributions; in fact, they can fail to represent
confidence sets for the unknown true distribution [30]. In contrast, the
latter contains the true distribution with high confidence for an appropriate
choice of $\varepsilon$, and hence the optimal value of the corresponding
distributionally robust problem provides an upper confidence bound on the true
out-of-sample cost [50]. Second, Kullback-Liebler and other $\phi$-divergence
based sets contain only those distributions that are absolutely continuous
with respect to reference distribution; that is, these distributions can
assign positive probability only to those realizations for which the empirical
distribution also assigns positive mass [7]. This situation is undesirable for
applications where uncertain events are rare, since the majority of possible
uncertain states are unlikely to have been observed empirically. In contrast,
the Wasserstein ball of any positive radius $\varepsilon>0$ and corresponding
to any finite-valued metric $d(\cdot,\cdot)$ includes distributions that
assign nonzero probability to any arbitrary realization. Indeed, this property
ensures that solutions of the distributionally robust problem (2) are robust
to the occurrence of rare events. We illustrate this via Example 1.

######  Example 1 (continued).

Consider again the network in Figure 1. In addition to the sample average
approximation, Figure 3(b) now compares the performance of solutions computed
by using the distributionally robust formulation (2) with a Wasserstein
ambiguity set $\mathcal{P}$ of radius $\varepsilon=10^{-3}$, centered around
the empirical distribution $\hat{\mathbb{P}}_{N}$ resulting from different
sample sizes $N$. Here, we use the metric
$d(\bm{\xi},\bm{\xi}^{\prime})=\lVert\bm{\xi}-\bm{\xi}^{\prime}\rVert_{1}$
induced by the 1-norm. The figure shows the difference in total expected costs
(computed under the true distribution) of the distributionally robust and
sample average solutions with respect to the true optimal solution. We observe
that the distributionally robust solution strongly outperforms the sample
average solution while being more stable to changes in the training data
(i.e., its performance has a smaller variance for a fixed $N$). The
performance of the solutions computed using a $\phi$-divergence ambiguity set
are presented in the next subsection.

1234567100100100100$\frac{2}{3}100$$\frac{1}{3}100$$\frac{1}{3}100$$\frac{2}{3}100$$\frac{2}{3}100$$\frac{1}{3}100$$\frac{1}{3}100$$\frac{2}{3}100$
(a) The optimal arc capacities determined by the distributionally robust
formulation when the empirical distribution $\hat{\mathbb{P}}_{N}$ puts all
its weight on $\bm{\xi}=\bm{0}$, where none of the nodes fail (see Appendix
A).

(b) Comparison of the optimality gap of the distributionally robust
formulation (dashed orange) with the sample average approximation (solid blue)
for increasing sample size.

Figure 3: Performance of the distributionally robust formulation with a
Wasserstein ambiguity set of radius $\varepsilon=10^{-3}$ on the network flow
instance from Figure 1.

###  2.2 Choice of the underlying metric and radius of the Wasserstein ball

The underlying metric $d(\cdot,\cdot)$ should ideally have the following
properties: (i) if $d(\bm{\xi}^{\prime},\bm{\xi}^{\prime\prime})$ is small,
then the probabilities of occurrence of the two realizations
$\bm{\xi}^{\prime},\bm{\xi}^{\prime\prime}\in\Xi$ should be similar; and (ii)
if $\bm{\xi}^{\prime\prime}$ is rarer than $\bm{\xi}^{\prime}$, then for some
fixed (say nominal) realization $\bm{\xi}\in\Xi$ (e.g., in a network, this
could be the realization where none of the elements fail), we must have
$d(\bm{\xi},\bm{\xi}^{\prime})\leq d(\bm{\xi},\bm{\xi}^{\prime\prime})$. These
properties can be satisfied whenever
$d(\bm{\xi},\bm{\xi^{\prime}})=\lVert\bm{\xi}-\bm{\xi^{\prime}}\rVert$ is
induced by a norm on $\mathbb{R}^{M}$ and by using an appropriate bit
representation of the sample space. Throughout the paper, we will therefore
assume that the metric $d(\cdot,\cdot)$ is induced by an arbitrary norm,
although our results also apply when $d(\cdot,\cdot)$ is any mixed-integer
conic-programming-representable metric.

A noteworthy example of a metric that does not satisfy the above requirement
is the _discrete metric_ , defined as
$d(\bm{\xi}^{\prime},\bm{\xi}^{\prime\prime})=1$ whenever
$\bm{\xi}^{\prime}\neq\bm{\xi}^{\prime\prime}$ and $0$ otherwise. In this
case, the Wasserstein distance is equivalent to the _total variation distance_
, and the distributions that solve the inner supremum in (2) can assign
positive probability only to the training samples and the worst-case
realization [36, 56]. In a network with rare failures, this means that only
past observed realizations and the realization corresponding to complete
network failure are taken into account, resulting in poor out-of-sample
performance.

######  Example 1 (continued).

Suppose that we define the metric $d(\cdot,\cdot)$ to be the discrete metric.
Then, the performance of solutions of the distributionally robust problem (2)
with a total variation ambiguity set $\mathcal{P}$ of _any_ radius
$\varepsilon\geq 0$ is equal to or worse than that of the sample average
solution.

For a given choice of the underlying metric, the radius $\varepsilon$ of the
ambiguity set allows us to control the level of conservatism of solutions of
(2). Specifically, given a confidence level $\beta\in(0,1)$, one can choose
the radius as a function of $\beta$ and the number of observations $N$ such
that the true distribution $\mathbb{P}$ is contained in the ambiguity set with
high probability:

$\mathbb{P}^{N}\left[d_{W}(\mathbb{P},\hat{\mathbb{P}}_{N})\leq\varepsilon_{N}(\beta)\right]\geq
1-\beta.$ (4)

Here, $\mathbb{P}^{N}$ is the product distribution that governs the observed
data $\\{\hat{\bm{\xi}}^{(1)},\ldots,\hat{\bm{\xi}}^{(N)}\\}$. It was shown in
[50] that (4) holds if we select
$\varepsilon_{N}(\beta)=c_{0}\left({N^{-1}\log\beta^{-1}}\right)^{1/M}$, where
$c_{0}$ is a problem-dependent constant. Since this choice can lead to
unnecessarily large values for the radius, the authors suggest solving the
two-stage problem (2) for several a priori fixed choices of $\varepsilon$ and
then using cross-validation (e.g., $k$-fold cross-validation) to pick a good
value. In the following, we show that Sanov’s theorem can be used to obtain
stronger finite sample guarantees; that is, less conservative values for
$\varepsilon_{N}(\beta)$, by exploiting the finiteness of the support $\Xi$.

######  Theorem 1 (Finite sample guarantee).

For every fixed sample size $N>0$, and confidence level $\beta\in(0,1)$, the
probabilistic guarantee (4) holds whenever

$\varepsilon_{N}(\beta)\geq
D\sqrt{(2N)^{-1}\left(\lvert\Xi\rvert\log(N+1)+\log\beta^{-1}\right)},$ (5)

where
$D\coloneqq\max_{\bm{\xi},\bm{\xi}^{\prime}\in\Xi}d(\bm{\xi},\bm{\xi}^{\prime})$
is the diameter of $\Xi$ with respect to the metric $d$.

###### Proof.

See Appendix G. ∎

The right-hand side of (5) indicates that for a fixed choice of the support
$\Xi$ and confidence level $\beta$, $\varepsilon_{N}(\beta)$ is roughly
proportional to $\sqrt{N^{-1}\log(N+1)}$. For such choices, the optimal value
of the corresponding distributionally robust two-stage problem (2) can be
expected to provide an upper bound on the true (unknown) out-of-sample cost.
We empirically verify this upper bound in Section 5, when we vary
$\varepsilon=\nu\sqrt{N^{-1}\log(N+1)}$ as a function of a scalar $\nu$.
Theorem 1 also indicates that $\varepsilon_{N}(\beta)\to 0$ as the sample size
becomes large, $N\to\infty$, and that $\varepsilon_{N}(\beta)\to\infty$ if we
are overly conservative, $\beta\to 0$. This observation also elucidates that
the distributionally robust formulation (2) generalizes both classical
stochastic and robust optimization.

######  Remark 1 (Reduction to two-stage stochastic and robust optimization).

The distributionally robust two-stage problem (2) reduces to the classical
sample average approximation whenever the radius of the ambiguity set
$\varepsilon=0$, since $\mathcal{P}=\\{\hat{\mathbb{P}}_{N}\\}$ reduces to a
singleton in this case. Similarly, it reduces to a classical two-stage robust
optimization problem whenever
$\varepsilon\geq\max_{\bm{\xi},\bm{\xi}^{\prime}\in\Xi}d(\bm{\xi},\bm{\xi}^{\prime})$
since $\Xi$ is compact and $\mathcal{P}$ contains all Dirac distributions,
$\delta_{\bm{\xi}}$, $\bm{\xi}\in\Xi$, in this case. Therefore, the worst-case
expectation in (2) reduces to
$\max_{\bm{\xi}\in\Xi}\mathcal{Q}(\bm{x},\bm{\xi})$.

###  2.3 Contributions

This paper addresses the relatively unexplored topic of rare high-impact
uncertainties through the lens of data-driven distributionally robust
optimization. Existing methods for addressing rare high-impact uncertainties
[6, 17] are few, and they are all based on variants of Monte Carlo methods
(e.g., importance sampling), which require the existence of a probability
distribution that can be sampled to generate additional observations.

In contrast, our techniques fall within the scope of distributionally robust
optimization [58, 57]. A surge in the popularity of Wasserstein ambiguity sets
has coccurred in this area because of recent results [30, 50, 11, 14, 40] that
have established (i) strong finite sample and asymptotic guarantees of their
formulations, (ii) connections with function regularization in machine
learning, and (iii) tractable reformulations for various classes of loss
functions $\mathcal{Q}(\bm{x},\bm{\xi})$ and metric spaces $(\Xi,d)$. Most of
the tractability results are for one-stage problems, however, with piecewise
linear loss functions and continuous random parameters. Existing tractability
results for two-stage problems are limited to the case where
$\mathcal{Q}(\bm{x},\bm{\xi})$ is the optimal value of a linear program, $\Xi$
is a polyhedron, and $d$ is induced by the $1$-norm, see [50, 35]. Sufficient
conditions for zero-one supports $\Xi$ and type-$\infty$ Wasserstein ambiguity
sets have been established in [68]. In the absence of sufficient conditions
that ensure tractability, one resorts either to iterative global optimization
methods (e.g., see [67, 71, 46]) or tractable approximations. The latter
commonly include discretization schemes (of which sample average approximation
is a special case) [45, 20] and decision rule methods [8, 32, 12].

In this paper, we extend the state of the art in data-driven optimization by
studying two-stage conic programs with a particular focus on high-dimensional
zero-one uncertainties. This is crucial because existing reformulations and
algorithms for distributionally robust optimization with finitely supported
distributions [53, 9, 11, 5] scale with the size of the support set
$\lvert\Xi\rvert$, which can grow exponentially large in such cases. We
circumvent this exponential growth by utilizing tractable conservative
approximations inspired by lift-and-project convexification techniques in
global optimization [43, 60, 41].

Closest in spirit to our work are the papers of [35, 69, 2] who consider the
case where the second-stage loss $\mathcal{Q}(\bm{x},\bm{\xi})$ is the optimal
value of a linear program with uncertain right-hand sides and the support set
$\Xi$ is a polytope. In this setting, [35, 69] reformulate (2) as a copositive
cone program and then approximate this using semidefinite programming, whereas
[2] provide approximations by leveraging reformulation-linearization
techniques from bilinear programming. Although some extensions of these
approaches to the case of zero-one support sets $\Xi$ have been made [49, 37],
problems where $\mathcal{Q}(\bm{x},\bm{\xi})$ is the optimal value of a conic
program have not been addressed. This is partly because such extensions lead
to so-called _generalized copsitive programs_ or _set-semidefinite programs_
(e.g., see [18, 28]), and relatively little is known about their tractable
approximations. In contrast, generalizations of lift-and-project techniques to
mixed zero-one conic problems are fairly well known (e.g., see [61, 21]), and
we exploit these to derive tractable approximations for distributionally
robust optimization. The relationship between convexification hierarchies
based on copositive programming and lift-and-project techniques (for specific
problem classes) has been explored in [54, 19].

We highlight the following main contributions:

  1. 1.

By exploiting ideas from penalty methods and bilinear programming, we develop
reformulations of the Wasserstein distributionally robust two-stage problem
(2) that reduce its solution to optimization problems over the convex hulls of
mixed-integer conic representable sets. Extensions to conditional value-at-
risk are also presented.

  2. 2.

We provide sufficient conditions for our convex hull reformulations and hence
the distributionally robust two-stage problem (2) to be tractable. We also
show that they are generically NP-hard, however, even if there are no first-
stage decisions and the second-stage loss function is the optimal value of a
two-dimensional linear program with uncertain objective coefficients.

  3. 3.

By using lift-and-project hierarchies to approximate the convex hull of the
mixed-integer conic representable sets, we derive tractable conservative
approximations of the distributionally robust two-stage problem (2), and we
provide practical guidelines to compute them efficiently. The approximations
are tractable irrespective of the aforementioned sufficient conditions, and
they become exact as the Wasserstein radius $\varepsilon$ shrinks to zero.

  4. 4.

We demonstrate the practical viability of our method and its out-of-sample
performance on challenging nonlinear optimal power flow and multi-commodity
network design problems that are affected by rare network contingencies, and
we study its behavior as a function of the rarity and impact of these
contingencies, illustrating improvements over classical sample average and
two-stage robust optimization formulations.

The rest of the paper is organized as follows. Section 3 derives the mixed-
integer conic programming representation of interest, Section 4 derives their
lift-and-project approximations, and Section 5 reports numerical results. For
ease of exposition, the complexity analysis as well as proofs of all
assertions is deferred to the appendix.

##### Notation.

Vectors and matrices are printed in bold lowercase and bold uppercase letters,
respectively, while scalars are printed in regular font. The set of non-
negative integers and reals is denoted by $\mathbb{Z}_{+}$ and
$\mathbb{R}_{+}$, respectively. For any positive integer $N$, we define $[N]$
as the index set $\left\\{1,\ldots,N\right\\}$. We use $\bm{\mathrm{e}}_{k}$
to denote the $k^{\text{th}}$ unit basis vector, $\bm{\mathrm{e}}$ to denote
the vector of ones, $\mathbf{I}$ to denote the identity matrix, and $\bm{0}$
to denote the vector or matrix of zeros, respectively; their dimensions will
be clear from the context. For a matrix $\bm{A}$, we use
$\mathop{\text{vec}}(\bm{A})$ to denote the vector obtained by stacking the
columns of $\bm{A}$ in order. The inner product between two matrices
$\bm{A},\bm{B}\in\mathbb{R}^{m\times n}$ is denoted by
$\langle\bm{A},\bm{B}\rangle\coloneqq\sum_{i\in[m]}\sum_{j\in[n]}A_{ij}B_{ij}$.
We use
$\mathcal{C}^{n}=\left\\{(\bm{x},t)\in\mathbb{R}^{n-1}\times\mathbb{R}:\lVert\bm{x}\rVert\leq
t\right\\}$ to denote the norm cone associated with the norm
$\lVert\cdot\rVert$. For a logical expression $\mathcal{E}$, we define
$\mathbb{I}[\mathcal{E}]$ as the indicator function which takes a value of $1$
if $\mathcal{E}$ is true and $0$ otherwise. Throughout the paper, we refer to
an optimization problem as _tractable_ if it can be solved in polynomial time
in the size of its input data, and _intractable_ if it is NP-hard.

##  3 Mixed-Integer Conic Representations

Throughout the paper, we assume that the distributionally robust two-stage
problem (2) satisfies the assumptions of _complete_ and _sufficiently
expensive recourse_.

  * (A1)

For every realization $\bm{\xi}\in\Xi$, there exists
$\bm{y}^{+}\in\mathrm{int}(\mathcal{Y})$ such that
$\bm{W}(\bm{\xi})\bm{y}^{+}>0$.

  * (A2)

For every first-stage decision $\bm{x}\in\mathcal{X}$ and every realization
$\bm{\xi}\in\Xi$, the second-stage loss function
$\mathcal{Q}(\bm{x},\bm{\xi})$ is bounded.

A natural way to ensure these assumptions is to add slack variables in the
formulation of the second-stage problem $\mathcal{Q}(\bm{x},\bm{\xi})$ and
penalize them in the objective function. Whenever the assumptions are
satisfied, they imply that _(i)_ $\mathcal{Q}(\bm{x},\bm{\xi})$ is always
strictly feasible and bounded, _(ii)_ the dual of
$\mathcal{Q}(\bm{x},\bm{\xi})$, given in the following, is always feasible,
and _(iii)_ strong conic duality holds between the second-stage problem and
its dual, $\mathcal{Q}(\bm{x},\bm{\xi})=\mathcal{Q}_{d}(\bm{x},\bm{\xi})$,
where

$\mathcal{Q}_{d}(\bm{x},\bm{\xi}):=\sup_{\bm{\lambda}\in\mathbb{R}^{L}_{+}}\left\\{\left[\bm{T}(\bm{x})\bm{\xi}+\bm{h}(\bm{x})\right]^{\top}\bm{\lambda}:\bm{q}(\bm{\xi})-\bm{W}(\bm{\xi})^{\top}\bm{\lambda}\in\mathcal{Y}^{*}\right\\}.$
(6)

Here, $\mathcal{Y}^{*}$ denotes the dual cone of $\mathcal{Y}$. We assume that
the uncertain vectors and matrices in (1) are affine and can be represented as
$\bm{q}(\bm{\xi})=\bm{q}_{0}+\bm{Q}\bm{\xi}$ and
$\bm{W}(\bm{\xi})=\bm{W}_{0}+\sum_{j\in[M]}\xi_{j}\bm{W}_{j},$ where
$\bm{Q}\in\mathbb{R}^{N_{2}\times M}$, and for each $j\in\\{0,1,\ldots,M\\}$,
we have $\bm{W}_{j}\in\mathbb{R}^{L\times N_{2}}$.

A consequence of the above assumptions is the following lemma, which states
that computing the worst-case expectation in the two-stage problem (2) is
equivalent to averaging $N$ worst-case values of the loss function
$\mathcal{Q}(\bm{x},\bm{\xi})$ over $\bm{\xi}\in\Xi$, each regularized by one
of the training samples. We omit the proof since it follows directly from the
compactness of $\Xi$ and the definition of the Wasserstein ambiguity set
$\mathcal{P}$ in (3); see [30, 15] for proofs in much more general settings.

######  Lemma 1.

The distributionally robust two-stage problem (2) admits the following
reformulation:

$\displaystyle\mathop{\text{minimize}}_{\bm{x}\in\mathcal{X},\alpha\geq
0}\;\displaystyle
c(\bm{x})+\alpha\varepsilon+\frac{1}{N}\sum_{i=1}^{N}\max_{\bm{\xi}\in\Xi}\left\\{\mathcal{Q}(\bm{x},\bm{\xi})-\alpha
d(\bm{\xi},\hat{\bm{\xi}}^{(i)})\right\\}.$ (7)

The remainder of this section establishes that the inner maximization in (7)
is equivalent to optimizing a linear function over the convex hull of the
feasible region of a _mixed-integer conic program_ (MICP). The following
theorem is key to establishing this result.

######  Theorem 2 (Convex hull reformulation).

The distributionally robust two-stage problem (2) admits the following convex
hull reformulation:

$\displaystyle\mathop{\text{minimize}}_{\bm{x}\in\mathcal{X},\alpha\geq
0}\;\;\displaystyle
c(\bm{x})+\alpha\varepsilon+\frac{1}{N}\sum_{i=1}^{N}Z_{i}(\bm{x},\alpha),$
(8)

where, for each $i\in[N]$, we define the function
$Z_{i}:\mathcal{X}\times\mathbb{R}_{+}\mapsto\mathbb{R}$ and the set
$\mathcal{Z}_{i}$ as follows:

$\displaystyle\displaystyle
Z_{i}(\bm{x},\alpha)=\mathop{\text{maximize}}_{(\bm{\xi},\bm{\lambda},\bm{\Lambda},\tau)\in\mathop{\mathrm{cl\,conv}}(\mathcal{Z}_{i})}\displaystyle\left\\{\langle\bm{T}(\bm{x}),\bm{\Lambda}\rangle+\bm{h}(\bm{x})^{\top}\bm{\lambda}-\alpha\tau\right\\}$
(9a)
$\displaystyle\displaystyle\mathcal{Z}_{i}=\left\\{(\bm{\xi},\bm{\lambda},\bm{\Lambda},\tau)\in\Xi\times\mathbb{R}^{L}_{+}\times\mathbb{R}^{L\times
M}\times\mathbb{R}_{+}:\begin{array}[]{l}\displaystyle\bm{\Lambda}=\bm{\lambda}\bm{\xi}^{\top},\;(\bm{\xi}-\hat{\bm{\xi}}^{(i)},\tau)\in\mathcal{C}^{M+1}\\\
\displaystyle\bm{q}_{0}+\bm{Q}\bm{\xi}-\bm{W}_{0}^{\top}\bm{\lambda}-\sum_{j\in[M]}\bm{W}_{j}^{\top}\bm{\Lambda}\bm{\mathrm{e}}_{j}\in\mathcal{Y}^{*}\end{array}\right\\}.$
(9d)

###### Proof.

See Appendix G. ∎

The inner optimization problem (9a) is over the closed convex hull of the set
$\mathcal{Z}_{i}$, which couples the binary uncertain parameters $\bm{\xi}$
with the continuous dual variables $\bm{\lambda}$ via the bilinear equation
$\bm{\Lambda}=\bm{\lambda}\bm{\xi}^{\top}$. This set is, therefore, not the
feasible region of an MICP. We propose two approaches to ensure MICP
representability. The first is to linearize the bilinear equation
$\bm{\Lambda}=\bm{\lambda}\bm{\xi}^{\top}$ using McCormick inequalities, which
requires _a priori_ upper bounds on the dual variables $\bm{\lambda}$; we
briefly discuss this in Section 3.1. The second is to reformulate the loss
function $\mathcal{Q}(\bm{x},\bm{\xi})$ using ideas from penalty methods that
circumvents any bilinear terms; this approach is the subject of Section 3.2.
We compare the two approaches and summarize their merits in Section 3.3. We
note that our results also apply if the risk-neutral expectation in the
objective function of the two-stage problem (2) is replaced with the
conditional value-at-risk; for ease of exposition, we defer this analysis to
the appendix.

###  3.1 Linearized reformulation

The decision variables $\bm{\lambda}$ of the dual problem
$\mathcal{Q}_{d}(\bm{x},\bm{\xi})$ must be necessarily bounded for any fixed
value of $\bm{x}\in\mathcal{X}$ and $\bm{\xi}\in\Xi$. Indeed, under
assumptions (A1) and (A2), the value of $\mathcal{Q}_{d}(\bm{x},\bm{\xi})$ is
bounded for any fixed $\bm{x}\in\mathcal{X}$ and $\bm{\xi}\in\Xi$; and since
$\mathcal{X}$ and $\Xi$ are compact sets, the variables $\bm{\lambda}$ must
also be necessarily bounded. Suppose that
$\bar{\bm{\lambda}}\in\mathbb{R}^{L}_{+}$ are _a priori_ known upper bounds
(independent of $\bm{x}$ and $\bm{\xi}$) on these variables. Such bounds may
be analytically known whenever we explicitly add slack variables to ensure
feasibility of the second-stage problem or if the latter has some structure
(e.g., see [29]). Whenever such bounds are known, we can exactly linearize the
bilinear equation $\bm{\Lambda}=\bm{\lambda}\bm{\xi}^{\top}$ using McCormick
inequalities since $\bm{\xi}$ is binary-valued (e.g., see [33]), and
reformulate the set $\mathcal{Z}_{i}$ in (9d) as the feasible region of an
MICP:

$\displaystyle\mathcal{Z}_{i}=\left\\{(\bm{\xi},\bm{\lambda},\bm{\Lambda},\tau)\in\Xi\times\mathbb{R}^{L}_{+}\times\mathbb{R}^{L\times
M}_{+}\times\mathbb{R}_{+}:\begin{array}[]{l}\displaystyle(\bm{\xi}-\hat{\bm{\xi}}^{(i)},\tau)\in\mathcal{C}^{M+1}\\\
\displaystyle\bm{\Lambda}-\bm{\lambda}\bm{\mathrm{e}}^{\top}+\bar{\bm{\lambda}}(\bm{\mathrm{e}}-\bm{\xi})^{\top}\in\mathbb{R}^{L\times
M}_{+}\\\
\displaystyle\bm{\lambda}\bm{\mathrm{e}}^{\top}-\bm{\Lambda}\in\mathbb{R}^{L\times
M}_{+},\;\bar{\bm{\lambda}}\bm{\xi}^{\top}-\bm{\Lambda}\in\mathbb{R}^{L\times
M}_{+}\\\
\displaystyle\bm{q}_{0}+\bm{Q}\bm{\xi}-\bm{W}_{0}^{\top}\bm{\lambda}-\sum_{j\in[M]}\bm{W}_{j}^{\top}\bm{\Lambda}\bm{\mathrm{e}}_{j}\in\mathcal{Y}^{*}\end{array}\right\\}.$
(9d-$\ell$)

The MICP representation (9d-$\ell$) adds $O(ML)$ variables and constraints for
each $\mathcal{Z}_{i}$, $i\in[N]$, which can be prohibitively large. Moreover,
analytical upper bounds $\bar{\bm{\lambda}}\in\mathbb{R}^{L}_{+}$ on the dual
variables, which are independent of (and valid for all) $\bm{x}$ and
$\bm{\xi}$, may be unavailable. The following section shows that we can ensure
MICP representability without _a priori_ knowledge of these bounds, and at the
expense of adding far fewer variables and constraints.

###  3.2 Penalty reformulation

Our goal in this section will be to move the uncertainty to the objective
function of the second-stage problem, $\mathcal{Q}(\bm{x},\bm{\xi})$. This
approach is motivated by the following corollary to Theorem 2 that shows that
MICP-representability is guaranteed if only the objective function of
$\mathcal{Q}(\bm{x},\bm{\xi})$ is uncertain.

######  Corollary 1 (Convex hull reformulation for objective uncertainty).

Suppose that only the objective function of the second-stage problem
$\mathcal{Q}(\bm{x},\bm{\xi})$ is uncertain: $\bm{W}(\bm{\xi})=\bm{W}_{0}$ and
$\bm{T}(\bm{x})=\bm{0}$. Then the distributionally robust two-stage problem
(2) admits reformulation (8) where, for each $i\in[N]$, we define the function
$Z_{i}:\mathcal{X}\times\mathbb{R}_{+}\mapsto\mathbb{R}$ and the MICP-
representable set $\mathcal{Z}_{i}$ as follows:

$\displaystyle\displaystyle
Z_{i}(\bm{x},\alpha)=\mathop{\text{maximize}}_{(\bm{\xi},\bm{\lambda},\tau)\in\mathop{\mathrm{cl\,conv}}(\mathcal{Z}_{i})}\displaystyle\left\\{\bm{h}(\bm{x})^{\top}\bm{\lambda}-\alpha\tau\right\\}$
(10a)
$\displaystyle\displaystyle\mathcal{Z}_{i}=\left\\{(\bm{\xi},\bm{\lambda},\tau)\in\Xi\times\mathbb{R}^{L}_{+}\times\mathbb{R}_{+}:\displaystyle(\bm{\xi}-\hat{\bm{\xi}}^{(i)},\tau)\in\mathcal{C}^{M+1},\;\;\displaystyle\bm{q}_{0}+\bm{Q}\bm{\xi}-\bm{W}_{0}^{\top}\bm{\lambda}\in\mathcal{Y}^{*}\right\\}.$
(10b)

For the remainder of the paper, we shall make the following additional
assumption of _fixed recourse_ which will be key to achieving our goal.

  * (A3)

For every realization $\bm{\xi}\in\Xi$ and every first-stage decision
$\bm{x}\in\mathcal{X}$, we have $\bm{W}(\bm{\xi})=\bm{W}_{0}$ and
$\bm{T}(\bm{x})=\bm{T}_{0}$, respectively.

We note that since $\bm{\xi}$ is binary valued, this assumption is without
loss of generality if the primal decision variables $\bm{x}$ and $\bm{y}$ are
bounded. Indeed, in such cases, we can exactly linearize any products of
uncertain parameters and decisions in the second-stage problem
$\mathcal{Q}(\bm{x},\bm{\xi})$ using McCormick inequalities [33], to ensure
that it satisfies the assumption of fixed recourse. Under this additional
assumption, the following theorem states that the second-stage problem
$\mathcal{Q}(\bm{x},\bm{\xi})$ with constraint uncertainty can be equivalently
reformulated as one with objective uncertainty.

######  Theorem 3 (Penalty reformulation of the loss function).

There exists a sufficiently large, yet _finite_ , penalty parameter $\rho>0$
such that the second-stage loss function $\mathcal{Q}(\bm{x},\bm{\xi})$ in (1)
is equivalent to

$\mathcal{Q}^{\rho}(\bm{x},\bm{\xi}):=\inf_{\bm{y}\in\mathcal{Y},\bm{z}\in[0,1]^{M}}\left\\{\bm{q}(\bm{\xi})^{\top}\bm{y}+\rho\left((\bm{\mathrm{e}}-2\bm{\xi})^{\top}\bm{z}+\bm{\mathrm{e}}^{\top}\bm{\xi}\right):\bm{W}_{0}\bm{y}\geq\bm{T}_{0}\bm{z}+\bm{h}(\bm{x})\right\\}.$
(11)

###### Proof.

See Appendix G. ∎

In conjunction with Corollary 1, Theorem 3 implies that the distributionally
robust two-stage problem (2) admits a convex hull reformulation of the form
(8), where the function
$Z_{i}:\mathcal{X}\times\mathbb{R}_{+}\mapsto\mathbb{R}$ and the MICP-
representable set $\mathcal{Z}_{i}$ for each $i\in[N]$ are given as follows:

$\displaystyle\displaystyle
Z_{i}(\bm{x},\alpha)=\mathop{\text{maximize}}_{(\bm{\xi},\bm{\lambda},\bm{\mu},\tau)\in\mathop{\mathrm{cl\,conv}}(\mathcal{Z}_{i})}\displaystyle\left\\{\bm{h}(\bm{x})^{\top}\bm{\lambda}+\rho(\bm{\mathrm{e}}^{\top}\bm{\xi})-\bm{\mathrm{e}}^{\top}\bm{\mu}-\alpha\tau\right\\}$
(9a-$\rho$)
$\displaystyle\displaystyle\mathcal{Z}_{i}=\left\\{(\bm{\xi},\bm{\lambda},\bm{\mu},\tau)\in\Xi\times\mathbb{R}^{L}_{+}\times\mathbb{R}^{M}_{+}\times\mathbb{R}_{+}:\begin{array}[]{l}\displaystyle(\bm{\xi}-\hat{\bm{\xi}}^{(i)},\tau)\in\mathcal{C}^{M+1}\\\
\displaystyle\bm{q}_{0}+\bm{Q}\bm{\xi}-\bm{W}_{0}^{\top}\bm{\lambda}\in\mathcal{Y}^{*}\\\
\displaystyle\rho(\bm{\mathrm{e}}-2\bm{\xi})+\bm{T}_{0}^{\top}\bm{\lambda}+\bm{\mu}\in\mathbb{R}^{M}_{+}\end{array}\right\\}.$
(9d-$\rho$)

In contrast to the linearized reformulation (9d-$\ell$), the MICP
representation (9d-$\rho$) adds only $O(M+L)$ variables and constraints for
each $\mathcal{Z}_{i}$, $i\in[N]$.

The reformulation (11) is reminiscent of penalty methods in nonlinear
programming. However, in contrast to the latter, which suffer from numerical
issues because the penalty parameter $\rho$ must be driven to $\infty$, a
finite value for $\rho$ can be precomputed in our case. This process only
requires solving the classical robust optimization formulation over any
support $\Xi^{0}\supseteq\Xi$; for example, some choices are $\Xi^{0}=\Xi$ or
$\Xi^{0}=\\{0,1\\}^{M}$. The procedure is as follows. We compute the classical
robust solution $\bm{x}^{r}$ and a corresponding worst-case realization, as
per (12a) shown below, and then set $\rho^{r}$ to be an optimal Lagrange
multiplier of the last inequality in (12b).

$\displaystyle\bm{x}^{r}\in\mathop{\arg\min}_{\bm{x}\in\mathcal{X}}\;\left\\{c(\bm{x})+\max_{\bm{\xi}\in\Xi^{0}}\mathcal{Q}(\bm{x},\bm{\xi})\right\\},\;\bm{\xi}^{r}\in\mathop{\arg\max}_{\bm{\xi}\in\Xi^{0}}\;\mathcal{Q}(\bm{x}^{r},\bm{\xi}),$
(12a)
$\displaystyle\inf_{\bm{y}\in\mathcal{Y},\bm{z}\in[0,1]^{M}}\left\\{\bm{q}(\bm{\xi}^{r})^{\top}\bm{y}:\bm{W}_{0}\bm{y}\geq\bm{T}_{0}\bm{z}+\bm{h}(\bm{x}^{r}),\;(\bm{\mathrm{e}}-2\bm{\xi}^{r})^{\top}\bm{z}+\bm{\mathrm{e}}^{\top}\bm{\xi}^{r}\leq
0\right\\}.$ (12b)

We next prove that our procedure is valid: the optimal Lagrange multiplier
$\rho^{r}$ is indeed an _exact value_ for the penalty parameter.

######  Theorem 4 (Exact penalty reformulation).

The optimal value of the distributionally robust two-stage problem (2) remains
unchanged if we replace the second-stage loss function
$\mathcal{Q}(\bm{x},\bm{\xi})$ with its penalty reformulation
$\mathcal{Q}^{\rho^{r}}(\bm{x},\bm{\xi})$ as defined in (11) with penalty
parameter $\rho^{r}$.

###### Proof.

See Appendix G. ∎

The classical robust formulation (12a) is tractable if we choose
$\Xi^{0}=\\{0,1\\}^{M}$ and the loss function $\mathcal{Q}(\bm{x},\bm{\xi})$
exhibits a down-monotone (or up-monotone) property with respect to the random
parameters $\bm{\xi}$; that is,
$\mathcal{Q}(\bm{x},\bm{\xi}^{\prime})\geq\mathcal{Q}(\bm{x},\bm{\xi})$
whenever $\bm{\xi}^{\prime}\geq\bm{\xi}$. In particular, this is the case for
network optimization problems where removing network elements is never
advantageous. In such cases, the classical robust formulation reduces to a
deterministic problem under the worst-case realization of the uncertain
parameters. Moreover, this worst-case realization is often independent of the
optimal robust solution $\bm{x}^{r}$. For example, suppose
$\bm{\xi}\in\\{0,1\\}^{M}$ is a random vector denoting which of $M$ arcs in a
network have been “disrupted” and arc-flow variables $\bm{y}$ satisfy $0\leq
y_{a}\leq\xi_{a}\bar{y}_{a}$ (i.e., the flow on arc $a$ is zero whenever it is
disrupted $\xi_{a}=0$, and is bounded between $0$ and $\bar{y}_{a}$
otherwise). The worst-case realization is to disrupt all arcs in the network,
$\bm{\xi}^{r}=\bm{0}$, independent of the optimal robust solution
$\bm{x}^{r}$. Notably, this also implies that we can circumvent the
computation of (12a) when determining the value of $\rho^{r}$. We generalize
this observation in the following proposition.

######  Proposition 1.

Assume without loss of generality that
$\mathcal{Y}\subseteq\mathbb{R}_{+}^{N_{2}}$. Suppose also that for all
$j\in[M]$, we have either
$\bm{T}_{0}\bm{\mathrm{e}}_{j},\bm{Q}\bm{\mathrm{e}}_{j}\in\mathbb{R}^{L}_{+}$
or
$-\bm{T}_{0}\bm{\mathrm{e}}_{j},-\bm{Q}\bm{\mathrm{e}}_{j}\in\mathbb{R}^{L}_{+}$.
Then, the classical robust optimization problem
$\min_{\bm{x}\in\mathcal{X}}\;\Big{\\{}c(\bm{x})+\max_{\bm{\xi}\in\\{0,1\\}^{M}}\mathcal{Q}(\bm{x},\bm{\xi})\Big{\\}}$
reduces to a deterministic problem
$\min_{\bm{x}\in\mathcal{X}}\;\left\\{c(\bm{x})+\mathcal{Q}(\bm{x},\bm{\xi}^{r})\right\\}$,
where for each $j\in[M]$ we have $\xi^{r}_{j}=1$ if
$\bm{T}_{0}\bm{\mathrm{e}}_{j},\bm{Q}\bm{\mathrm{e}}_{j}\in\mathbb{R}^{L}_{+}$
and $\xi^{r}_{j}=0$ otherwise.

###### Proof.

See Appendix G. ∎

####  3.2.1 Penalty reformulation for indicator constraints

The penalty reformulation $\mathcal{Q}^{\rho}(\bm{x},\bm{\xi})$ in (11)
introduces $M$ additional variables $\bm{z}\in[0,1]^{M}$ in the second-stage
loss function. This can be avoided by exploiting a structure that is common in
network optimization problems. In several such applications, an uncertain
parameter $\xi_{j}\in\\{0,1\\}$ may switch on and off a single constraint
$f_{j}(\bm{y})\geq 0$, leading to the following definition of the loss
function:

$\mathcal{Q}_{\mathrm{ind}}(\bm{x},\bm{\xi})=\inf_{\bm{y}\in\mathcal{Y}}\left\\{\bm{q}(\bm{\xi})^{\top}\bm{y}:\begin{array}[]{l}\displaystyle\bm{W}_{0}\bm{y}\geq\bm{h}(\bm{x})\\\
\displaystyle\xi_{j}=1\implies\left[f_{j}(\bm{y})=0\right],\;j\in[M]\\\
\displaystyle\xi_{j}=0\implies\left[f_{j}(\bm{y})\geq
0\right],\;j\in[M]\end{array}\right\\},$ (13)

where $f_{j}:\mathcal{Y}\mapsto\mathbb{R}$ is an affine function for each
$j\in[M]$. For example, in a network, $\xi_{j}=1$ may indicate that link $j$
has failed and the corresponding flow variable $f_{j}(\bm{y})=y_{j}$ must be
set to $0$, whereas $\xi_{j}=0$ may indicate that the flow variable $y_{j}$
can take nonzero values. Such constraints are generally written as
$f_{j}(\bm{y})\leq\bar{f}(1-\xi_{j})$, where $\bar{f}$ is a big-M upper bound
on $f_{j}(\bm{y})$. Obtaining tight estimates on $\bar{f}$ may be non-trivial
and at the same time, an overestimation of $\bar{f}$ can lead to numerical
issues. In such cases, we can avoid both _(i)_ estimating $\bar{f}$, and
_(ii)_ introducing auxiliary variables $\bm{z}$ in the penalty reformulation
(11), by simply retaining the constraint $f_{j}(\bm{y})\geq 0$ and adding
$+\rho\xi_{j}f(\bm{y}_{j})$ to the objective, leading to a tighter penalty
reformulation.

######  Corollary 2 (Penalty reformulation of the loss function with indicator
constraints).

There exists a sufficiently large, yet _finite_ , penalty parameter $\rho>0$
such that the second-stage loss function
$\mathcal{Q}_{\mathrm{ind}}(\bm{x},\bm{\xi})$ in (13) is equivalent to

$\mathcal{Q}^{\rho}_{\mathrm{ind}}(\bm{x},\bm{\xi})=\inf_{\bm{y}\in\mathcal{Y}}\left\\{\bm{q}(\bm{\xi})^{\top}\bm{y}+\rho\bm{\xi}^{\top}\bm{f}(\bm{y}):\bm{W}_{0}\bm{y}\geq\bm{h}(\bm{x}),\;\bm{f}(\bm{y})\geq\bm{0}\right\\}.$
(14)

The penalty term $+\rho\xi_{j}f(\bm{y}_{j})$ ensures that $f_{j}(\bm{y})$ is
driven to $0$ whenever $\xi_{j}=1$ for large values of $\rho$. Note that we no
longer need to estimate the big-M upper bounds $\bar{f}$ since the constraint
$f_{j}(\bm{y})\leq\bar{f}(1-\xi_{j})$ is not required anymore.

Our previous results continue to be valid as long as we replace each
occurrence of the penalty term
$(\bm{\mathrm{e}}-2\bm{\xi})^{\top}\bm{z}+\bm{\mathrm{e}}^{\top}\bm{\xi}$ that
multiples $\rho$ in the objective function of (11) with this modification. For
example, suppose that $\bm{f}(\bm{y})=\bm{f}_{0}+\bm{F}\bm{y}$. Then,
Corollaries 1 and 2 imply that the distributionally robust two-stage problem
(2) admits a convex hull reformulation of the form (8), where the function
$Z_{i}:\mathcal{X}\times\mathbb{R}_{+}\mapsto\mathbb{R}$ and the MICP-
representable set $\mathcal{Z}_{i}$ for each $i\in[N]$ are given as follows:

$\displaystyle\displaystyle
Z_{i}(\bm{x},\alpha)=\mathop{\text{maximize}}_{(\bm{\xi},\bm{\lambda},\bm{\mu},\tau)\in\mathop{\mathrm{cl\,conv}}(\mathcal{Z}_{i})}\displaystyle\left\\{\bm{h}(\bm{x})^{\top}\bm{\lambda}+\rho\bm{f}_{0}^{\top}\bm{\xi}-\bm{f}_{0}^{\top}\bm{\mu}-\alpha\tau\right\\}$
$\displaystyle\displaystyle\mathcal{Z}_{i}=\left\\{(\bm{\xi},\bm{\lambda},\bm{\mu},\tau)\in\Xi\times\mathbb{R}^{L}_{+}\times\mathbb{R}^{M}_{+}\times\mathbb{R}_{+}:\begin{array}[]{l}\displaystyle(\bm{\xi}-\hat{\bm{\xi}}^{(i)},\tau)\in\mathcal{C}^{M+1}\\\
\displaystyle\bm{q}_{0}+\bm{Q}\bm{\xi}+\rho\bm{F}^{\top}\bm{\xi}-\bm{W}_{0}^{\top}\bm{\lambda}-\bm{F}^{\top}\bm{\mu}\in\mathcal{Y}^{*}\end{array}\right\\}.$

Similarly, a value for the penalty parameter $\rho$ can be computed as
follows. First, we compute the classical robust solution $\bm{x}^{r}$ and a
corresponding worst-case realization by solving the problems:

$\bm{x}^{r}\in\mathop{\arg\min}_{\bm{x}\in\mathcal{X}}\;\left\\{c(\bm{x})+\max_{\bm{\xi}\in\Xi^{0}}\mathcal{Q}_{\mathrm{ind}}(\bm{x},\bm{\xi})\right\\},\;\bm{\xi}^{r}\in\mathop{\arg\max}_{\bm{\xi}\in\Xi^{0}}\;\mathcal{Q}_{\mathrm{ind}}(\bm{x}^{r},\bm{\xi}).$

Next, we set $\rho^{r}$ to be an optimal Lagrange multiplier of the last
inequality of the following problem:

$\inf_{\bm{y}\in\mathcal{Y}}\left\\{\bm{q}(\bm{\xi}^{r})^{\top}\bm{y}:\bm{W}_{0}\bm{y}\geq\bm{h}(\bm{x}^{r}),\;\bm{f}(\bm{y})\geq\bm{0},\;(\bm{\xi}^{r})^{\top}\bm{f}(\bm{y})\leq
0\right\\}$

The proof for the validity of
$\mathcal{Q}_{\mathrm{ind}}=\mathcal{Q}^{\rho^{r}}_{\mathrm{ind}}$ is similar
to that of Theorem 4 and we omit it for the sake of brevity. Finally, note
that one can avoid the computation of $\bm{x}^{r}$ and $\bm{\xi}^{r}$ (as
before) by examining the structure of $\bm{f}(\bm{y})$ and $\bm{q}(\bm{\xi})$,
see Proposition 1.

###  3.3 Summary and comparison

Table 1 summarizes the main differences between the linearized and penalty-
based MICP reformulations of the sets $\mathcal{Z}_{i}$, $i\in[N]$ appearing
in the convex hull reformulation (9a)–(9d). Notably, the penalty reformulation
adds far fewer variables and constraints. However, it also requires additional
assumptions and computations. In particular, it requires computing a value for
the penalty parameter $\rho$, which may further entail the solution of a
classical robust optimization problem; see (12a). We do not expect this to be
a limitation, however, because the latter will likely reduce to a
deterministic optimization problem for most practical applications.

Table 1: Summary of the MICP representations of $\mathcal{Z}_{i}$ based on the linearized and penalty reformulations. Reformulation | Size | $\bm{W}(\bm{\xi})$ | $\bm{T}(\bm{x})$ |  Requirements
---|---|---|---|---
Linearized (9d-$\ell$) | $O(ML)$ | affine | convex |  _a priori_ bounds on $\bm{\lambda}$ in $\mathcal{Q}_{d}(\bm{x},\bm{\xi})$
Penalty (9d-$\rho$) | $O(M+L)^{*}$ |  constant† |  constant† |  computation of worst-case realization over any superset of $\Xi^{\sharp}$
∗Can be further reduced; see Corollary 2
†Can be relaxed; see discussion following assumption (A3).
♯Can be done in closed form; see Lemma 1 and preceding discussion.

##  4 Lift-and-Project Approximations

The key challenge in solving the convex hull reformulation (8) is the inner
optimization (9a) over the convex hull of the MICP-representable set
$\mathcal{Z}_{i}$, $i\in[N]$. Appendix B shows that although one can tractably
compute these convex hulls in several cases, the problem remains NP-hard even
in benign settings. Therefore, Appendix C presents a Benders scheme, similar
to the ones proposed in [62, 71], to tackle the convex hull constraints. This
scheme iteratively refines an _inner approximation_ of the MICP representation
of $\mathcal{Z}_{i}$. An alternative to solving the convex hull reformulation
(8) is direct solution of the original reformulation (7) using a column-and-
constraint generation scheme [70]. In contrast to the Benders scheme, the
latter models the second-stage loss function $\mathcal{Q}(\bm{x},\bm{\xi})$
via explicit second-stage variables and constraints by implicitly enumerating
$\bm{\xi}\in\Xi$. In both schemes, however, each iteration requires the
solution of $N$ global MICP problems and therefore, they can become
computationally prohibitive. Moreover, intermediate solutions obtained from
early termination provide no guarantees whatsoever since they bound the
optimal value of the distributionally robust problem (2) from below, which
itself is an upper bound on the true (unknown) optimal value.

These observations motivate the development of _tractable outer
approximations_ of the convex hulls of $\mathcal{Z}_{i}$, $i\in[N]$, that
provide not only (i) upper bounds on the optimal value of (2) but also (ii)
guarantees of polynomial time solvability. Our approximations are based on the
following key observations. First, if we suppose that
$\Xi=\\{\bm{\xi}\in\mathbb{Z}^{M}_{+}:\bm{E}\bm{\xi}\leq\bm{f}\\}$ has a given
outer description, then Section 3 establishes that each of the sets
$\mathcal{Z}_{i}$, $i\in[N]$, can be represented as the feasible region of an
MICP as follows:

$\mathcal{Z}=\left\\{\bm{z}\in\mathbb{R}^{n}:\bm{A}\bm{z}-\bm{b}\in\mathcal{K},\;\;z_{j}\in\\{0,1\\}\;j\in[M]\right\\},$
(15)

where, for ease of exposition, we have dropped the subscript $i$ and included
the bounds, $\bm{z}\geq\bm{0}$ and $1\geq z_{j}\coloneqq\xi_{j}$, $j\in[M]$ in
$\bm{A}\bm{z}-\bm{b}\in\mathcal{K}$. For example, in case of (9d-$\ell$), we
have $\bm{z}=(\bm{\xi},\bm{\lambda},\mathop{\text{vec}}(\bm{\Lambda}),\tau)$,
$n=M+L+ML+1$, $\mathcal{K}=\tilde{\mathcal{K}}\times\mathbb{R}^{n}_{+}$, where
$\tilde{\mathcal{K}}=\mathbb{R}_{+}^{F}\times_{i=1}^{3}\mathbb{R}^{LM}_{+}\times\mathcal{C}^{M+1}\times\mathcal{Y}^{*}$
and $F$ is the dimension of $\bm{f}\in\mathbb{R}^{F}$, and
$\bm{A}=[\tilde{\bm{A}}^{\top}\,\mathbf{I}]^{\top}$ and
$\bm{b}=[\tilde{\bm{b}}^{\top}\,\bm{0}^{\top}]^{\top}$ model the right-hand
half of (9d-$\ell$). Similarly, in case of (9d-$\rho$), we have
$\bm{z}=(\bm{\xi},\bm{\lambda},\bm{\mu},\tau)$, $n=2M+L+1$,
$\mathcal{K}=\tilde{\mathcal{K}}\times\mathbb{R}^{n}_{+}$, where
$\tilde{\mathcal{K}}=\mathbb{R}_{+}^{F}\times\mathbb{R}^{M}_{+}\times\mathcal{C}^{M+1}\times\mathcal{Y}^{*}$,
and $\bm{A}=[\tilde{\bm{A}}^{\top}\,\mathbf{I}]^{\top}$,
$\bm{b}=[\tilde{\bm{b}}^{\top}\,\bm{0}^{\top}]^{\top}$ for suitable matrices
$\tilde{\bm{A}}$ and $\tilde{\bm{b}}$.

Second, given an MICP representation such as the above, its convex hull can be
approximated by a hierarchy of increasingly tight _convex relaxations_ ,

$\mathcal{Z}^{0}\supseteq\mathcal{Z}^{1}\supseteq\ldots\supseteq\mathcal{Z}^{M}=\mathop{\mathrm{cl\,conv}}(\mathcal{Z}),$

that converge to the convex hull in $M$ iterations. Here,
$\mathcal{Z}^{0}\coloneqq\left\\{\bm{z}\in\mathbb{R}^{n}:\bm{A}\bm{z}-\bm{b}\in\mathcal{K}\right\\}$
is the continuous relaxation of $\mathcal{Z}$. Several such _sequential
convexification_ hierarchies are known, the most popular ones being those of
[43, 4, 60, 41]. They are based on the concept of _lift-and-project_ and
represent $\mathop{\mathrm{cl\,conv}}(\mathcal{Z})$ as the _projection_ of
another convex set lying in a higher-dimensional space. These hierarchies were
originally proposed for (pure or mixed-) integer linear sets and later
extended to mixed-integer convex sets in [61, 21]. Our proposal is to use an
intermediate relaxation $\mathcal{Z}^{t}$ of any such hierarchy to outer
approximate $\mathop{\mathrm{cl\,conv}}(\mathcal{Z})$, which results in an
outer approximation of the convex hull reformulation (8) and, hence, a
conservative approximation of the distributionally robust two-stage problem
(2). Notably, since we can optimize an objective function over the level-$t$
relaxation in time $n^{O(t)}$ (which is polynomial for fixed $t$), we can also
obtain tight approximations of the original problem (2) in polynomial time.
The approximation can be refined, if desired, by using higher values of $t$.

Third, the approximation of $\mathop{\mathrm{cl\,conv}}(\mathcal{Z})$ when
used in the convex hull reformulation allows us to dualize the inner
optimization in (9a) using conic duality. The result is a single-stage convex
conic optimization model that can be solved with off-the-shelf solvers.
Notably, we can prove that the resulting approximations of the
distributionally robust two-stage problem (2), obtained by replacing
$\mathop{\mathrm{cl\,conv}}(\mathcal{Z})$ with any of the relaxations
$\mathcal{Z}^{0},\ldots,\mathcal{Z}^{M}$, become exact if the radius
$\varepsilon$ of the Wasserstein ambiguity set $\mathcal{P}$ shrinks to zero
with increasing sample size $N$.

######  Theorem 5 (Lift-and-project approximation quality).

Suppose that $N\varepsilon_{N}\to 0$ as $N\to\infty$. Then, the optimal value
of the distributionally robust two-stage problem (2) coincides with that of
the convex hull reformulation (8) even if we approximate each
$\mathop{\mathrm{cl\,conv}}(\mathcal{Z}_{i})$, $i\in[N]$, with
$\mathcal{Z}_{i}^{0}$ in (9a).

###### Proof.

See Appendix G. ∎

We emphasize that our method is not tied to any particular convexification
technique. This feature is important because each technique has its advantages
and disadvantages. For example, in the linear case (i.e.,
$\mathcal{K}=\mathbb{R}^{n^{\prime}}_{+}$), it is known [42] that the
approximations in order of decreasing tightness are those of [41], [60], [43],
and [4]; however, this ranking is reversed when they are ordered based on
increasing computational complexity. For its simplicity and tradeoff between
tightness and tractability, we focus on the Lovász-Schrijver approximation
[43] in the remainder of this section. We show how it can be used to obtain a
single-stage approximation of the distributionally robust two-stage problem
(2), and we provide practical guidelines for its efficient computation.

###  4.1 Lovász-Schrijver approximation

The level-$1$ Lovász-Schrijver approximation $\mathcal{Z}^{1}$ is defined as a
set-valued mapping, and the level-$t$ approximation $\mathcal{Z}^{t}$ is
defined as an iterated application of this mapping. For any
$u\in\mathbb{R}_{+}$ and any conic representable set such as the continuous
relaxation
$\mathcal{Z}^{0}=\left\\{\bm{z}\in\mathbb{R}^{n}:\bm{A}\bm{z}-\bm{b}\in\mathcal{K}\right\\}$,
we denote by
$\mathcal{Z}^{0}(u)=\left\\{\bm{z}\in\mathbb{R}^{n}:\bm{A}\bm{z}-\bm{b}u\in\mathcal{K}\right\\}$
to be the homogenization of $\mathcal{Z}^{0}$ with respect to $u$. Next, we
define the following lifted set:

$\mathcal{L}(\mathcal{Z}^{0})=\left\\{\bm{z},\\{\bm{z}^{j0}\\}_{j\in[M]},\\{\bm{z}^{j1}\\}_{j\in[M]}\in\mathbb{R}^{n}:\begin{array}[]{l@{\;}l}\displaystyle\exists
u^{j0},u^{j1}\geq 0,\;u^{j0}+u^{j1}=1,&j\in[M]\\\
\displaystyle\bm{z}^{j0}\in\mathcal{Z}^{0}(u^{j0}),\;\bm{z}^{j1}\in\mathcal{Z}^{0}(u^{j1}),&j\in[M]\\\
\displaystyle\bm{z}=\bm{z}^{j0}+\bm{z}^{j1},&j\in[M]\\\ \displaystyle
z^{j0}_{j}=0,\;z^{j1}_{j}=u^{j1},&j\in[M]\\\ \displaystyle
z^{j1}_{k}=z^{k1}_{j},\;\;k\in[M]:k>j,&j\in[M]\end{array}\right\\}.$ (16)

Consider now the following set-valued map, which is the projection of
$\mathcal{L}(\mathcal{Z}^{0})$ onto $\mathbb{R}^{n}$:

$\mathcal{P}(\mathcal{Z}^{0})=\left\\{\bm{z}\in\mathbb{R}^{n}:\exists\bm{z}^{j0},\bm{z}^{j1},\;j\in[M]\text{
such that
}\left(\bm{z},\\{\bm{z}^{j0}\\}_{j\in[M]},\\{\bm{z}^{j1}\\}_{j\in[M]}\right)\in\mathcal{L}(\mathcal{Z}^{0})\right\\}.$
(17)

One can easily verify that $\mathcal{P}(\mathcal{Z}^{0})$ is a convex
relaxation of $\mathop{\mathrm{cl\,conv}}(\mathcal{Z})$. In fact, we have the
following relationship [21, Theorem 1]:
$\mathop{\mathrm{cl\,conv}}(\mathcal{Z})\subseteq\mathcal{P}(\mathcal{Z}^{0})\subseteq\bigcap_{j\in[M]}\mathop{\mathrm{cl\,conv}}(\left\\{\bm{z}\in\mathcal{Z}^{0}:z_{j}\in\\{0,1\\}\right\\})\subseteq\mathcal{Z}^{0}.$
The set $\mathcal{P}(\mathcal{Z}^{0})$ corresponds to the level-1 relaxations
of the Lovász-Schrijver hierarchy. For any $t\geq 1$, the level-$t$ relaxation
is given by $\mathcal{Z}^{t}=\mathcal{P}(\mathcal{Z}^{t-1})$, and one can show
that $\mathcal{Z}^{M}=\mathop{\mathrm{cl\,conv}}(\mathcal{Z})$. This is known
as the _linear_ Lovász-Schrijver hierarchy. One can obtain stronger
relaxations by imposing positive semidefiniteness on the submatrix of
$\\{\bm{z}^{j1}\\}_{j\in[M]}$ corresponding to the binary variables.

######  Remark 2 (Positive semidefinite Lovász-Schrijver hierarchy).

Consider the following set:

$\mathcal{P}_{+}(\mathcal{Z}^{0})=\left\\{\bm{z}\in\mathbb{R}^{n}:\begin{array}[]{l}\exists\bm{z}^{j0},\bm{z}^{j1},\;j\in[M]\text{
such that
}\left(\bm{z},\\{\bm{z}^{j0}\\}_{j\in[M]},\\{\bm{z}^{j1}\\}_{j\in[M]}\right)\in\mathcal{L}(\mathcal{Z}^{0})\\\
\text{and
}\left[\bm{\xi}^{11}\ldots\bm{\xi}^{M1}\right]\succeq\bm{\xi}\bm{\xi}^{\top},\text{where
}\bm{\xi}=[z_{1}\ldots z_{M}]^{\top},\bm{\xi}^{j1}=[z^{j1}_{1}\ldots
z^{j1}_{M}]^{\top}\end{array}\right\\}.$

This set corresponds to the level-1 relaxation of the _positive semidefinite_
Lovász-Schrijver hierarchy and is related to its linear counterpart as
follows:
$\mathop{\mathrm{cl\,conv}}(\mathcal{Z})\subseteq\mathcal{P}_{+}(\mathcal{Z}^{0})\subseteq\mathcal{P}(\mathcal{Z}^{0})$.
For any $t\geq 1$, the level-$t$ relaxation in this hierarchy is given by
$\mathcal{Z}^{t}=\mathcal{P}_{+}(\mathcal{Z}^{t-1})$. As before, we have
$\mathcal{Z}^{M}=\mathop{\mathrm{cl\,conv}}(\mathcal{Z})$.

For a given MICP representation of $\mathcal{Z}_{i}$ and any $t\geq 1$, we can
use the level-$t$ Lovász-Schrijver relaxation to approximate
$\mathop{\mathrm{cl\,conv}}(\mathcal{Z}_{i})$ in (9a)–(9d). We can then
convert the inner maximization (9a) to a minimization and embed it in the
first-stage problem. The following lemma illustrates this for $t=1$; we omit
the proof since it follows from a straightforward application of conic
duality. In this lemma, $\bm{\gamma}$ represents the objective function of the
inner optimization (9a); in the linearized MICP representation (9d-$\ell$) of
$\mathcal{Z}$, it is given by
$\bm{\gamma}=\left[\bm{0}^{\top}\;\bm{h}(\bm{x})^{\top}\,\mathop{\text{vec}}(\bm{T}(\bm{x}))^{\top}\,-\alpha\right]^{\top}$,
whereas in the penalty MICP representation (9d-$\rho$), we have
$\bm{\gamma}=\left[\rho\bm{\mathrm{e}}^{\top}\;\bm{h}(\bm{x})^{\top}\,-\bm{\mathrm{e}}^{\top}\,-\alpha\right]^{\top}$.

######  Lemma 2.

Suppose that $\mathcal{Z}$ is defined as in (15) with
$\mathcal{K}=\tilde{\mathcal{K}}\times\mathbb{R}^{n}_{+}$,
$\bm{A}=[\tilde{\bm{A}}^{\top}\,\mathbf{I}]^{\top}$ and
$\bm{b}=[\tilde{\bm{b}}^{\top}\,\bm{0}^{\top}]^{\top}$ for suitable matrices
$\tilde{\bm{A}}$ and $\tilde{\bm{b}}$. If $\mathcal{Z}^{1}$ denotes the
level-1 Lovász-Schrijver relaxation of
$\mathop{\mathrm{cl\,conv}}(\mathcal{Z})$, then for any
$\bm{\gamma}\in\mathbb{R}^{n}$, we have the following strong duality result:

$\begin{array}[]{r@{\;\;}l}\displaystyle\sup_{\bm{z}\in\mathcal{Z}^{1}}\,\bm{\gamma}^{\top}\bm{z}\;=\;\mathop{\text{minimize}}&\displaystyle\sum_{j\in[M]}\max\left\\{-\tilde{\bm{b}}^{\top}\tilde{\bm{\zeta}}^{j0},-\tilde{\bm{b}}^{\top}\tilde{\bm{\zeta}}^{j1}-\beta^{j1}\right\\}\\\
\text{subject
to}&\displaystyle\sum_{j\in[M]}\bm{\delta}^{j}+\bm{\gamma}\leq\bm{0},\\\
&\displaystyle\begin{drcases}\tilde{\bm{A}}^{\top}\tilde{\bm{\zeta}}^{j\ell}+\beta^{j\ell}\bm{\mathrm{e}}_{j}-\mathbb{I}[\ell=1]\bm{Y}\bm{\mathrm{e}}_{j}\leq\bm{\delta}^{j}\\\
\bm{\delta}^{j}\in\mathbb{R}^{n},\;\tilde{\bm{\zeta}}^{j\ell}\in\tilde{\mathcal{K}}^{*},\;\beta^{j\ell}\in\mathbb{R}\end{drcases}\;\forall\ell\in\\{0,1\\},j\in[M],\\\
&\bm{Y}\in\mathbb{R}^{M\times M}:\bm{Y}=-\bm{Y}^{\top}.\end{array}$ (18)

######  Remark 3 (Relationship to approximations obtained by relaxing the
support).

An alternative outer approximation of the distributionally robust two-stage
problem (2) can be obtained by simply relaxing the zero-one constraints on the
support in reformulation (7); i.e., by replacing
$\Xi=\\{\bm{\xi}\in\mathbb{Z}^{M}_{+}:\bm{E}\bm{\xi}\leq\bm{f}\\}$ in (7) with
its continuous relaxation
$\\{\bm{\xi}\in\mathbb{R}^{M}_{+}:\bm{E}\bm{\xi}\leq\bm{f}\\}$. The resulting
approximation is intractable in general, unless the uncertainty $\bm{\xi}$
appears only in the objective of the loss function
$\mathcal{Q}(\bm{x},\bm{\xi})$ (see discussion in Section 3). In the latter
case, it can be easily seen that the resulting approximation coincides
precisely with that obtained by replacing
$\mathop{\mathrm{cl\,conv}}(\mathcal{Z}_{i})$ in (10a)–(10b) with the
continuous relaxation $\mathcal{Z}_{i}^{0}$. Therefore, our lift-and-project
approximation technique can be viewed as a generalization of this approach.
Unlike the former, however, a crucial difference is that this approach
provides no formal mechanism to improve the quality of the final
approximation; we illustrate this empirically in Section 5.

###  4.2 Numerical considerations

Several factors can impact the numerical solution of the approximations
obtained by replacing $\mathop{\mathrm{cl\,conv}}(\mathcal{Z}_{i})$,
$i\in[N]$, in (9a), with their level-1 Lovász-Schrijver relaxations
$\mathcal{Z}^{1}_{i}$. First, if we use the linearized MICP representation
(9d-$\ell$) to reformulate $\mathcal{Z}_{i}$, then formulation (18) has
$O(M(F+N_{2}+LM))$ variables, $O(M^{2}L)$ linear constraints, and
$O(M(M+N_{2}))$ conic constraints. If we use the penalty representation
(9d-$\rho$), then these are reduced to $O(M(F+N_{2}+L+M))$ variables,
$O(M(L+M))$ linear constraints, and $O(M(M+N_{2}))$ conic constraints. In
either case, the number of conic constraints can be reduced to $O(MN_{2})$
whenever the metric
$d(\bm{\xi},\bm{\xi}^{\prime})=\lVert\bm{\xi}-\bm{\xi}^{\prime}\rVert_{1}$ is
induced by the 1-norm (see argument in proof of Proposition 3 in Appendix B).

Second, observe that setting $\bm{Y}=\bm{0}$ also reduces the number of
variables in (18) for only a minor loss in approximation quality; indeed, the
resulting relaxation is still equal to
$\bigcap_{j\in[M]}\mathop{\mathrm{cl\,conv}}(\\{\bm{z}\in\mathcal{Z}^{0}_{i}:z_{j}\in\\{0,1\\}\\})$.
In fact, it is equal to
$\bigcap_{j\in\mathcal{J}_{i}}\mathop{\mathrm{cl\,conv}}(\\{\bm{z}\in\mathcal{Z}^{0}_{i}:z_{j}\in\\{0,1\\}\\})$,
where $\mathcal{J}_{i}\subseteq[M]$ is the index set of binary parameters
whose optimal values in the left-hand side of (18) are fractional. We expect
$\lvert\mathcal{J}_{i}\rvert$ to be small since the optimal value of
$\bm{\xi}$ in any convex relaxation of $\mathcal{Z}_{i}$ is unlikely to be far
from the binary-valued $\hat{\bm{\xi}}^{(i)}$ (see argument in proof of
Theorem 5). This motivates the following iterative heuristic to identify the
index sets $\mathcal{J}_{i}$. Note that this procedure is independent of the
MICP representation used for each $\mathcal{Z}_{i}$.

  1. 1.

Select $\texttt{tol}\in(0,0.5)$, $\texttt{niter}\in\mathbb{Z}_{+}$. Set
$\texttt{iter}\leftarrow 1$. For each $i\in[N]$, set
$\tilde{\mathcal{Z}}_{i}^{1}\leftarrow\mathcal{Z}^{0}_{i}$ and
$\mathcal{J}_{i}\leftarrow\emptyset$.

  2. 2.

For each $i\in[N]$, replace $\mathop{\mathrm{cl\,conv}}(\mathcal{Z}_{i})$ with
its current approximation $\tilde{\mathcal{Z}}_{i}^{1}$ and dualize the
corresponding problem (9a). Solve the resulting convex hull approximation (8).

  3. 3.

For each $i\in[N]$, let $\bar{\bm{\xi}}^{[i]}$ be the optimal value of
$\bm{\xi}$ in (9a), recovered as scaled dual multipliers. For each
$j\in[M]\setminus\mathcal{J}_{i}$, if
$\bar{{\xi}}^{[i]}_{j}\in[\texttt{tol},1-\texttt{tol}]$, update
$\mathcal{J}_{i}\leftarrow\mathcal{J}_{i}\cup\\{j\\}$ and
$\tilde{\mathcal{Z}}_{i}^{1}$ as follows:

$\tilde{\mathcal{Z}}_{i}^{1}\leftarrow\left\\{\bm{z}\in\mathbb{R}^{n}:\begin{array}[]{l@{\;}l}\displaystyle\exists
u^{j0},u^{j1}\geq 0,\;u^{j0}+u^{j1}=1,&j\in\mathcal{J}_{i}\\\
\displaystyle\exists\bm{z}^{j0}\in\mathcal{Z}^{0}_{i}(u^{j0}),\;\bm{z}^{j1}\in\mathcal{Z}^{0}_{i}(u^{j1}),&j\in\mathcal{J}_{i}\\\
\displaystyle\bm{z}=\bm{z}^{j0}+\bm{z}^{j1},&j\in\mathcal{J}_{i}\\\
\displaystyle
z^{j0}_{j}=0,\;z^{j1}_{j}=u^{j1},&j\in\mathcal{J}_{i}\end{array}\right\\}.$

  4. 4.

If none of the index sets $\mathcal{J}_{1},\ldots,\mathcal{J}_{N}$ were
updated or if $\texttt{iter}\geq\texttt{niter}$, stop. Otherwise, update
$\texttt{iter}\leftarrow\texttt{iter}+1$ and go to Step 2.

Note that the successive optimizations in Step 2 can benefit from an efficient
initialization of their variables by using the optimal solution from the
previous solve. Moreover, the size of these problems can be controlled by
using smaller values of niter and larger values of tol, since they directly
influence the size of $\mathcal{J}_{i}$ and $\tilde{\mathcal{Z}}_{i}^{1}$,
albeit at the expense of coarser approximations. In our implementation, we
found that a setting of $\texttt{iterlim}=5$ and $\texttt{tol}=10^{-2}$
achieved a good tradeoff between approximation quality and computational
effort.

##  5 Computational Experiments

We illustrate the applicability of our method to operational problems in
electric power systems in Section 5.1, and to design problems in multi-
commodity flow networks in Section 5.2. Our goals are to: (i) study the lift-
and-project approximations $\mathcal{Z}^{0}$ and $\mathcal{Z}^{1}$ in terms of
their computational effort and ability to approximate
$\mathop{\mathrm{cl\,conv}}(\mathcal{Z})$; (ii) compare their out-of-sample
performance with the standard sample average approximation and with classical
two-stage robust optimization; and, (iii) elucidate the effect of two key
parameters on the relative benefits of the distributionally robust two-stage
problem (2) over these classical formulations: the “rareness” of network
failures and the relative magnitude of “impact” when failures occur.

Our code was implemented in Julia 1.5.3, using JuMP 0.21.4. We used Mosek 9.2
for solving our lift-and-project approximations, and Gurobi 9.1.1 as the
solver for the Benders and column-and-constraint generation schemes (which we
compare in Sections 5.1 and 5.2 respectively), since the latter performed
better than the former in solving the mixed-integer subproblems in those
schemes; whereas Mosek performed better in solving the conic programming
relaxations. All runs were conducted on an Intel Xeon 2.3 GHz computer, with a
limit of four cores per run.

###  5.1 Optimal power flow

We use our method to address the security-constrained optimal power flow
problem that is fundamental to the secure operation of electric power grids
and solved every fifteen minutes or so by grid operators (e.g.,see [1, 22]).
The goal is to determine voltages and generation levels of available
generators so as to satisfy power demand in the network, while adhering to
various physical and engineering constraints. For example, electric power
between network nodes (also known as _buses_) can flow only along capacitated
edges or transmission lines. As such, the latter are failure prone, and
transmission line outages can lead to an unstable power network or even
complete system failure, resulting in costly blackouts. However, such high-
impact failure events are rare. For example, between the years 2000 and 2014,
fewer than 1,500 power outages have occurred that affected 50,000 or more
residents in the entire United States, which is fewer than 100 events per year
[64]. This rarity complicates the accurate estimation of their underlying
distribution.

Because electric power is governed by complex physical laws, optimal power
flow is a highly nonlinear optimization problem. Nevertheless, the underlying
physics can be approximated well by using second-order cone or semidefinite
programming relaxations [44]. Although our method generalizes to any convex
cone relaxation, we focus on the standard second-order cone relaxation [38],
where $\mathcal{X}$ is second-order cone representable and
$\mathcal{Q}(\bm{x},\bm{\xi})$ is the optimal value of a second-order cone
program.

Our two-stage optimization model is inspired by [63] and presented in Appendix
E. Conceptually, the first-stage problem determines minimum cost power
generation levels assuming no line outages. Upon line failure, the second-
stage model seeks to adjust the power generation levels subject to physical
constraints where failed lines cannot be used, with a goal of minimizing the
total penalty cost of violating power balances. This model satisfies
assumptions (A1), (A2), and (A3) and allows the use of the penalty
reformulation (9d-$\rho$), which also has the advantage of using fewer
variables and constraints compared with the linearized reformulation (see
Section 3.3).

The operational state of transmission lines is modeled as a random binary
vector ${\bm{\xi}}$ with support $\Xi=\\{0,1\\}^{M}$, where ${\xi}_{l}=1$
indicates that line $l$ has failed. In particular, since ${\bm{\xi}}$
represent on/off switches, we can use Corollary 2 to get not only a smaller
MICP formulation but also tighter values of the penalty parameter
$\rho=\rho^{r}$. The latter is computed by using Theorem 4, where the
classical robust counterpart reduces to a deterministic problem (see Lemma 1);
indeed, the second-stage loss function trivially attains its worst-case value
when each component of $\bm{\xi}$ is one, that is, when all transmission lines
fail.

####  5.1.1 Test instances

We conduct our experiments on the standard IEEE 14-, 30- and 118-bus test
cases from the PGLib-OPF library [3]. In each case, the second-stage per-unit
penalty cost for violating the power balance equations is set to be $\phi$
times the maximum per-unit first-stage generation cost. Note that this choice
depends on the economic cost of failure to meet power demand. Since loss of
power and blackouts tend to be costly and the associated penalty costs much
larger than the cost of generation, we set $\phi=100$ in Sections 5.1.2 and
5.1.3 and analyze its sensitivity in Section 5.1.4.

We generate empirical data using a Bernoulli model. Specifically, we model
each component of $\tilde{\bm{\xi}}$ as independent and identically
distributed Bernoulli random variables with parameter $\psi M^{-1}$, where $M$
denotes the number of transmission lines. Note that this choice reflects the
rare nature of line failures; in particular, it implies that only $\psi\cdot
100\%$ of training samples record a line failure. We set $\psi=0.1$ in
Sections 5.1.2 and 5.1.3 and analyze its impact in Section 5.1.4.

In all our experiments, for a fixed sample size $N$ and radius $\varepsilon$,
we report average results using 100 statistically independent sets of training
samples, and we estimate the variance by reporting the standard deviation over
these 100 runs. In Sections 5.1.3 and 5.1.4, the out-of-sample performances of
a candidate solution are estimated by using 1,000 statistically independent
sets of testing samples.

####  5.1.2 Approximation quality and computational effort

To study the quality of our lift-and-project approximations, we compute the
following quantities for each sample size $N\in\\{10,100,1000\\}$ and radius
$\varepsilon=\nu\sqrt{N^{-1}\log(N+1)}$, where
$\nu\in\\{0,10^{-3},10^{-2},10^{-1}\\}$: _(i)_ the optimal value $v^{\star}$
of the convex hull reformulation (8) using the Benders scheme described in
Appendix C, and _(ii)_ the optimal values $v^{0}$, $\tilde{v}^{1}$ and $v^{1}$
of formulation (8) when the convex hulls
$\mathop{\mathrm{cl\,conv}}(\mathcal{Z}_{i})$ in (9a)–(9d) are approximated by
using the continuous relaxation $\mathcal{Z}^{0}$ and heuristically and
exactly computed level-1 Lovász-Schrijver relaxations
$\tilde{\mathcal{Z}}^{1}$ and ${\mathcal{Z}}^{1}$, respectively (see Section
4.2). The choice of the radius $\varepsilon=\nu\sqrt{N^{-1}\log(N+1)}$ is
motivated from Theorem 1, and we elaborate on it further in the next
subsection.

Figure 4 reports the average (line plot) and standard deviation (error bar) of
the optimality gaps, defined as $\left(v-v^{\star}\right)/{v^{\star}}\times
100\%$, where $v\in\\{v^{0},\tilde{v}^{1},v^{1}\\}$, based on 100
statistically independent sets of training samples. We make the following
observations.

  * •

The exact level-1 relaxation ${\mathcal{Z}}^{1}$ is near optimal, with
optimality gaps never exceeding $10\%$, whereas the continuous relaxation
$\mathcal{Z}^{0}$ is less accurate, especially for larger radii (_e.g._ , 50%
gap for $\nu=0.1$). The heuristically computed level-1 relaxation
$\tilde{\mathcal{Z}}^{1}$ is also near optimal for small and large radii but
performs relatively poorly for intermediate values of $\nu$.

  * •

For a fixed sample size $N$ and decreasing radius $\nu$, the optimality gaps
of all approximations decrease to $0$. For increasing $\nu$, the gaps of the
level-1 relaxations increase far less rapidly than that of the continuous
relaxation.

  * •

For a fixed radius $\nu$ and increasing sample size $N$, the gaps of all
approximations decrease.

(a) 14-bus, $N=10$

(b) 14-bus, $N=100$

(c) 14-bus, $N=1000$

(d) 30-bus, $N=10$

(e) 30-bus, $N=100$

(f) 30-bus, $N=1000$

Figure 4: Optimality gaps using the continuous relaxation $\mathcal{Z}^{0}$
and the heuristically and exactly computed level-1 Lovász-Schrijver
relaxations $\tilde{\mathcal{Z}}^{1}$ and ${\mathcal{Z}}^{1}$ as a function of
$\nu$ and $N$, where $\varepsilon=\nu\sqrt{N^{-1}\log(N+1)}$.

Figure 5 reports the average computation time to solve the various
approximations and compares them with that of the Benders decomposition
scheme. We offer the following comments.

  * •

The continuous relaxation $\mathcal{Z}^{0}$ and heuristically computed level-1
relaxations $\tilde{\mathcal{Z}}^{1}$ have the smallest computation times,
with the former being faster for larger values of $N$ and $\nu$ and, in
particular, for the larger 118-bus case where it is more than 10 times faster.
When compared with the Benders scheme, the relative difference in their
computation times is minor for small sample sizes $N$ but increases
significantly for large sample sizes. For $N=1000$, both approximations run 10
times faster than the Benders scheme for the 14-bus case, while
$\tilde{\mathcal{Z}}^{1}$ runs 4 times faster and $\mathcal{Z}^{0}$ runs
almost 100 times faster for the 30-bus case.

  * •

The exact level-1 relaxation ${\mathcal{Z}}^{1}$ and the Benders scheme appear
to be the most difficult to solve. Although not shown, for the 30-bus case the
former took about 1, 2 and 10 minutes for $N=10,100$, and $1,000$,
respectively, whereas for the larger 118-bus case neither scheme terminated
within 10 minutes for $N=10,100$ or within 1 hour for $N=1000$. Moreover, some
of the MICP subproblems within the Benders scheme can cause slow convergence
(e.g., due to search tree enumeration). This is evidenced by the fact that
about 1% of the Benders runs did not terminate within 10 minutes even for the
smaller 14-bus and 30-bus cases, respectively.

  * •

In conjunction with Figure 5, the heuristically computed level-1 relaxation
$\tilde{\mathcal{Z}}^{1}$ appears to offer the best tradeoff in terms of
approximation quality and computational effort.

  * •

The run times of all approximations, and in particular
$\tilde{\mathcal{Z}}^{1}$, can be significantly improved by using an efficient
initialization of their variables (see Section 4.2), which we did not
implement.

(a) 14-bus, $N=10$

(b) 14-bus, $N=100$

(c) 14-bus, $N=1000$

(d) 30-bus, $N=10$

(e) 30-bus, $N=100$

(f) 30-bus, $N=1000$

(g) 118-bus, $N=10$

(h) 118-bus, $N=100$

(i) 118-bus, $N=1000$

Figure 5: Computation times for solving formulation (8) using the continuous
relaxation $\mathcal{Z}^{0}$, the heuristically computed and exact level-1
Lovász-Schrijver relaxations $\tilde{\mathcal{Z}}^{1}$ and
${\mathcal{Z}}^{1}$, and the Benders decomposition scheme, as a function of
$\nu$ and $N$, where $\varepsilon=\nu\sqrt{N^{-1}\log(N+1)}$.

####  5.1.3 Out-of-sample performance and finite sample guarantee

To understand the potential benefits of a distributionally robust approach, we
evaluate its out-of-sample performance. For a given training dataset of size
$N$ and a given choice of $\nu$, we obtain a candidate first-stage solution
$\bm{x}^{\nu}$ by solving formulation (8) with the level-1 Lovász-Schrijver
relaxation $\tilde{\mathcal{Z}}^{1}$ and Wasserstein radius
$\varepsilon=\nu\sqrt{N^{-1}\log(N+1)}$. We then estimate the out-of-sample
performance of $\bm{x}^{\nu}$ by recording
$z^{\nu}=\bm{c}(\bm{x}^{\nu})+1000^{-1}\sum_{i=1}^{1000}\mathcal{Q}(\bm{x}^{\nu},\hat{\bm{\xi}}^{(i)})$,
where $\hat{\bm{\xi}}^{(1)},\ldots,\hat{\bm{\xi}}^{(1000)}$ are $1,000$
independently generated testing samples. This entire process is repeated 100
times for statistically independent sets of $N$ training samples and 1,000
testing samples.

We first justify the choice of the radius
$\varepsilon=\nu\sqrt{N^{-1}\log(N+1)}$ as a function of the training sample
size $N$. This dependence is motivated by inequality (5) in Theorem 1.
However, the latter inequality can be loose, especially when accounting for
problem-dependent constants such as the size and diameter of the support
$\Xi$. Therefore, we empirically verify the finite sample guarantee of the
first-stage solutions $\bm{x}^{\nu}$ under the tighter parameterization
$\varepsilon=\nu\sqrt{N^{-1}\log(N+1)}$ for several choices of the coefficient
$\nu\in\\{0,10^{-3},10^{-2},10^{-1}\\}$. Figure 6 reports the _reliability_ of
$\bm{x}^{\nu}$, which we define as the empirical probability (over the 100
sets of training samples) that the optimal value $\tilde{v}^{1}$ of the
level-1 Lovász-Schrijver relaxation $\tilde{\mathcal{Z}}^{1}$ is an upper
bound on the out-of-sample cost $z^{\nu}$.

Figure 6 shows that, for fixed values of $N$, the reliability of
$\bm{x}^{\nu}$ increases with increasing values of $\nu$, and this can be used
to guide the choice of $\nu$. For example, depending on their risk level,
decision-makers can select the smallest value of $\nu$ with sufficiently high
reliability. In particular, for training datasets with small sample size $N$,
observe that $\nu=10^{-2}$ yields an upper bound on the out-of-sample cost
with probability more than 0.5, whereas $\nu=10^{-1}$ yields an upper bound
with probability 1.0. However, note that we cannot always access the true out-
of-sample cost (and hence, the true reliability). Nevertheless, one could
estimate the out-of-sample cost by using cross validation or the holdout
method (e.g, see [50, 37]), and then select a value for the coefficient $\nu$
that offers a good trade-off between low out-of-sample cost and high
reliability.

(a) 14-bus

(b) 30-bus

(c) 118-bus

Figure 6: Reliability of the level-1 Lovász-Schrijver relaxation
$\tilde{\mathcal{Z}}^{1}$, as a function of training sample size $N$.

We now evaluate the benefits of our distributionally robust model over the
sample average approximation, by computing the relative improvement in out-of-
sample cost, which we define as $(z^{0}-z^{\nu})/z^{0}\times 100\%$. Figure 7
reports the mean (solid line) and standard deviation (shaded ribbon) of the
relative improvement over the 100 independent sets of training samples. We
make the following observations from Figure 7.

  * •

The distributionally robust model (2) consistently outperforms the sample
average approximation, particularly for small sample sizes $N$. The magnitude
of the relative improvement is instance dependent (roughly 15%, 10%, and 5%
for the 14-, 30-, and 118-bus cases, respectively) but consistently decreases
for large values of $N$ as expected. The magnitude of the radius that leads to
the best possible improvement also is instance dependent.

  * •

The larger variances in improvement for smaller $N$ and for larger instances
can be partially explained by the combinatorial growth in the number of truly
distinct training datasets of size $N$ (i.e., those that lead to distinct
first-stage solutions) that are possible under the rare event model of line
outages. The large variances for $\nu=10^{-1}$ can also be similarly explained
by the larger number of truly distinct first-stage solutions that can result
from slight variations in the training dataset.

(a) 14-bus

(b) 30-bus

(c) 118-bus

Figure 7: Relative improvement in the out-of-sample performance of the
distributionally robust two-stage model (2) when compared with the sample
average approximation, as a function of sample size $N$.

We now compare the out-of-sample performance of the level-1 Lovász-Schrijver
relaxation $\tilde{\mathcal{Z}}^{1}$ corresponding to the best choice of $\nu$
(which is $10^{-1}$ for the 14- and 30-bus cases, and $10^{-2}$ for the
118-bus case) with (i) the continuous relaxation $\mathcal{Z}^{0}$ for the
same $\nu$, and (ii) a classical two-stage robust optimization model of the
following form:

$\mathop{\text{minimize}}_{\bm{x}\in\mathcal{X}}\;c(\bm{x})+\max_{\bm{\xi}\in\Xi_{K}}\mathcal{Q}(\bm{x},\bm{\xi}),$
(19)

where $\Xi_{K}=\\{\bm{\xi}\in\\{0,1\\}^{M}:\xi_{1}+\ldots+\xi_{M}\leq K\\}$ is
the uncertainty set with $K$ being the _budget_ of uncertainty. In particular,
for each instance, we let $K\in\\{0,1,2,5,10\\}$ and obtain optimal first-
stage solutions $\bm{x}^{K}$ of problem (19) with the Benders decomposition
scheme. As before, we estimate the out-of-sample performance of $\bm{x}^{K}$
as
$z^{K}=\bm{c}(\bm{x}^{K})+1000^{-1}\sum_{i=1}^{1000}\mathcal{Q}(\bm{x}^{K},\hat{\bm{\xi}}^{(i)})$,
where $\hat{\bm{\xi}}^{(1)},\ldots,\hat{\bm{\xi}}^{(1000)}$ are the same
$1,000$ testing samples used to estimate $z^{\nu}$. We then record the best
possible $K$ yielding the lowest $z^{K}$, and compute the relative improvement
in out-of-sample cost, which we define as $(z^{K}-z^{\nu})/z^{\nu}\times
100\%$. Figure 8 reports the mean (solid line) and standard deviation (shaded
ribbon) of the relative improvement over the 100 independent sets of training
samples. We make the following observations from Figure 8.

  * •

The distributionally robust model strongly outperforms its classical robust
counterpart, across all instances, with relative improvements of 5%, 1% and
15% for the 14-, 30- and 118-bus cases, respectively. In contrast to Figure 7,
the relative improvements increase with increasing values of $N$; this is
expected since the classical robust model (19) ignores all sample data and
therefore, it becomes overly conservative in the presence of a moderate amount
of data. Thus, we observe that for small to moderate values of $N$, the
distributionally robust model improves upon both the sample average
approximation and classical robust optimization. Finally, although not shown,
solving the classical robust model with the Benders scheme took longer than
solving the level-1 relaxation, especially for the larger 118-bus case.

  * •

The relative improvement of the level-1 relaxation $\tilde{\mathcal{Z}}^{1}$
over the continuous relaxation $\mathcal{Z}^{0}$ is smaller, but can be as
high as 4% as seen in the 14-bus case. It should be noted that these
quantities are necessarily upper bounded by the improvements over the sample
average approximations reported in Figure 7, and can be significant for
applications including optimal power flow, which are executed several times
each day of the year. Finally, we note that the tradeoff between tighter in-
sample optimality gaps (see Figure 4) and hence stronger finite sample
reliability guarantees (see Figure 6) offered by the level-1 relaxation
$\tilde{\mathcal{Z}}^{1}$, with the faster computation times for solving the
continuous relaxation $\mathcal{Z}^{0}$ (see Figure 5), can guide the design
of an algorithmic scheme.

(a) 14-bus

(b) 30-bus

(c) 118-bus

Figure 8: Relative improvement in the out-of-sample performance of the
distributionally robust two-stage model (2) solved using the level-1 Lovász-
Schrijver relaxation $\tilde{\mathcal{Z}}^{1}$, when compared with the
classical two-stage robust optimization model (19), and the continuous
relaxation $\mathcal{Z}^{0}$.

####  5.1.4 Sensitivity analysis

The instance-dependent behavior of the out-of-sample performance from the
previous subsection suggests that it might also be influenced by other
parameters. Here, we investigate the effect of the “rareness” $\psi$ of
transmission line failures and the relative magnitude $\phi$ of “impact” when
failures occur. Recall from Section 5.1.1 that higher values of $\psi$
increase the probability of line failures, whereas higher values of $\phi$
increase the penalty cost for failing to satisfy power demand due to
transmission line failures. Figure 9 shows the relative improvement of the
distributionally robust two-stage model (2) over the sample average
approximation for various choices of $\psi$ and $\phi$. For brevity, we report
results only for the 30-bus instance; the high-level insights do not change
for other instances. We make the following observations from Figure 9.

  * •

For fixed values of the line failure probability $\psi$, Figures 9(a)–9(c)
show that as the impact due to failure $\phi$ increases, the relative benefits
of a distributionally robust approach strongly increase. In other words,
benefits increase with higher impacts of failures. Interestingly, Figure 9(a)
also shows that if failures are rare but low impact, then ignoring them (as in
the sample average approximation) may not incur high out-of-sample costs, even
for small values of $N$.

  * •

For fixed values of the magnitude of impact $\phi$, Figures 9(d)–9(f) show
that as the probability of failures $\psi$ increases, the relative benefits of
a distributionally robust approach increases. However, observe that this does
not necessarily imply that the relative benefits are small when line failure
probabilities are small. Indeed, we observe that the relative benefits remain
as high as 10% even when individual line failure probabilities are less than
$0.05M^{-1}$.

(a) $(\psi,\phi)=(0.1,50)$

(b) $(\psi,\phi)=(0.1,100)$

(c) $(\psi,\phi)=(0.1,200)$

(d) $(\psi,\phi)=(0.05,200)$

(e) $(\psi,\phi)=(0.1,200)$

(f) $(\psi,\phi)=(0.2,200)$

Figure 9: Relative improvement in the out-of-sample performance of the
distributionally robust two-stage model (2) when compared with the sample
average approximation, for various values of $(\psi,\phi)$.

###  5.2 Multi-commodity network design

We now consider multi-commodity network design problems that have applications
in telecommunications, transportation, logistics and production planning,
among others (e.g., see [48, 23]). In several of these applications, it is
required to send flows to satisfy known demands between multiple origin-
destination pairs or commodities. The goal is to minimize the total cost,
which is the sum of fixed costs of installing arc capacities and variable
costs of routing flows for each commodity. As such, failures of network
elements can lead to a reduction in its available flow capacity and a
subsequent failure to meet demands. This is particularly true in
telecommunication networks where the loss of even a single (typically high-
capacity) fiber-optic cable or router equipment can cause a substantial
fraction of the overall flow (e.g., internet traffic) to be lost, leading to
potentially huge economic impacts [47]. Fortunately, these networks are
typically well-engineered and therefore, such high-impact failures are rare.
At the same time, this general lack of failure data in real networks
complicates the accurate estimation of their underlying distribution.

The two-stage optimization model we consider is presented in Appendix F, and
can be described as follows. The first-stage problem determines the arc
capacities that can be used for routing flows. Upon failure, the second-stage
model determines the routing of each commodity along the degraded network
topology constrained by the first-stage arc capacities, and with the objective
of minimizing the sum of variable routing costs and penalty costs for not
satisfying demands.

For ease of exposition, we model only node failures, where ${\bm{\xi}}$ is
supported on $\Xi=\\{0,1\\}^{M}$ and ${\xi}_{i}=1$ indicates that node $i$ has
failed. Since ${\bm{\xi}}$ represent on/off switches, we can employ Corollary
2, and the penalty parameter $\rho=\rho^{r}$ can be computed using Theorem 4.
Here, the classical robust counterpart reduces to a deterministic problem (see
Lemma 1), since the second-stage loss function trivially attains its worst-
case value when each component of $\bm{\xi}$ is one; that is, when all nodes
fail.

We conduct our experiments on the $(20,230,40,V,L)$ instance from the so-
called Class I set of instances in [23]. As the name indicates, the instance
has 20 nodes, 230 arcs, 40 commodities, and the letters $V$ and $L$ indicate
that fixed-costs are relatively low compared to variable costs, and that the
problem is loosely capacitated. We generate empirical data by modeling each
component of $\tilde{\bm{\xi}}$ as independent and identically distributed
Bernoulli random variables with parameter $\psi M^{-1}$, where $\psi=0.1$. As
before, for a fixed sample size $N$ and radius $\varepsilon$, we report
average results using 100 statistically independent sets of training samples,
and we estimate the variance by reporting the standard deviation over these
100 runs. The out-of-sample performances of candidate solutions are estimated
by using 1,000 statistically independent sets of testing samples.

####  5.2.1 Approximation quality and computational effort

Similar to Section 5.1.2, we compute the optimality gaps of the convex hull
reformulation (8) when the convex hulls
$\mathop{\mathrm{cl\,conv}}(\mathcal{Z}_{i})$ in (9a)–(9d) are approximated by
using the continuous relaxation $\mathcal{Z}^{0}$ and heuristically computed
level-1 Lovász-Schrijver relaxation $\tilde{\mathcal{Z}}^{1}$, for sample
sizes $N\in\\{10,100,1000\\}$ and radii
$\varepsilon=\nu\sqrt{N^{-1}\log(N+1)}$,
$\nu\in\\{0,10^{-3},10^{-2},10^{-1}\\}$. In doing so, the true optimal value
of each instance is computed by solving formulation (7) using the column-and-
constraint generation scheme. The mixed-integer subproblems in this scheme are
solved by dualizing the second-stage loss function as opposed to using its
Karush-Kuhn-Tucker conditions, since it results in fewer additional variables
and constraints, see [70]. Any bilinear expressions involving dual variables
and uncertain parameters are reformulated using indicator constraints since
the lack of a priori known upper bounds on the dual variables prohibits direct
linearization using McCormick inequalities.

Figure 10 reports the average (line plot) and standard deviation (error bar)
of the optimality gaps, whereas Figure 11 reports the corresponding
computation times, based on 100 statistically independent sets of training
samples. Figure 10 shows that both the continuous and level-1 relaxations
provide very similar and near-optimal approximations with optimality gaps
never exceeding 3% for $N=10$ and 0.5% for $N\geq 100$. Interestingly, the
optimal first-stage decisions are different, and this can be seen from their
out-of-sample performance that we present in the next subsection.

(a) $N=10$

(b) $N=100$

(c) $N=1000$

Figure 10: Optimality gaps using the continuous relaxation $\mathcal{Z}^{0}$
and the heuristically computed level-1 Lovász-Schrijver relaxation
$\tilde{\mathcal{Z}}^{1}$, as a function of $\nu$ and $N$ where
$\varepsilon=\nu\sqrt{N^{-1}\log(N+1)}$.

Figure 11 shows that both relaxations have small computation times ($<60$
seconds on average) across all $N$ and $\nu$. When compared with the column-
and-constraint generation scheme, the relative difference in their computation
times is minor for small sample sizes $N$ but increases significantly for
large sample sizes, where the column-and-constrain generation method can be
slower by more than a factor of 10. Similar to the Benders scheme, the mixed-
integer subproblems in the column-and-constraint generation scheme can cause
slow convergence, and roughly 3% of its runs did not terminate within 10
minutes.

(a) $N=10$

(b) $N=100$

(c) $N=1000$

Figure 11: Computation times using the continuous $\mathcal{Z}^{0}$ and
heuristically computed level-1 Lovász-Schrijver relaxation
$\tilde{\mathcal{Z}}^{1}$ for solving formulation (8), and using the column-
and-constraint generation scheme for solving formulation (7), as a function of
$\nu$ and $N$, where $\varepsilon=\nu\sqrt{N^{-1}\log(N+1)}$.

####  5.2.2 Out-of-sample performance and finite sample guarantee

Similar to Section 5.1.3, we estimate the out-of-sample performance of the
first-stage solutions of the lift-and-project and sample average
approximations for different sample sizes $N$ and Wasserstein radii
$\varepsilon=\nu\sqrt{N^{-1}\log(N+1)}$, where
$\nu\in\\{0,10^{-3},10^{-2},10^{-1}\\}$. The results are summarized in Figure
12.

First, Figure 12(a) reports the reliability of the level-1 Lovász-Schrijver
relaxation $\tilde{\mathcal{Z}}^{1}$, which is the empirical probability (over
the 100 sets of training samples) that its optimal value is an upper bound on
its out-of-sample cost. We observe that the reliability increases not only
with increasing values of $\nu$ (for fixed values of $N$) but also with
increasing values of $N$ (for fixed values of $\nu$). The choice $\nu=10^{-1}$
is reliable with probability 0.8 for training datasets of small size $N$ and
this increases to $>0.9$ for large values of $N$.

Figure 12(b) reports the relative improvement in out-of-sample cost of the
distributionally robust model over the sample average approximation, over the
100 independent sets of training samples. As before, we find that the
distributionally robust model consistently outperforms the sample average
approximation, particularly for small sample sizes $N$, where the magnitude of
the relative improvement can be as high as 7%, but decreases for large values
of $N$.

Finally, Figure 12(c) reports the improvement in the out-of-sample performance
of the level-1 Lovász-Schrijver relaxation $\tilde{\mathcal{Z}}^{1}$ for the
best choice of $\nu=10^{-1}$, relative to the continuous relaxation
$\mathcal{Z}^{0}$ for the same $\nu$. The relative improvement of the level-1
relaxation $\tilde{\mathcal{Z}}^{1}$ over the continuous relaxation
$\mathcal{Z}^{0}$ is small yet consistently non-negative. This is not
unexpected since the latter already has tight in-sample optimality gaps, as
can be seen from Figure 10. Nevertheless, we expect the relative improvements
to be instance-dependent similar to optimal power flow, and even a few
percentage points can result in long-term economic benefits.

(a) Reliability

(b) Out-of-sample performance

(c) Relative improvement

Figure 12: Reliability (left plot), and relative improvements in the out-of-
sample performance of the level-1 Lovász-Schrijver relaxation
$\tilde{\mathcal{Z}}^{1}$ when compared with the sample average approximation
(middle plot), and the continuous relaxation $\mathcal{Z}^{0}$ (right plot),
as a function of training sample size $N$.

We conclude this section by noting that [68] address problems where
$\mathcal{Q}(\bm{x},\bm{\xi})$ is the optimal value of a linear program, and
the ambiguity set $\mathcal{P}$ is the type-${\infty}$ Wasserstein ball which
is closely related to yet distinct from (3). They provide formulations that
can serve as an alternative means to generate first-stage decisions.
Interestingly, when applied to the multi-commodity network design problem,
they reduce to one of the following: (i) sample average approximation
($\varepsilon<1$) or deterministic problem under the worst-case realization
where all nodes have failed ($\varepsilon\geq 1$), when $d$ is induced by the
$\infty$-norm [68, Theorem 4], or (ii) classical robust optimization (19) with
$K=1$ when $d$ is induced by the $p$-norm with $p\in[1,\infty)$ and
$\varepsilon<\sqrt[p]{2}$ [68, Theorem 6].

##  6 Conclusions

Despite their ubiquity in real-world networks, optimization problems affected
by rare high-impact uncertainties have not received much attention. This is
partly because of the lack of available data given their rare nature, and
partly because of the incapability of classical sample average approximations
to address them effectively. This paper takes a step toward addressing these
limitations by motivating a distributionally robust approach to the problem
using Wasserstein ambiguity sets. Notably, we extend the state of the art in
data-driven optimization by encompassing not only two-stage conic problems but
also high-dimensional discrete uncertainties. By exploiting ideas from
nonlinear penalty methods and lift-and-project techniques in global
optimization, we present a simple, tractable, and tight approximation of the
problem that can be efficiently computed and iteratively improved. We use our
method to tackle optimal power flow problems with random transmission line
outages and multi-commodity network design problems with random node failures.
We find that the method can strongly outperform classical sample average and
robust optimization approaches, especially when failures are rare but can lead
to high costs associated with loss of electric power or commodity flows.

## Acknowledgments

This material is based upon work supported in part by the U.S. Department of
Energy, Office of Science and Office of Electricity Delivery & Energy
Reliability, Advanced Grid Research and Development (under contract number DE-
AC02-06CH11357), and in part by the Office of Naval Research (under Award
number N00014-18-1-2075).

## References

  * [1] Ongun Alsac and Brian Stott “Optimal load flow with steady-state security” In _IEEE Transactions on Power Apparatus and Systems_ IEEE, 1974, pp. 745–751
  * [2] Amir Ardestani-Jaafari and Erick Delage “Linearized Robust Counterparts of Two-stage Robust Optimization Problems with Applications in Operations Management” In _Available at Optimization Online_ , 2017
  * [3] Sogol Babaeinejadsarookolaee et al. “The power grid library for benchmarking ac optimal power flow algorithms” In _arXiv preprint arXiv:1908.02788_ , 2019
  * [4] Egon Balas, Sebastián Ceria and Gérard Cornuéjols “A lift-and-project cutting plane algorithm for mixed 0–1 programs” In _Mathematical Programming_ 58.1-3 Springer, 1993, pp. 295–324
  * [5] Manish Bansal, Kuo-Ling Huang and Sanjay Mehrotra “Decomposition algorithms for two-stage distributionally robust mixed binary programs” In _SIAM Journal on Optimization_ 28.3, 2018, pp. 2360–2383
  * [6] Javiera Barrera et al. “Chance-constrained problems and rare events: an importance sampling approach” In _Mathematical Programming_ 157.1 Springer, 2016, pp. 153–189
  * [7] Güzin Bayraksan and David K Love “Data-driven stochastic programming using phi-divergences” In _The Operations Research Revolution_ INFORMS, 2015, pp. 1–19
  * [8] Aharon Ben-Tal, Alexander Goryashko, Elana Guslitzer and Arkadi Nemirovski “Adjustable robust solutions of uncertain linear programs” In _Mathematical Programming_ 99.2, 2004, pp. 351–376
  * [9] Aharon Ben-Tal et al. “Robust solutions of optimization problems affected by uncertain probabilities” In _Management Science_ 59.2, 2013, pp. 341–357
  * [10] Bertrand Berche, Christian Von Ferber, Taras Holovatch and Yu Holovatch “Resilience of public transport networks against attacks” In _The European Physical Journal B_ 71.1, 2009, pp. 125–137
  * [11] Dimitris Bertsimas, Vishal Gupta and Nathan Kallus “Robust sample average approximation” In _Mathematical Programming_ 171.1-2 Springer, 2018, pp. 217–282
  * [12] Dimitris Bertsimas, Shimrit Shtern and Bradley Sturt “Two-stage sample robust optimization” In _arXiv preprint arXiv:1907.07142_ , 2019
  * [13] Daniel Bienstock and Abhinav Verma “The $N-k$ Problem in Power Grids: New Models, Formulations, and Numerical Experiments” In _SIAM Journal on Optimization_ 20.5, 2010, pp. 2352–2380
  * [14] Jose Blanchet, Yang Kang and Karthyek Murthy “Robust Wasserstein profile inference and applications to machine learning” In _Journal of Applied Probability_ 56.3, 2019, pp. 830–857
  * [15] Jose Blanchet and Karthyek Murthy “Quantifying distributional model risk via optimal transport” In _Mathematics of Operations Research_ 44.2, 2019, pp. 565–600
  * [16] Natashia Boland, Matteo Fischetti, Michele Monaci and Martin Savelsbergh “Proximity Benders: a decomposition heuristic for stochastic programs” In _Journal of Heuristics_ 22.2 Springer, 2016, pp. 181–198
  * [17] Amarjit Budhiraja, Shu Lu, Yang Yu and Quoc Tran-Dinh “Minimization of a Class of Rare Event Probabilities and Buffer Probabilities of Exceedance” In _arXiv preprint arXiv:1902.07829_ , 2019
  * [18] Samuel Burer and Hongbo Dong “Representing quadratically constrained quadratic programs as generalized copositive programs” In _Operations Research Letters_ 40.3 Elsevier, 2012, pp. 203–206
  * [19] Samuel Burer and Adam N Letchford “Non-convex mixed-integer nonlinear programming: A survey” In _Surveys in Operations Research and Management Science_ 17.2 Elsevier, 2012, pp. 97–106
  * [20] Giuseppe Calafiore and Marco C Campi “Uncertain convex programs: randomized solutions and confidence levels” In _Mathematical Programming_ 102.1 Springer, 2005, pp. 25–46
  * [21] M.T. Çezik and G. Iyengar “Cuts for mixed 0-1 conic programming” In _Mathematical Programming_ 104.1, 2005, pp. 179–202
  * [22] Naiyuan Chiang and Andreas Grothey “Solving security constrained optimal power flow problems by a structure exploiting interior point method” In _Optimization and Engineering_ 16.1 Springer, 2015, pp. 49–71
  * [23] Teodor Gabriel Crainic, Antonio Frangioni and Bernard Gendron “Bundle-based relaxation methods for multicommodity capacitated fixed charge network design” In _Discrete Applied Mathematics_ 112.1, 2001, pp. 73–99
  * [24] Etienne De Klerk “The complexity of optimizing over a simplex, hypercube or sphere: a short survey” In _Central European Journal of Operations Research_ 16.2 Springer, 2008, pp. 111–125
  * [25] Erick Delage and Yinyu Ye “Distributionally robust optimization under moment uncertainty with application to data-driven problems” In _Operations research_ 58.3, 2010, pp. 595–612
  * [26] Amir Dembo and Ofer Zeitouni “Large Deviations Techniques and Applications” Springer Berlin Heidelberg, 2010
  * [27] Bharat T Doshi et al. “Optical network design and restoration” In _Bell Labs Technical Journal_ 4.1, 1999, pp. 58–84
  * [28] Gabriele Eichfelder and Johannes Jahn “Set-semidefinite optimization” In _Journal of Convex Analysis_ 15.4, 2008, pp. 767–801
  * [29] Virginie Gabrel, Mathieu Lacroix, Cécile Murat and Nabila Remli “Robust location transportation problems under uncertain demands” In _Discrete Applied Mathematics_ 164 Elsevier, 2014, pp. 100–111
  * [30] Rui Gao and Anton J Kleywegt “Distributionally robust stochastic optimization with Wasserstein distance” In _arXiv preprint arXiv:1604.02199_ , 2016
  * [31] Michael R. Garey and David S. Johnson “Computers and Intractability: A Guide to the Theory of NP-Completeness” New York, NY, USA: W. H. Freeman & Co., 1979
  * [32] Angelos Georghiou, Daniel Kuhn and Wolfram Wiesemann “The decision rule approach to optimization under uncertainty: methodology and applications” In _Computational Management Science_ , 2018
  * [33] Fred Glover “Improved linear integer programming formulations of nonlinear integer problems” In _Management Science_ 22.4 INFORMS, 1975, pp. 455–460
  * [34] Chung Ha “A noncompact minimax theorem” In _Pacific Journal of Mathematics_ 97.1 Mathematical Sciences Publishers, 1981, pp. 115–117
  * [35] Grani A Hanasusanto and Daniel Kuhn “Conic programming reformulations of two-stage distributionally robust linear programs over Wasserstein balls” In _Operations Research_ 66.3, 2018, pp. 849–869
  * [36] Ruiwei Jiang and Yongpei Guan “Risk-Averse Two-Stage Stochastic Program with Distributional Ambiguity” In _Operations Research_ 66.5, 2018, pp. 1390–1405
  * [37] Ruiwei Jiang, Minseok Ryu and Guanglin Xu “Data-Driven Distributionally Robust Appointment Scheduling over Wasserstein Balls” In _arXiv preprint arXiv:1907.03219_ , 2019
  * [38] Burak Kocuk, Santanu S Dey and X Andy Sun “Strong SOCP relaxations for the optimal power flow problem” In _Operations Research_ 64.6 INFORMS, 2016, pp. 1177–1196
  * [39] Hiroshi Konno “Maximization of a convex quadratic function under linear constraints” In _Mathematical Programming_ 11.1 Springer, 1976, pp. 117–127
  * [40] Daniel Kuhn, Peyman Mohajerin Esfahani, Viet Anh Nguyen and Soroosh Shafieezadeh-Abadeh “Wasserstein distributionally robust optimization: Theory and applications in machine learning” In _Operations Research & Management Science in the Age of Analytics_ INFORMS, 2019, pp. 130–166
  * [41] Jean B Lasserre “Global optimization with polynomials and the problem of moments” In _SIAM Journal on optimization_ 11.3 SIAM, 2001, pp. 796–817
  * [42] Monique Laurent “A comparison of the Sherali-Adams, Lovász-Schrijver, and Lasserre relaxations for 0–1 programming” In _Mathematics of Operations Research_ 28.3 INFORMS, 2003, pp. 470–496
  * [43] László Lovász and Alexander Schrijver “Cones of matrices and set-functions and 0–1 optimization” In _SIAM Journal on Optimization_ 1.2 SIAM, 1991, pp. 166–190
  * [44] Steven H Low “Convex relaxation of optimal power flow—Part I: Formulations and equivalence” In _IEEE Transactions on Control of Network Systems_ 1.1 IEEE, 2014, pp. 15–27
  * [45] James Luedtke and Shabbir Ahmed “A sample approximation approach for optimization with probabilistic constraints” In _SIAM Journal on Optimization_ 19.2 SIAM, 2008, pp. 674–699
  * [46] Fengqiao Luo and Sanjay Mehrotra “Decomposition algorithm for distributionally robust optimization using Wasserstein metric with an application to a class of regression models” In _European Journal of Operational Research_ 278.1, 2019, pp. 20–35
  * [47] Athina Markopoulou et al. “Characterization of Failures in an Operational IP Backbone Network” In _IEEE/ACM Transactions on Networking_ 16.4, 2008, pp. 749–762
  * [48] Michel Minoux “Networks synthesis and optimum network design problems: Models, solution methods and applications” In _Networks_ 19.3 Wiley Online Library, 1989, pp. 313–360
  * [49] Areesh Mittal, Can Gokalp and Grani A Hanasusanto “Robust quadratic programming with mixed-integer uncertainty” In _INFORMS Journal on Computing_ 32.2 INFORMS, 2020, pp. 201–218
  * [50] Peyman Mohajerin Esfahani and Daniel Kuhn “Data-driven distributionally robust optimization using the Wasserstein metric: Performance guarantees and tractable reformulations” In _Mathematical Programming_ 171.1-2, 2018, pp. 115–166
  * [51] North American Electric Reliability Corporation “Transmission System Planning Performance Requirements” In _TPL-001-4_ , 2017
  * [52] Manfred Padberg “Approximating separable nonlinear functions via mixed zero-one programs” In _Operations Research Letters_ 27.1 Elsevier, 2000, pp. 1–5
  * [53] Krzysztof Postek, Dick Hertog and Bertrand Melenberg “Computationally tractable counterparts of distributionally robust constraints on risk measures” In _SIAM Review_ 58.4 SIAM, 2016, pp. 603–650
  * [54] Janez Povh and Franz Rendl “Copositive and semidefinite relaxations of the quadratic assignment problem” In _Discrete Optimization_ 6.3 Elsevier, 2009, pp. 231–241
  * [55] Pavel Praks, Vytis Kopustinskas and Marcelo Masera “Monte-Carlo-based reliability and vulnerability assessment of a natural gas transmission system due to random network component failures” In _Sustainable and Resilient Infrastructure_ 2.3, 2017, pp. 97–107
  * [56] Hamed Rahimian, Güzin Bayraksan and Tito Homem-de-Mello “Identifying effective scenarios in distributionally robust stochastic programs with total variation distance” In _Mathematical Programming_ 173.1-2, 2019, pp. 393–430
  * [57] Hamed Rahimian and Sanjay Mehrotra “Distributionally robust optimization: A review” In _arXiv preprint arXiv:1908.05659_ , 2019
  * [58] Alexander Shapiro “Tutorial on risk neutral, distributionally robust and risk averse multistage stochastic programming” In _Available at Optimization Online_ , 2018
  * [59] Alexander Shapiro, Darinka Dentcheva and Andrzej Ruszczyński “Lectures on stochastic programming: modeling and theory” SIAM, 2009
  * [60] Hanif D Sherali and Warren P Adams “A hierarchy of relaxations between the continuous and convex hull representations for zero-one programming problems” In _SIAM Journal on Discrete Mathematics_ 3.3 SIAM, 1990, pp. 411–430
  * [61] Robert A. Stubbs and Sanjay Mehrotra “A branch-and-cut method for 0-1 mixed convex programming” In _Mathematical Programming_ 86, 1999, pp. 515–532
  * [62] Aurélie Thiele, Tara Terry and Marina Epelman “Robust linear optimization with recourse”, 2009
  * [63] U.S. Department of Energy Advanced Research Projects Agency-Energy “SCOPF Problem Formulation: Challenge 1” In _Grid Optimization Competition_ , 2019
  * [64] U.S. Energy Information Administration “Annual Summaries” See also http://insideenergy.org/2014/08/18/data-explore-15-years-of-power-outages/ In _Electric Disturbance Events (OE-417)_ , 2017
  * [65] J. Vielma “Mixed Integer Linear Programming Formulation Techniques” In _SIAM Review_ 57.1, 2015, pp. 3–57
  * [66] Wolfram Wiesemann, Daniel Kuhn and Melvyn Sim “Distributionally robust convex optimization” In _Operations Research_ 62.6, 2014, pp. 1358–1376
  * [67] David Wozabal “A framework for optimization under ambiguity” In _Annals of Operations Research_ 193.1, 2012, pp. 21–47
  * [68] Weijun Xie “Tractable reformulations of two-stage distributionally robust linear programs over the type-$\infty$ Wasserstein ball” In _Operations Research Letters_ 48.4, 2020, pp. 513–523
  * [69] Guanglin Xu and Samuel Burer “A copositive approach for two-stage adjustable robust optimization with uncertain right-hand sides” In _Computational Optimization and Applications_ 70.1, 2018, pp. 33–59
  * [70] 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
  * [71] Chaoyue Zhao “Data-driven risk-averse stochastic program and renewable energy integration”, 2014
  * [72] Chaoyue Zhao and Yongpei Guan “Data-driven risk-averse stochastic optimization with Wasserstein metric” In _Operations Research Letters_ 46.2 Elsevier, 2018, pp. 262–267
  * [73] Ray D Zimmerman and Carlos E Murillo-Sánchez “MATPOWER 6.0 user’s manual” In _PSERC: Tempe, AZ, USA_ , 2016

##  Appendix A Extended discussion of Example 1

In Section A.1, we analyze the optimal solution of the network optimization
problem in Example 1 under the true failure distribution. In Section A.2, we
analyze the solutions of the sample average and distributionally robust
formulations when the empirical distribution puts all its weight on the
realization $\bm{\xi}=\bm{0}$, where none of the nodes fail.

###  A.1 Optimal solution under the true distribution

Problem symmetry implies that the optimal arc capacities have the structure
shown in Figure 13. Table 2 presents a calculation of the second-stage loss
function $\mathcal{Q}(\bm{x},\bm{\xi})$ for each $\bm{\xi}\in\\{0,1\\}^{3}$.
The optimal arc capacities can then be determined by solving the two-
dimensional piecewise linear convex optimization problem (note that
$0<\epsilon\ll 1$ can be any small positive constant):

$\mathop{\text{minimize}}_{\bm{x}\in\mathbb{R}^{2}_{+}}\;4x_{1}+(4-\epsilon)x_{2}+\sum_{\bm{\xi}\in\\{0,1\\}^{3}}\mathbb{P}[\tilde{\bm{\xi}}=\bm{\xi}]\mathcal{Q}(\bm{x},\bm{\xi})$

A straightforward calculation shows that the objective function minimized at
$(x_{1},x_{2})=(100,100)$.

121$0.25\%$2$1\%$3$0.25\%$67100100100100$x_{1}$$x_{2}$$x_{2}$$x_{1}$$x_{1}$$x_{2}$$x_{2}$$x_{1}$ Figure 13: The structure of the optimal solution, with arc capacities $x_{1}$ and $x_{2}$ indicated above each arc. Table 2: The second-stage loss function $\mathcal{Q}(\bm{x},\bm{\xi})$ for each $\bm{\xi}\in\\{0,1\\}^{3}$. Here, we use $q=1000$ to denote the penalty cost per unit of supply shortfall, and we use $[\cdot]_{+}\coloneqq\max\\{\cdot,0\\}$. $\bm{\xi}=(\xi_{1},\xi_{2},\xi_{3})$ | $\lVert\bm{\xi}\rVert_{1}$ | $\mathbb{P}[\tilde{\bm{\xi}}=\bm{\xi}]$ |  $\mathcal{Q}(\bm{x},\bm{\xi})$
---|---|---|---
$(1,1,1)$ | $3$ | $6.3\times 10^{-8}$ |  $200q$
$(1,0,1)$ | $2$ | $6.2\times 10^{-6}$ |  $2[100-x_{2}]_{+}q$
$(1,1,0)$ | $2$ | $2.5\times 10^{-5}$ |  $\left(100+[100-x_{1}]_{+}\right)q$
$(0,1,1)$ | $2$ | $2.5\times 10^{-5}$ |  $\left(100+[100-x_{1}]_{+}\right)q$
$(1,0,0)$ | $1$ | $2.5\times 10^{-3}$ |  $\left(100-x_{2}+[100-x_{1}-x_{2}]_{+}\right)q$
$(0,0,1)$ | $1$ | $2.5\times 10^{-3}$ |  $\left(100-x_{2}+[100-x_{1}-x_{2}]_{+}\right)q$
$(0,1,0)$ | $1$ | $1.0\times 10^{-2}$ |  $2[100-x_{1}]_{+}q$
$(0,0,0)$ | $0$ | $9.6\times 10^{-1}$ |  $2\left[100-x_{1}-x_{2}\right]_{+}q$

###  A.2 Optimal solutions using the sample average and distributionally
robust formulations

Suppose now that the empirical distribution
$\hat{\mathbb{P}}_{N}=\delta_{\bm{0}}$; that is, it puts all its weight on the
realization $\bm{\xi}=(0,0,0)$, where none of the nodes fail. The problem
symmetry carries forth to the sample average approximation, and its optimal
solution has the same symmetry as Figure 13. The solution can be determined by
solving the following problem:

$\mathop{\text{minimize}}_{\bm{x}\in\mathbb{R}^{2}_{+}}\;4x_{1}+(4-\epsilon)x_{2}+\mathcal{Q}(\bm{x},\bm{0})\;=\;\mathop{\text{minimize}}_{\bm{x}\in\mathbb{R}^{2}_{+}}\;4x_{1}+(4-\epsilon)x_{2}+2\left[100-x_{1}-x_{2}\right]_{+}q$

Since $\epsilon>0$, the objective function is minimized at
$(x_{1},x_{2})=(0,100)$; this is plotted in Figure 2(a).

Consider now the distributionally robust formulation (2) using a Wasserstein
ambiguity set (3) of radius $\varepsilon>0$, and defined with respect to the
metric
$d(\bm{\xi},\bm{\xi^{\prime}})=\lVert\bm{\xi}-\bm{\xi^{\prime}}\rVert_{1}$.
Again, problem symmetry leads to the same optimal solution structure as Figure
13. The optimal solution can be determined by solving formulation (7):

$\mathop{\text{minimize}}_{(x_{1},x_{2},\alpha)\in\mathbb{R}^{3}_{+}}\;4x_{1}+(4-\epsilon)x_{2}+\alpha\varepsilon+\max_{\bm{\xi}\in\\{0,1\\}^{3}}\left\\{\mathcal{Q}(\bm{x},\bm{\xi})-\alpha\lVert\bm{\xi}\rVert_{1}\right\\}$

For any $\varepsilon\in(0,3)$, the minimizer is
$(x_{1},x_{2},\alpha)=(\frac{2}{3}100,\frac{1}{3}100,\frac{2}{3}100q)$; this
is plotted in Figure 3(a).

If we change the definition of the metric to be
$d(\bm{\xi}^{\prime},\bm{\xi}^{\prime\prime})=1$ whenever
$\bm{\xi}^{\prime}\neq\bm{\xi}^{\prime\prime}$ and $0$ otherwise, then the
ambiguity set is equivalent to the $\phi$-divergence ambiguity set based on
the total variation distance. The corresponding optimal solution can be
determined by solving formulation (7):

$\mathop{\text{minimize}}_{(x_{1},x_{2},\alpha)\in\mathbb{R}^{3}_{+}}\;4x_{1}+(4-\epsilon)x_{2}+\alpha\varepsilon+\max\left\\{\mathcal{Q}(\bm{x},\bm{0}),\max_{\bm{\xi}\in\\{0,1\\}^{3}\setminus\\{\bm{0}\\}}\mathcal{Q}(\bm{x},\bm{\xi})-\alpha\right\\}$

The innermost maximum is attained at $\bm{\xi}=\bm{\mathrm{e}}=(1,1,1)$
irrespective of $\varepsilon,\bm{x},\alpha$. Therefore, the outer maximum
becomes
$\max\left\\{\mathcal{Q}(\bm{x},\bm{0}),\mathcal{Q}(\bm{x},\bm{\mathrm{e}})-\alpha\right\\}$
and it is minimized at its break-point where
$\mathcal{Q}(\bm{x},\bm{0})=\mathcal{Q}(\bm{x},\bm{\mathrm{e}})-\alpha$. The
problem then becomes equivalent to the sample average approximation

$\mathop{\text{minimize}}_{(x_{1},x_{2},\alpha)\in\mathbb{R}^{3}_{+}}\;4x_{1}+(4-\epsilon)x_{2}+\alpha\varepsilon+\mathcal{Q}(\bm{x},\bm{0})$

and its minimum is attained at $(x_{1},x_{2},\alpha)=(0,100,0)$ for any
$\varepsilon\geq 0$.

##  Appendix B Complexity analysis

The feasible region of the inner optimization problem in Theorem 2 is the
convex hull of an MICP-representable set. This allows us to exploit existing
tools on characterizing the convex hulls of MICP sets. Section B.1 highlights
several problem instances (i.e., sufficient conditions) where this convex hull
is efficiently computable and hence the distributionally robust two-stage
problem (2) is tractable. In general, however, the presence of discrete random
parameters makes the problem intractable even when the second-stage problem
$\mathcal{Q}(\bm{x},\bm{\xi})$ is benign. This is proved in Section B.2.

###  B.1 Tractable cases

We preface this subsection by noting that any notion of tractability of the
distributionally robust two-stage problem (2) requires that optimizing
$c(\bm{x})$ over $\bm{x}\in\mathcal{X}$ can be done in a tractable manner. We
therefore assume this to be the case throughout this subsection.

Suppose that $\Xi=\\{\bm{v}^{(1)},\ldots,\bm{v}^{(K)}\\}$ has a known inner
description, where the number of vertices $K$ of $\Xi$ grows polynomially with
the dimension $M$. For example, this is the case of budget supports of the
form $\Xi=\\{\bm{\xi}\in\\{0,1\\}^{M}:\bm{e}^{\top}\bm{\xi}\leq k\\}$, where
$k$ is a small, fixed input (say, $\leq 3$) and $K=O(M^{k})$. Now, if
$\bm{\xi}$ is fixed to one of $\bm{v}^{(k)}$, $k\in[K]$, then the set
$\mathcal{Z}_{i}$ in (9d) becomes a closed convex set. In other words,
$\mathcal{Z}_{i}$ is the union of $K$ closed convex sets, and its convex hull
can be described compactly.

######  Proposition 2.

Assume $\Xi=\\{\bm{v}^{(1)},\ldots,\bm{v}^{(K)}\\}$, where $K$ is a polynomial
function of $M$. Then, the distributionally robust two-stage problem (2) is
equivalent to the following tractable problem:

$\begin{array}[]{l@{\;\;}l}\displaystyle\mathop{\text{minimize}}&\displaystyle
c(\bm{x})+\alpha\varepsilon+\frac{1}{N}\sum_{i=1}^{N}\sigma_{i}\\\
\text{subject to}&\bm{x}\in\mathcal{X},\;\alpha\geq 0\\\
&\begin{rcases}\sigma_{i}\in\mathbb{R},\;\bm{y}^{(k)}\in\mathcal{Y}\\\
\sigma_{i}\geq\bm{q}(\bm{v}^{(k)})^{\top}\bm{y}^{(k)}-\alpha
d(\bm{v}^{(k)},\hat{\bm{\xi}}^{(i)})\\\
\bm{W}(\bm{v}^{(k)})\bm{y}^{(k)}\geq\bm{T}(\bm{x})\bm{v}^{(k)}+\bm{h}(\bm{x})\end{rcases}\;\forall
k\in[K],i\in[N].\end{array}$

###### Proof.

From Lemma 1 and the stated hypothesis, problem (2) is equivalent to

$\begin{array}[]{l@{\;\;\;}l}\phantom{=}\displaystyle\mathop{\text{minimize}}_{\bm{x}\in\mathcal{X},\alpha\geq
0}&\displaystyle
c(\bm{x})+\alpha\varepsilon+\frac{1}{N}\sum_{i=1}^{N}\max_{k\in[K]}\left\\{\mathcal{Q}(\bm{x},\bm{v}^{(k)})-\alpha
d(\bm{v}^{(k)},\hat{\bm{\xi}}^{(i)})\right\\}.\end{array}$

We can then introduce the epigraphical variables $\sigma_{i}$, $i\in[N]$ to
model the $i$th inner maximization in the above sum. The stated reformulation
then follows after introducing second-stage variables $\bm{y}^{(k)}$ to
represent a minimizer of $\mathcal{Q}(\bm{x},\bm{v}^{(k)})$, $k\in[K]$. ∎

Suppose now that we have an outer description of
$\Xi=\\{\bm{\xi}\in\mathbb{Z}^{M}:\bm{E}\bm{\xi}\leq\bm{f}\\}$ and the second-
stage problem $\mathcal{Q}(\bm{x},\bm{\xi})$ is a linear program
$\mathcal{Y}=\mathbb{R}^{N_{2}}_{+}$ with objective uncertainty so that
Corollary 1 is applicable. The next observation relies on the concept of an
ideal formulation [52, 65]. A mixed-integer linear set
$\\{\bm{x}\in\mathbb{Z}^{p}\times\mathbb{R}^{n-p}:\bm{A}\bm{x}\leq\bm{b}\\}$
is said to be ideal if its convex hull has at least one vertex and it is equal
to its linear relaxation. We show that the convex hull reformulation (8)
becomes tractable if the set $\mathcal{Z}_{i}$ defined in (10b) is mixed-
integer linear and ideal.

######  Proposition 3.

Suppose that $\Xi=\\{\xi\in\mathbb{Z}^{M}:\bm{E}\bm{\xi}\leq\bm{f}\\}$ and
$\mathcal{Q}(\bm{x},\bm{\xi})$ is the optimal value of a linear program with
uncertain objective coefficients: $\bm{W}(\bm{\xi})=\bm{W}_{0}$ and
$\bm{T}(\bm{x})=\bm{0}$. Suppose also that the mixed-integer linear
formulation
$\\{(\bm{\xi},\bm{\lambda})\in\mathbb{Z}^{M}\times\mathbb{R}^{L}_{+}:\bm{E}\bm{\xi}\leq\bm{f},\;\bm{W}_{0}^{\top}\bm{\lambda}-\bm{Q}\bm{\xi}\leq\bm{q}_{0}\\}$
is ideal.

  1. (i)

If the reference metric
$d(\bm{\xi},\bm{\xi^{\prime}})=\lVert\bm{\xi}-\bm{\xi^{\prime}}\rVert_{1}$ is
the 1-norm, then the two-stage problem (2) is equivalent to the following
tractable optimization problem:

$\begin{array}[]{l@{\;\;}l}\displaystyle\mathop{\text{minimize}}&\displaystyle
c(\bm{x})+\alpha\varepsilon+\frac{1}{N}\sum_{i=1}^{N}\left[\bm{q}_{0}^{\top}\bm{y}^{(i)}+\bm{f}^{\top}\bm{\eta}^{(i)}-\alpha\bm{\mathrm{e}}^{\top}\hat{\bm{\xi}}^{(i)}\right]\\\
\text{subject to}&\bm{x}\in\mathcal{X},\;\alpha\geq 0\\\
&\begin{rcases}\bm{y}^{(i)}\in\mathcal{Y},\;\bm{\eta}^{(i)}\in\mathbb{R}^{F}_{+}\\\
-\bm{Q}^{\top}\bm{y}^{(i)}+\bm{E}^{\top}\bm{\eta}^{(i)}\geq-\alpha(\bm{\mathrm{e}}-2\hat{\bm{\xi}}^{(i)})\\\
\bm{W}_{0}\bm{y}^{(i)}\geq\bm{h}(\bm{x})\end{rcases}\;\forall
i\in[N].\end{array}$

  2. (ii)

If
$\varepsilon\geq\max_{\bm{\xi},\bm{\xi}^{\prime}\in\Xi}d(\bm{\xi},\bm{\xi}^{\prime})$,
then the two-stage (classical robust optimization) problem (2) is equivalent
to the following tractable optimization problem:

$\begin{array}[]{l@{\;\;}l}\displaystyle\mathop{\text{minimize}}&\displaystyle
c(\bm{x})+\bm{q}_{0}^{\top}\bm{y}+\bm{f}^{\top}\bm{\eta}\\\ \text{subject
to}&\displaystyle\bm{x}\in\mathcal{X},\;\bm{y}\in\mathcal{Y},\;\bm{\eta}\in\mathbb{R}^{F}_{+}\\\
&\displaystyle-\bm{Q}^{\top}\bm{y}+\bm{E}^{\top}\bm{\eta}\geq\bm{0}\\\
&\displaystyle\bm{W}_{0}\bm{y}\geq\bm{h}(\bm{x})\end{array}$

###### Proof.

The conditions of this proposition allow us to use the result of Corollary 1.

(i) When the reference metric
$d(\bm{\xi},\bm{\xi^{\prime}})=\lVert\bm{\xi}-\bm{\xi^{\prime}}\rVert_{1}$ is
the 1-norm, the condition
$(\bm{\xi}-\hat{\bm{\xi}}^{(i)},\tau)\in\mathcal{C}^{M+1}$ appearing in (10b)
is equivalent to $\tau\geq\lVert\bm{\xi}-\hat{\bm{\xi}}^{(i)}\rVert_{1}$. If
we exploit the fact that $\bm{\xi}$ and $\hat{\bm{\xi}}^{(i)}$ are both
binary-valued, then the right-hand side of this inequality is equal to
$(\bm{\mathrm{e}}-2\hat{\bm{\xi}}^{(i)})^{\top}\bm{\xi}+\bm{\mathrm{e}}^{\top}\hat{\bm{\xi}}^{(i)}$
and hence, linear in $\bm{\xi}$. Since the objective function of (10a) is
linear in $\tau$ and $\tau$ only appears in a single constraint in
$\mathcal{Z}_{i}$ that is linear in $\bm{\xi}$, we can reformulate (10a)–(10b)
as follows:

$\displaystyle\displaystyle
Z_{i}(\bm{x},\alpha)=\mathop{\text{maximize}}_{(\bm{\xi},\bm{\lambda})\in\mathop{\mathrm{cl\,conv}}(\hat{\mathcal{Z}})}\displaystyle\left\\{\bm{h}(\bm{x})^{\top}\bm{\lambda}-\alpha\left((\bm{\mathrm{e}}-2\hat{\bm{\xi}}^{(i)})^{\top}\bm{\xi}+\bm{\mathrm{e}}^{\top}\hat{\bm{\xi}}^{(i)}\right)\right\\},$
$\displaystyle\displaystyle\hat{\mathcal{Z}}=\left\\{(\bm{\xi},\bm{\lambda})\in\mathbb{Z}^{M}\times\mathbb{R}^{L}_{+}:\displaystyle\bm{E}\bm{\xi}\leq\bm{d},\;\bm{q}_{0}+\bm{Q}\bm{\xi}-\bm{W}_{0}^{\top}\bm{\lambda}\geq\bm{0}\right\\},$

where $\Xi=\\{\xi\in\mathbb{Z}^{M}:\bm{E}\bm{\xi}\leq\bm{d}\\}$ and
$\mathcal{Y}=\mathbb{R}^{N_{2}}_{+}$. By hypothesis, this formulation of
$\hat{\mathcal{Z}}$ is ideal, and hence
$\mathop{\mathrm{cl\,conv}}(\hat{\mathcal{Z}})=\\{(\bm{\xi},\bm{\lambda})\in\mathbb{R}^{M}\times\mathbb{R}^{L}_{+}:\bm{E}\bm{\xi}\leq\bm{f},\;\bm{W}_{0}^{\top}\bm{\lambda}-\bm{Q}\bm{\xi}\leq\bm{q}_{0}\\}$.
We then obtain the stated result by exploiting strong linear programming
duality.

(ii) When
$\varepsilon\geq\max_{\bm{\xi},\bm{\xi}^{\prime}\in\Xi}d(\bm{\xi},\bm{\xi}^{\prime})$,
problem (2) reduces to a classical two-stage robust optimization as per Remark
1. In this case, the optimal value of $\alpha$ in the convex hull
reformulation (8) is $0$ and the two-stage problem (2) is equivalent to

$\displaystyle\mathop{\text{minimize}}_{\bm{x}\in\mathcal{X}}\;\;\displaystyle
c(\bm{x})+\hat{Z}(\bm{x}),$

where the set $\hat{\mathcal{Z}}$ is defined as above and we define the
function $\hat{Z}:\mathcal{X}\mapsto\mathbb{R}$ as follows:

$\displaystyle\displaystyle\hat{Z}(\bm{x})=\mathop{\text{maximize}}_{(\bm{\xi},\bm{\lambda})\in\mathop{\mathrm{cl\,conv}}(\hat{\mathcal{Z}})}\displaystyle\bm{h}_{0}(\bm{x})^{\top}\bm{\lambda}.$

By the same argument as in part (i), the above formulation of
$\hat{\mathcal{Z}}$ is ideal, and we obtain the stated result using strong
linear programming duality. ∎

A well-known sufficient condition that guarantees idealness of a mixed-integer
linear formulation is total unimodularity of the constraint matrices. In
particular, if $\bm{W}_{0}$, $\bm{Q}$ and $\bm{E}$ are totally unimodular
(e.g., they are network matrices), then the conditions of Proposition 3 can be
satisfied. We refer to [65] for a more general overview of ideal formulations.

###  B.2 Computational complexity

Even though the preceding subsection presented some sufficient conditions for