# Tight and Compact Sample Average Approximation for Joint Chance-constrained
Problems with Applications to Optimal Power Flow

Álvaro Porras  Concepción Domínguez  Juan Miguel Morales<EMAIL_ADDRESS>Salvador Pineda  Optimization and Analytics for Sustainable energY Systems
(OASYS), University of Málaga, 29071 Málaga (Spain)

###### Abstract

In this paper, we tackle the resolution of chance-constrained problems
reformulated via Sample Average Approximation. The resulting data-driven
deterministic reformulation takes the form of a large-scale mixed-integer
program cursed with Big-Ms. We introduce an exact resolution method for the
MIP that combines the addition of a set of valid inequalities to tighten the
linear relaxation bound with coefficient strengthening and constraint
screening algorithms to improve its Big-Ms and considerably reduce its size.
The proposed valid inequalities are based on the notion of _k-envelopes_ , can
be computed offline using polynomial-time algorithms, and added to the MIP
program all at once. Furthermore, they are equally useful to boost the
strengthening of the Big-Ms and the screening rate of superfluous constraints.
We apply our procedures to a probabilistically-constrained version of the DC
Optimal Power Flow problem with uncertain demand. The chance constraint
requires that the probability of violating _any_ of the power system’s
constraints be lower than some parameter $\epsilon>0$. In a series of
numerical experiments which involve five power systems of different size, we
show the efficiency of the proposed methodology and compare it with some of
the best-performing convex inner approximations currently available in the
literature.

###### keywords:

Chance Constraints , Probabilistic Constraints , $k$-Violation Problems ,
Combinatorial Optimization , Mixed-Integer Programming , Valid Inequalities ,
Optimal Power Flow

††journal: arXiv

##  1 Introduction

Chance-constrained programming suits applications in areas where decisions
have to be made dealing with random parameters (Miller and Wagner (1965)). In
these situations, it is desirable to ensure feasibility of the system almost
surely, but there is hardly any decision which would guarantee it under
extreme events or unexpected random circumstances. Within this context,
chance-constrained programming offers a powerful modeling framework to
identify cost-efficient decisions that satisfy the problem’s constraints with
a high level of confidence.

In this paper, we focus on chance-constrained problems (CCPs) with joint
linear chance constraints and random RHS and LHS, i.e.,

$\displaystyle\min_{x}\quad$ $\displaystyle f(x)$ (1a) s.t. $\displaystyle
x\in X$ (1b) $\displaystyle\mathbb{P}\left\\{a_{j}(\omega)^{\top}x\leq
b_{j}(\omega),\ \forall j\in\mathcal{J}\right\\}\geq 1-\epsilon.$ (1c)

In (1), $x\in\mathbb{R}^{|\mathcal{I}|}$ is a vector of continuous decision
variables, $X\subseteq\mathbb{R}^{|\mathcal{I}|}$ represents a set of
deterministic constraints, and
$f:\mathbb{R}^{|\mathcal{I}|}\longrightarrow\mathbb{R}$ is a convex function.
Uncertainty is represented through the random vector $\omega$ taking values in
$\mathbb{R}^{d}$ and giving rise to a technology matrix with random rows
$a_{j}(\omega)\in\mathbb{R}^{|\mathcal{I}|}$, $j\in\mathcal{J}$ and random
$b_{j}(\omega)\in\mathbb{R}$, $j\in\mathcal{J}$. $\mathbb{P}$ is a probability
measure, and $\epsilon$ is a confidence or risk parameter, typically near
zero, so that the set of constraints (1c) are satisfied with probability at
least $(1-\epsilon)$. Apart from power systems, applications of CCPs include
supply chain, location and logistics (Taleizadeh et al. (2012), Shaw et al.
(2016), Elçi et al. (2018)), risk control in finance (Daníelsson et al.
(2008), Natarajan et al. (2008)), and healthcare problems such as operating
room planning (Najjarbashi and Lim (2020)) or vaccine allocation (Tanner et
al. (2008)), among others.

When the probabilistic constraint corresponds to (1c), the CCP has joint
chance constraints (JCC) and is hence classified as a joint CCP (JCCP), in
contrast with single CCPs (SCCPs), i.e. CCPs with individual or single chance
constraints (SCC) of the form
$\mathbb{P}\left\\{a_{j}(\omega)^{\top}x-b_{j}(\omega)\leq 0\right\\}\geq
1-\epsilon_{j}$, $\forall j\in\mathcal{J}$. JCCPs are suitable for contexts
where all constraints need to be simultaneously satisfied with a high
probability, and the dependence between random variables makes them clearly
harder. Both SCCPs and JCCPs have been extensively studied (see Prékopa
(2003), Van Ackooij et al. (2011) and the references therein).

There are a number of reasons why general CCPs are challenging. The first one
is the non-convexity of the feasible set. In general, the feasible region of a
CCP is not convex in the original space even when $x$ is continuous, there is
only RHS uncertainty and the constraints inside the probability in (1c) define
a polyhedral region (Küçükyavuz and Jiang (2022)). To circumvent this problem,
several approaches have been proposed. Some methods (e.g. Lagoa et al. (2005),
Henrion (2007), Henrion and Strugarek (2011)) give convexity results and
investigate the conditions under which the feasible region of problem (1) is
convex. In another line of research, various convex approximation schemes such
as quadratic (Ben-Tal and Nemirovski (2000)) or Bernstein approximation
(Nemirovski and Shapiro (2007)), have been proposed in the literature. The
CVaR approximation has gained a lot of popularity since its introduction
(Rockafellar and Uryasev (2000), Sun et al. (2014)), and remains one of the
most used methods to deal with stochastic problems. Nonetheless, the solutions
to the approximated problems err on the side of over-conservatism. In this
context, some iterative schemes such as ALSO-X have been recently proposed to
identify tighter inner convex approximations of the CCP at the expense of a
higher computational cost (Ahmed et al. (2017), Jiang and Xie (2022)).
Finally, other works suggest convex approximations for non-linear CCPs. For
instance, Hong et al. (2011) propose to solve the JCCP by a sequence of convex
approximations followed by a gradient-based Monte Carlo method, whereas Peña-
Ordieres et al. (2020) introduce a smooth sampling-based approximation.

The second difficulty of CCPs is that checking the feasibility of a given
solution is not, in general, an easy task. For instance, even if the
uncertainty follows a known continuous distribution, calculating the joint
probability requires a multi-dimensional integration, which becomes
increasingly difficult with the dimension of the random vector $\omega$. On
top of that, in most cases the distribution $\mathbb{P}$ is not fully known.

A commonly used approach to cope with these obstacles is to replace
$\mathbb{P}$ with its empirical estimate based on a set of random samples
(Shapiro (2003)). This approach is known in the technical literature as Sample
Average Approximation (SAA) and may also be seen as a data-driven approach
that works with observations of $\omega$ that are available to the decision-
maker even if the true data-generating distribution $\mathbb{P}$ is unknown.
The main advantage of SAA is that it allows to reformulate CCPs as
deterministic optimization problems that can be solved using standard
optimization techniques.

Indeed, the deterministic reformulation of (1) using SAA results in a mixed-
integer problem (MIP). For the resolution of CCPs using MIP reformulations, we
refer the reader to Küçükyavuz and Jiang (2022). When there is only RHS
uncertainty (i.e. the technology matrix is fixed), a reformulation of the
problem leads to a MIP with a set of constraints that form a _mixing set_ and
that has been extensively studied, alone or in combination with the knapsack
constraint that also appears in the formulation (Günlük and Pochet (2001),
Luedtke et al. (2010), Abdi and Fukasawa (2016)). Alternative reformulations
like the ones proposed by Dentcheva et al. (2000), Nair and Miller-Hooks
(2011) rely on the concept of $(1-\epsilon)$-efficient points. The case when
the technology matrix is random, while the RHS is not, has been studied e.g.
in Tayur et al. (1995). As for the general case, it has also been addressed in
the literature. Specifically, a large line of research has focused on the
development of _quantile cuts_ , a particular type of valid inequality that
can be viewed as a projection of a set of mixing inequalities for the MIP onto
the original problem space. These cuts and the associated quantile closure
have been recently studied in Qiu et al. (2014), Xie and Ahmed (2016, 2018),
Ahmed and Xie (2018) and successfully applied to computational experiments of
CCPs in Song et al. (2014), Ahmed et al. (2017), among others.

Regarding its disadvantages, the SAA solution tends to be optimistic and
underestimate the probability of violating the constraints (Nemirovski and
Shapiro (2006)). In fact, the performance of SAA is directly contingent on the
number of samples available. We refer the reader to the article by Luedtke and
Ahmed (2008), where relations are established between the empirical acceptable
probability of violation and the number of samples such that the SAA-based
solution be feasible in the JCCP with a predefined confidence level. In this
paper, however, we do not investigate how well the SAA-based solution
approximates the solution of the original JCCP. Statistical considerations
apart, the main focus of this paper is on proposing an efficient methodology
to solve the deterministic optimization model that results from reformulating
JCCPs using SAA. In particular, the method we propose solves the MIP to
optimality, and is based on the combination of a tightening-and-screening
procedure with the development of valid inequalities to obtain a formulation
which is compact and tight at a time.

We begin with a description of an iterative algorithm to strengthen the Big-Ms
present in the mixed-integer reformulation of the JCC. Interestingly, although
the procedure is not new (see Qiu et al. (2014)), we complement it with a
screening procedure that allows us to eliminate inequalities of the MIP.

Next, we introduce a new set of valid inequalities designed to strengthen the
linear relaxation of the model. As for their structure, each valid inequality
is developed using quantile information and avoiding the use of Big-M
constants, which are known to lead to weak linear relaxation bounds.
Furthermore, from a theoretical point of view, they are closely related to the
aforementioned quantile cuts. In fact, as proven in Appendix A, the addition
of our inequalities yields a feasible region equivalent to the quantile
closure of a specific relaxation of the MIP problem. The main advantage of our
inequalities is that, unlike quantile cuts, they are not NP-hard to compute,
and neither do they require a specific separation algorithm to include them
dynamically (that is, they can all be included in the model from the outset).
As many other techniques developed for JCCPs, our cuts also apply to the
SCCPs. Finally, our valid inequalities extend the results introduced in Roos
and Widmayer (1994) for the $k$-_violation problem_ , which can be seen as a
mixed-integer linear problem (MILP) reformulation of a SCCP. For ease of
explanation, the valid inequalities are proposed for a particular class of
linear chance constraints, but our results are also applicable to more general
types of CCPs, as we detail in Appendix B.

We apply our resolution method to the Optimal Power Flow (OPF) problem, for
which we propose a chance-constrained formulation and extensive computational
experiments featuring results for five standard power systems available in the
literature. The combination of the valid inequalities with the tightening of
the Big-Ms and the screening procedure allows us to effectively solve to
optimality instances that are not solved with the initial MIP formulation,
since the combination of both techniques notably reduces their size and
difficulty. We also compare our resolution approach with state-of-the-art
convex inner approximations of CCPs, in particular, the CVaR-based
approximation, ALSO-X, and ALSO-X+ (Jiang and Xie (2022)). Finally, we compare
our results with an exact approach for solving CCPs assuming SAA, namely, the
branch-and-cut decomposition algorithm proposed in Luedtke (2014). This
algorithm is based on the addition of quantile cuts strengthened using the
star inequalities of Atamtürk et al. (2000) to a relaxation of the original
problem.

The remainder of the paper is organized as follows. Section 2 involves the
reformulation of the problem into a MIP using SAA and the proposed
methodology: Subsection 2.1 describes the tightening and screening procedures,
whereas Subsection 2.2 introduces the valid inequalities and the necessary
algorithms to compute them. Section 3 is devoted to the application of the
results to the OPF problem. The section begins with an introduction to the
problem and the associated literature review, and Subsection 3.1 formalizes
the notation and the mathematical model solved in the computational study. In
Subsection 3.2 we present the case study, testing our results to solve
instances of the DC OPF available in the literature. Section 4 points further
research topics and includes some concluding remarks. Finally, the main text
is accompanied by two appendices: Appendix A establishes the relationship
between our valid inequalities and the quantile cuts present in the
literature, and Appendix B details the complexity of the algorithms proposed
in Subsection 2.2 and discusses the generalization of our methodology and its
possible application to other types of CCPs.

##  2 Solving JCCPs via Sample Average Approximation

Using SAA, problem (1) can be reformulated into a MIP. To this aim, assume
that $\omega$ has a finite discrete support defined by a collection of points
$\\{\omega_{s}\in\mathbb{R}^{|\mathcal{N}|},s\in\mathcal{S}\\}$ and respective
probability masses $\mathbb{P}(\omega=\omega_{s})=\frac{1}{|\mathcal{S}|}$,
$\forall s\in\mathcal{S}=\\{1,\dots,|\mathcal{S}|\\}$. Consider the rows
$a_{js}=a_{j}(\omega_{s})$, and the random $b_{js}=b_{j}(\omega_{s})$. Define
$p:=\left\lfloor\epsilon|\mathcal{S}|\right\rfloor$, the vector $y$ of binary
variables $y_{s}$, $\forall s\in\mathcal{S}$, and the large enough constants
$M_{js}$. Then, the MIP reformulation of problem (1) can be stated as:

$\displaystyle\min_{x}\quad$ $\displaystyle f(x)$ (2a) s.t. $\displaystyle
x\in X$ (2b) $\displaystyle a_{js}^{\top}x\leq b_{js}+M_{js}y_{s},\quad\forall
j\in\mathcal{J},s\in\mathcal{S}$ (2c)
$\displaystyle\sum_{s\in\mathcal{S}}y_{s}\leq p$ (2d) $\displaystyle
y_{s}\in\\{0,1\\},\quad\forall s\in\mathcal{S}.$ (2e)

Using the generic formulation (2), we present in Subsection 2.1 a procedure to
properly tune the values of the large constants $M_{js}$. We also explain in
this subsection how the intermediate results of the tightening procedure can
be efficiently used to remove constraints from set (2c) that are superfluous,
thus making model (2) more compact. Finally, we introduce in Subsection 2.2 a
set of valid inequalities that makes the linear relaxation of a particular
case of (2) remarkably tighter, and we relate our inequalities to the quantile
closure studied in Xie and Ahmed (2018).

###  2.1 Tightening and Screening

It is well-known that the linear relaxation of a Big-M formulation tends to
provide weak lower bounds in general (Conforti et al. (2014)). This is even
more so when the Big-Ms are chosen too loose. Constants $M_{js}$ in (2c)
should be set large enough for the corresponding constraints to be redundant
when the associated binary variables $y_{s}$ are equal to 1, and _as small as
possible_ to tighten the MIP formulation. To this end, Qiu et al. (2014)
provide an algorithm called “Iterative Coefficient Strengthening” that has
been successfully applied to other CCPs (see, e.g., Song et al. (2014)). A
customization of this procedure for the joint chance-constrained formulation
(2) is detailed in Algorithm 1.

Algorithm 1 Iterative Coefficient Strengthening

Input: The LHS and RHS vectors
$\\{a_{js}\\}_{j\in\mathcal{J},s\in\mathcal{S}}$ and coefficients
$\\{b_{js}\\}_{j\in\mathcal{J},s\in\mathcal{S}}$, respectively, and the
maximum allowed violation probability $p$ that determine the joint chance-
constraint system (2d), the deterministic feasible set $X$, and the total
number of iterations $\kappa$. Output: The large constants $M_{js}$, $\forall
j\in\mathcal{J}$, $s\in\mathcal{S}$.

  1. 1.

Initialization, $k=0$, $M_{js}^{0}=\infty$, $\forall j\in\mathcal{J}$,
$s\in\mathcal{S}$.

  2. 2.

For each $j\in\mathcal{J}$ and $s\in\mathcal{S}$ update $M_{js}^{k+1}$ as
follows: If $M_{js}^{k}<0$, then $M_{js}^{k+1}=M^{k}_{js}$. Otherwise,

$\displaystyle M^{k+1}_{js}=\max_{x,y_{s}}\quad$ $\displaystyle
a_{js}^{\top}x-b_{js}$ (3a) s.t. $\displaystyle x\in X$ (3b) $\displaystyle
a_{js}^{\top}x-b_{js}\leq M^{k}_{js}y_{s},\quad\forall
j\in\mathcal{J},s\in\mathcal{S}$ (3c)
$\displaystyle\sum_{s\in\mathcal{S}}y_{s}\leq p$ (3d) $\displaystyle 0\leq
y_{s}\leq 1,\quad\forall s\in\mathcal{S}.$ (3e)

  3. 3.

If $k+1<\kappa$, then $k=k+1$ and go to Step 2. Otherwise, stop.

The output of Algorithm 1 is the tuned Big-Ms that are input to the MIP
reformulation (2). This algorithm produces Big-Ms whose value either decreases
or remains equal at each iteration. There is, in fact, a number of iterations
beyond which the resulting Big-Ms converge. To reduce the computational burden
of running Algorithm 1, all problems (3) in Step 2 can be solved in parallel.
Hereinafter, we use the short name “T($\kappa$)” (from “Tightening”) to refer
to the solution of model (2) using the Big-M values given by the “Iterative
Coefficient Strengthening” algorithm with $\kappa$ iterations.

Equally important, Algorithm 1 can be easily upgraded to delete constraints
$(j,s)$ in (2c) that are redundant, and therefore, can be removed from problem
(2). Indeed, if Algorithm 1 delivers a large constant $M_{js}\leq 0$, then
constraint $j$ in scenario $s$ can be deleted from (2) without altering its
feasible region or its optimal solution. This is so because a non-positive
$M_{js}$ means that there is no $x$ satisfying (3b)–(3e) such that the
constraint takes on a value strictly greater than zero. Consequently, the
constraint is redundant in (2), since the feasibility region of (3) is a
relaxation of (2). This upgrade of method T not only makes formulation (2)
tighter through coefficient strengthening, but also more compact by screening
out redundant constraints. Naturally, the tightening and screening power of
algorithm T increases at each iteration. From now on, we use the short name
“TS($\kappa$)” (from “Tightening and Screening”) to refer to the strategy
whereby model (2) is solved without the constraints (2c) for which the value
of $M_{js}$ provided by Algorithm 1 after $\kappa$ iterations is lower than or
equal to 0.

While strengthening the parameters $M_{js}$ is a common strategy in the
technical literature to reduce the computational burden of CCPs, this is the
first time, to our knowledge, that intermediate results of Algorithm 1 are
used to eliminate superfluous constraints from model (2). We stress that the
screening process itself comes at no cost from Algorithm 1, while removing
superfluous constraints from (2) may substantially facilitate its solution. As
we show in Section 3.2, this is particularly true for the JCC-OPF.

Apart from making formulation (2) more compact, the screening of superfluous
constraints can also be used to accelerate the “Iterative Coefficient
Strengthening” process at each iteration. To do so, it suffices to modify
Algorithm 1 so that (3c) only includes the constraints for which
$M^{k}_{js}>0$. Thus, the number of constraints of model (3) is significantly
reduced and so is its solution time.

###  2.2 Valid inequalities

In this subsection, we propose valid inequalities that apply to JCCPs with a
particular structure: each row $a_{j}(\omega)$ in the technology matrix of (1)
can be rewritten as $a_{j}(\omega)=a^{0}_{j}+\Omega_{j}(\omega)\hat{a}_{j}$,
with $a_{j}^{0},\hat{a}_{j}\in\mathbb{R}^{|\mathcal{I}|}$ and where
$\Omega_{j}(\omega_{s})=\Omega_{js}$ is a real-valued function whose domain
includes the support of $\omega$. In other words, we derive a set of valid
inequalities to make the linear relaxation of problem

$\displaystyle\min_{x,y_{s}}\quad$ $\displaystyle f(x)$ (4a) s.t.
$\displaystyle x\in X$ (4b)
$\displaystyle\Omega_{js}\hat{a}_{j}^{\top}x-b_{js}+x^{\top}a_{j}^{0}\leq
M_{js}y_{s},\quad\forall j\in\mathcal{J},s\in\mathcal{S}$ (4c)
$\displaystyle\sum_{s\in\mathcal{S}}y_{s}\leq p$ (4d) $\displaystyle
y_{s}\in\\{0,1\\},\quad\forall s\in\mathcal{S}$ (4e)

tighter. Note that we do not make any assumptions on the sign of
$\hat{a}_{j}^{\top}$, $b_{js}$ and $a_{j}^{0}$. Furthermore, these valid
inequalities can also be added to the constraint set (3) of Algorithm 1, thus
dramatically increasing the tightening and screening power of TS. To
facilitate the comparative analysis carried out in Section 3.2, the so
upgraded algorithm is named “TS+V($\kappa$)” (from “Tightening and Screening
with Valid inequalities”).

To derive the set of valid inequalities for problem (4), we define the real
variables $z_{j}\in\mathbb{R}$ as $z_{j}:=\hat{a}_{j}^{\top}x$, and let
$z_{j}^{d}=\inf_{x\in X}\hat{a}_{j}^{\top}x$, $z_{j}^{u}=\sup_{x\in
X}\hat{a}_{j}^{\top}x$ denote the lower and upper bounds on $z$ induced by the
feasibility set (4b). Let us also define the function
$L_{js}:f_{js}(z_{j})=\Omega_{js}z_{j}-b_{js}$ for
$z_{j}\in[z_{j}^{d},z_{j}^{u}]$ and the set of functions
$\mathcal{L}_{j}:=\left\\{L_{js},\;\forall s\in\mathcal{S}\right\\}$. The
valid inequalities we propose are heavily supported by the concepts of _$k$
-lower_ and _$k$ -upper envelopes_, which we define in the following.

######  Definition 2.1.

For a given line $L_{js}$, we say that the point
$(\tilde{z},\tilde{t})\in\mathbb{R}^{2}$ _lies below_ , _on_ or _above_
function $L_{js}$ depending on whether
$\tilde{t}<\Omega_{js}\tilde{z}-b_{js}$,
$\tilde{t}=\Omega_{js}\tilde{z}-b_{js}$ or
$\tilde{t}>\Omega_{js}\tilde{z}-b_{js}$, respectively. Naturally, we also say
that function $L_{js}$ _lies above_ , _contains_ , or _lies below_ point
$(\tilde{z},\tilde{t})$ in these cases. We also say that a point
$(\tilde{z},\tilde{t})$ belongs to the set of lines $\mathcal{L}_{j}$ if there
exists a line $L_{js}\in\mathcal{L}_{j}$ that contains the point
$(\tilde{z},\tilde{t})$.

######  Definition 2.2.

For a set of lines $\mathcal{L}_{j}$, the _lower_ (resp.  _upper_) _score_ of
a point is the number of lines in $\mathcal{L}_{j}$ that lie below (resp.
above) that point. The _$k$ -lower_ (resp.  _$k$ -upper_) _envelope_ of a set
of lines $\mathcal{L}_{j}$ is the closure of the set of points that belong to
$\mathcal{L}_{j}$ and that have lower (resp. upper) score equal to $k-1$. The
$k$-lower envelope is also known as _$k$ -level_.

For the sake of illustration, Figure 1 shows in bold the $5$-upper envelope of
a set of 8 lines. Clearly, the $k$-envelopes of sets $\mathcal{L}_{j}$ can be
seen as piece-wise linear functions on $z_{j}$.

Figure 1: In bold, the 5-upper envelope (4-lower envelope, 4-level) of a set
$\boldsymbol{\mathcal{L}_{j}}$ of 8 lines in the plane. In dotted, the lower
hull of the 5-upper envelope

######  Proposition 2.3.

For a fixed $j\in\mathcal{J}$, let
$U_{j}^{p+1}(\mathchoice{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\displaystyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\textstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptstyle\bullet$}}}}}{\mathbin{\vbox{\hbox{\scalebox{0.5}{$\scriptscriptstyle\bullet$}}}}})$
be the (p+1)-upper envelope of the set of lines $\mathcal{L}_{j}$, with
$p:=\left\lfloor\epsilon|\mathcal{S}|\right\rfloor$. Then the inequality

$U_{j}^{p+1}(\hat{a}_{j}^{\top}x)+x^{\top}a_{j}^{0}\leq 0,\quad x\in X$ (5)

is valid for problem (4).

###### Proof.

Let $\bar{x},\bar{y}$ be any feasible solution of problem (4) with
$\bar{z}_{j}=\hat{a}_{j}^{\top}\bar{x}\in[z_{j}^{d},z_{j}^{u}]$, and suppose
that $U^{p+1}_{j}(\bar{z}_{j})+\bar{x}^{\top}a_{j}^{0}>0$. By definition of
$k$-upper envelope, there exist $p+1$ lines in $\mathcal{L}_{j}$,
$L_{js_{i}}:f_{js_{i}}(z_{j})=\Omega_{js_{i}}z_{j}-b_{js_{i}}$ $\forall
i\in\\{1,\dots,p+1\\}$ such that
$f_{js_{i}}(\bar{z}_{j})=\Omega_{js_{i}}\bar{z}_{j}-b_{js_{i}}\geq
U_{j}^{p+1}(\bar{z}_{j})$. Then for $i\in\\{1,\dots,p+1\\}$ it holds that
$\Omega_{js_{i}}\bar{z}_{j}-b_{js_{i}}+\bar{x}^{\top}a_{j}^{0}>0$.
Substituting in constraint (4c), we obtain that $\bar{y}_{s_{i}}M_{s_{i}}>0$
$\Rightarrow$ $\bar{y}_{s_{i}}=1$, for $i\in\\{1,\dots,p+1\\}$. But then
constraint (4d) is not satisfied, since
$\sum_{s\in\mathcal{S}}\bar{y}_{s}\geq\sum_{i=1}^{p+1}\bar{y}_{s_{i}}=p+1>p$.
This is, however, in contradiction with our initial statement that
$\bar{x},\bar{y}$ is a feasible solution of problem (4). ∎

The technical literature already includes references that propose
methodologies to determine the $k$-envelope of a set of linear functions. For
instance, the authors of Cheema et al. (2014) give a basic algorithm for
constructing $k$-envelopes called the Rider Algorithm. Essentially, this
algorithm is based on the fact that the $k$-envelope is an unbounded polygonal
chain that can be described by a sequence of vertices, which are intersections
of lines of the set. In fact, every point in the $k$-envelope is contained in
a given line. In this paper, we adapt the algorithm proposed in Cheema et al.
(2014) to the particular case in which $z_{j}^{d},z_{j}^{u}$ are finite.
Algorithm 2 describes our procedure in detail for a general set of the type
$\mathcal{L}_{j}$.

Algorithm 2 Rider Algorithm to build the $k$-upper envelope of
$\mathcal{L}_{j}$ (adapted to the bounded case)

To describe the $k$-upper envelope of a set of lines $\mathcal{L}_{j}$, we
derive the sequence of vertices of the polygonal chain (intersections of lines
in $\mathcal{L}_{j}$) $((z^{0},t^{0}),\dots,(z^{R},t^{R}))$ with
$(z_{j}^{d}=z^{0}<\dots<z^{R}=z_{j}^{u})$.

  1. 1.

Set $r=0$. Let $\\{(z^{d}_{j},\Omega_{js}z_{j}^{d}-b_{js}),\forall
s\in\mathcal{S}\\}$, be the intersections of the lines
$L_{js}\in\mathcal{L}_{j}$ with the vertical line $z=z_{j}^{d}$, and assume
w.l.o.g. that all the points have different upper scores. Let
$s_{0}\in\mathcal{S}$ be such that the intersection
$(z^{d}_{j},\Omega_{js_{0}}z_{j}^{d}-b_{js_{0}})$ has upper score equal to
$k-1$. Then $(z^{0},t^{0})=(z^{d}_{j},\Omega_{js_{0}}z_{j}^{d}-b_{js_{0}})$.

  2. 2.

Compute the value of $z_{j}$ for which the line $L_{js_{r}}$ intersects the
rest of the lines $L_{js}\in\mathcal{L}_{j}$ as
$\text{int}_{z}(s,s_{r})=\frac{b_{js}-b_{js_{r}}}{\Omega_{js}-\Omega_{js_{r}}}$.
If $\exists$ $s^{\prime}$ : $z^{r}<\text{int}_{z}(s^{\prime},s_{r})<z^{u}$, go
to Step 3. Otherwise, go to Step 4.

  3. 3.

Find the line that intersects $L_{js_{r}}$ at the leftmost point to the right
of $z^{r}$, that is, find
$s_{r+1}=\arg\min_{s}\left\\{\text{int}_{z}(s,s_{r}):\text{int}_{z}(s,s_{r})>z^{r}\right\\}$.
Set
$(z^{r+1},t^{r+1})=(\text{int}_{z}(s_{r+1},s_{r}),\Omega_{js_{r}}\text{int}_{z}(s_{r+1},s_{r})-b_{js_{r}})$,
update $r=r+1$ and go to Step 2.

  4. 4.

Set $(z^{R},t^{R})=(z^{u}_{j},\Omega_{js_{r}}z_{j}^{u}-b_{js_{r}})$.

The LHSs of the valid inequalities (5) are piece-wise linear functions not
necessarily convex. Therefore, the inclusion of these inequalities into model
(4) would require a significant amount of additional binary variables, which,
in turn, is expected to increase the computational burden of this problem even
further. Alternatively, we compute the lower hull of the $k$-upper envelope of
$\mathcal{L}_{j}$, which takes the form of a convex piece-wise linear
function. From this lower hull, we can extract a set of linear valid
inequalities (the linear extensions of the pieces) that can be seamlessly
inserted into model (4) without the need of any extra binary variables. In
doing so, we are able to tighten model (4), which can thus be solved more
efficiently by available optimization software. Before presenting the set of
valid inequalities, the following definition is required.

######  Definition 2.4.

Let $Z$ be the convex hull of a set of points $P$. The _upper_ (resp.
_lower_) _hull_ of $P$ is the set of edges of $Z$ that lie on or above (resp.
on or below) every point in $P$.

Corollary 2.5 presents the set of linear valid inequalities (6) given by the
lower hull of the $k$-upper envelope of $\mathcal{L}_{j}$.

######  Corollary 2.5.

Let $\\{(z^{r},t^{r})\\}$, $r\in\\{0,\dots,R\\}$, be the ordered set of
vertices obtained by applying Algorithm 2 to set $\mathcal{L}_{j}$, and let
$\\{(z^{r^{\prime}},t^{r^{\prime}})\\}$,
$r^{\prime}\in\\{0,\dots,R^{\prime}\\}\subseteq\\{0,\dots,R\\}$, be the
ordered subset of vertices such that the associated polygonal chain is the
lower hull. Then the following linear inequalities

$\frac{t^{r^{\prime}+1}-t^{r^{\prime}}}{z^{r^{\prime}+1}-z^{r^{\prime}}}(\hat{a}_{j}^{\top}x-z^{r^{\prime}})+t^{r^{\prime}}\leq-x^{\top}a_{j}^{0},\quad
x\in X,r^{\prime}\in\\{0,\dots,R^{\prime}-1\\}$ (6)

are valid for problem (4).

###### Proof.

The proof is straightforward, since for each $x\in X$ it holds
$\frac{t^{r^{\prime}+1}-t^{r^{\prime}}}{z^{r^{\prime}+1}-z^{r^{\prime}}}(\hat{a}_{j}^{\top}x-z^{r^{\prime}})+t^{r^{\prime}}\leq
U_{j}^{p+1}(\hat{a}_{j}^{\top}x)\leq-x^{\top}a_{j}^{0}$, by hypothesis and
using (5). ∎

There exist plenty of algorithms of convexification of a set of points in the
plane. Two of the most well-known are the _Jarvis march_ and the _Graham scan_
(Tóth et al. (2017)). Here, we give a simplified version of the former where
we make use of special features of our set of points and only compute the
lower hull. In particular, we assume that we have a set of presorted points
$\\{(z^{r},t^{r})\\}$, $r\in\\{0,\dots,R\\}$ whose first and last points
always belong to the hull. To speed up the process, we can find the point
$(z^{\bar{r}},t^{\bar{r}})$ from the previous set with the lowest
$t$-coordinate (which always belongs to the lower hull), and then apply the
algorithm to the subsets $\\{(z^{0},t^{0}),\dots,(z^{\bar{r}},t^{\bar{r}})\\}$
and $\\{(z^{\bar{r}},t^{\bar{r}}),\dots,(z^{R},t^{R})\\}$. Algorithm 3
describes in detail the proposed convexification procedure.

Algorithm 3 Jarvis March Algorithm to obtain the lower hull of a set of
presorted points $P$

Let $P=\\{(z^{r},t^{r})\\}$, $r\in\\{0,\dots,R\\}$ be a set of points with
$z^{0}<\dots<z^{R}$. The algorithm derives a subset
$P^{\prime}=\\{(z^{r^{\prime}},t^{r^{\prime}})\\}$,
$r^{\prime}\in\\{0,\dots,R^{\prime}\\}\subseteq\\{0,\dots,R\\}$, which
constitutes the lower hull of the first set.

  1. 1.

Initially, set $P^{\prime}=\\{(z^{0},t^{0})\\}$.

  2. 2.

Assume $(z^{r^{\prime}},t^{r^{\prime}})$ is the last point included in
$P^{\prime}$. If $r^{\prime}=R$, the lower hull of $P$ is given by the set of
points $P^{\prime}$. Otherwise, let
$r^{\prime}+1:=\arg\min_{r}\left\\{\frac{t^{r}-t^{r^{\prime}}}{z^{r}-z^{r^{\prime}}}:(z^{r},t^{r})\in
P\text{ with }z^{r}>z^{r^{\prime}}\right\\}$. Update
$P^{\prime}=P^{\prime}\cup\\{(z^{r^{\prime}+1},t^{r^{\prime}+1})\\}$. Repeat
Step 2.

Algorithm 2 constructs the $k$-upper envelope of $\mathcal{L}$ in
$O(n_{k}\log^{2}|\mathcal{L}|)$ time and $O(n_{k}|\mathcal{L}|+|\mathcal{L}|)$
space, where $n_{k}$ denotes the maximum number of edges of the $k$-upper
envelope of any set with $|\mathcal{L}|$ lines (see Edelsbrunner and Welzl
(1986)). An upper bound on $n_{k}=O(|\mathcal{L}|k^{\frac{1}{2}})$ can be
found in Edelsbrunner and Welzl (1985). Dey (1998), Tóth (2000), Nivasch
(2008) have subsequently improved the upper bound of the complexity of this
problem by studying a closely related problem, the $k$-set problem, in the
two-dimensional case. As for the Jarvis March Algorithm, in the two-
dimensional case presented here it has running time $O((R+2)(R^{\prime}+1))$
(Tóth et al. (2017)). For detailed results on the complexity of related
algorithms and for extensions to the multidimensional case, we refer the
reader to Appendix B.

Valid inequalities (5) and Algorithms 2 and 3 can be generalized to work with
finite discrete distributions with unequal probabilities. Observe that the
definition of the $(p+1)$-lower (resp. $(p+1)$-upper) envelope is based on
satisfying constraint (4d). In this extension, constraint (4d) is replaced
with $\sum_{s\in\mathcal{S}}p_{s}y_{s}\leq\epsilon$, where $p_{s}$ is the
probability of scenario $s$. Hence, it suffices to redefine the lower (resp.
upper) score of a point as the probability of the lines that lie below (resp.
above) that point, and modify the definition of $k$-lower (resp. upper)
envelope and Step 3 of Algorithm 2 accordingly, so that this more general
constraint is satisfied.

The following theorem sheds light on the strength of our valid inequalities.

######  Theorem 2.6.

Let $\mathcal{FR}_{j}$ be the feasible region of the relaxation of problem (4)
associated to a given row $j\in\mathcal{J}$, i.e.,

$\displaystyle\min_{x,y_{s}}\quad$ $\displaystyle f(x)$ (7a) s.t.
$\displaystyle\eqref{MIPGenCC_xinX},\eqref{MIPGenCC_violation},\eqref{MIPGenCC_ybinaria}$
(7b) $\displaystyle\Omega_{js}\hat{a}_{j}^{\top}x-b_{js}+x^{\top}a_{j}^{0}\leq
M_{js}y_{s},\quad\forall s\in\mathcal{S}.$ (7c)

Assuming that the set $X$ of deterministic constraints given in (4b) is non-
empty and compact, the set (6) of valid inequalities for all $j\in\mathcal{J}$
provides $\bigcap_{j\in\mathcal{J}}\mathrm{conv}(\mathcal{FR}_{j})$.

The proof of this theorem is a direct consequence of Proposition A.3 in
Appendix A.

Furthermore, it is possible to relate our valid inequalities with the
_quantile cuts_ described in Qiu et al. (2014) and Xie and Ahmed (2018). To
see this, consider problem (4) with feasible region $\mathcal{FR}$ and note
that $\mathcal{FR}=\bigcap_{j\in\mathcal{J}}\mathcal{FR}_{j}$. The quantile
cuts can be viewed as a projection of the well-known family of mixing
inequalities in the $(x,y)$-space onto the original space $x$. Quantile cuts
represent an infinite family of inequalities, and the intersection of all
quantile cuts is called the quantile closure. If another round of quantile
cuts is applied to a stronger formulation, specifically the original
formulation with the (first) quantile closure, a subsequent closure can be
derived (the second closure), and so on. In Xie and Ahmed (2018) it is proved
that the sequence of sets obtained by successive quantile closure operations
converges to the convex hull of the feasible region of the original problem
with respect to the Hausdorff metric, $\mathrm{conv}(\mathcal{FR})$. It is
also shown that the separation over the first quantile closure is NP-hard.
Corollary A.12 in Appendix A shows that the set (6) of valid inequalities (per
constraint $j\in\mathcal{J}$) coincides with the limit of the succession of
quantile closures of $\mathcal{FR}_{j}$ as the order of the closure grows to
infinity. This limit is precisely $\mathrm{conv}(\mathcal{FR}_{j})$.

##  3 Optimal Power Flow under Uncertainty: A Joint Chance-Constrained
Modeling Approach

The OPF is a routine at the core of important tools for power system
operations (Frank et al. (2012)). In its deterministic version, the OPF
problem seeks to determine the least-costly dispatch of thermal generating
units to satisfy the system’s _net_ demand (i.e., demand minus renewable
generation), while complying with the technical limits of production and
transmission network equipment. The main challenge of the OPF is that it is a
non-linear and non-convex optimization problem, due to the power flow
equations that govern the (static) behavior of power systems. For this reason,
the _direct current_ approximation (DC) of the power flow equations, which
transforms the problem into a linear program, is frequently used. The demand
and renewable generation are factors that increase the uncertainty in power
systems, and ignoring it can lead to unsafe operating conditions.

However, given the inherently uncertain nature of the electricity net demand,
the probabilistic version of the OPF problem can be formulated as a JCCP that
aims at minimizing the expected production cost while enforcing that the
technical constraints are satisfied with a given (high) probability. In this
context, chance-constrained programming can be used to minimize the expected
operating cost whilst guaranteeing that the system withstands unforeseen peeks
of electrical load due to stochastic demand or uncertainty in power generation
(Van Ackooij et al. (2011)). The chance-constrained OPF (CC-OPF) problem
addresses this uncertainty and pursues to ensure the safe operation of a power
system with a high level of probability.

To address the CC-OPF problem, several papers in the literature (e.g., Roald
and Andersson (2018)) directly work with SCCs. However, the main drawback of
this modeling approach is that, even in those cases where the probability of
violating each individual constraint seems more than tolerable, the resulting
_joint risk_ (that is, the probability that _any_ of the technical constraints
be violated) may still be excessive and inadmissible. This is the key
motivation behind the use of JCCs to tackle the CC-OPF problem (see, e.g.,
Peña-Ordieres et al. (2021)).

In this vein, there are several approaches in the literature to solve the JCC-
OPF problem. Vrakopoulou et al. (2013) adopt the scenario approach (SA) to
approximate the solution of the JCC-OPF, while Chen et al. (2021) propose a
heuristic data-driven method that involves enforcing the satisfaction of the
technical constraints for a box of the uncertainty. This box is inferred using
one-class support vector clustering and its size is contingent on the system’s
desired reliability. Hou and Roald (2020) propose an iterative tuning
algorithm to solve a robust reformulation of the JCC-OPF problem. Esteban-
Pérez and Morales (2023) introduce a distributionally robust JCC-OPF model
that considers contextual information using an ambiguity set based on
probability trimmings. To make their model tractable, they resort to the
widely known CVaR-based approximation of the JCC.

The aforementioned SAA method is another effective way to solve JCCPs and has
the potential to identify OPF solutions with a better cost performance than
that of the more conservative solutions delivered by the previous approaches.
However, solving the JCC-OPF problem using SAA is challenging due to the
presence of binary variables, the number of scenarios required and the size of
the power systems. Lejeune and Dehghanian (2020) propose a methodology to
solve the SAA of the JCC-OPF without including the power flow equations into
the joint chance-constrained system. To the best of our knowledge, we are the
first to efficiently solve the JCC-OPF problem by means of the SAA approach,
using a MIP reformulation and including the arduous power flow constraints.
Furthermore, unlike the sample-based approach introduced in Peña-Ordieres et
al. (2021), which relies on a smooth nonlinear approximation of the JCC-OPF,
ours offers optimality guarantees.

###  3.1 Notation and problem formulation

In this subsection, we present the formulation of the JCC-OPF problem that we
consider in the case study. To do so, we introduce the following notation and
modeling choices, widely used in the power system domain:

  1. 1.

_Power system:_ A power system can be represented by a directed graph denoted
by $G(\mathcal{N},\mathcal{L})$, where $\mathcal{N}$ is the set of buses
(vertices), indexed by $n$, and $\mathcal{L}$ is the set of transmissions
lines (edges), whose direction is chosen arbitrarily and are indexed by $l$.
To supply the net demand at each bus, the power system is equipped with a set
of generators $\mathcal{G}$, indexed by $g$, where $\mathcal{G}_{n}$
represents the set of generators at bus $n$. Note that the energy injected at
each bus, derived from the balance between generation and net demand, is
transported among the buses through the transmission lines. For illustration,
Figure 2 depicts a representation of a power system with three buses, three
transmission lines, two generators located at nodes 1 and 2 and one
electricity net demand at node 3.

$n=1$$n=2$$n=3$$g=1$$g=2$$l=1$ $l=2$$l=3$$d_{3}$ Figure 2: Illustrative three-
node power system

  2. 2.

_Nodal net loads_ : The (uncertain) electricity net demand at node $n$,
$\tilde{d}_{n}$, is given by $\tilde{d}_{n}=d_{n}-\omega_{n}$, where $d_{n}$
is the predicted value and $\omega_{n}$ is the forecast error with a change of
sign. This error is modeled as a random variable with zero mean which follows
an unknown continuous probability distribution.

  3. 3.

_Generation_ : To cope with the forecast errors
$(\omega_{n})_{n\in\mathcal{N}}$, generators’ power outputs are adjusted
according to the following affine control policy:

$\displaystyle\tilde{p}_{g}=p_{g}-\beta_{g}\Omega,\quad\forall
g\in\mathcal{G},$

where $\Omega:=\sum_{n\in\mathcal{N}}\omega_{n}$ is the system-wise aggregated
forecast error, and $p_{g}$ and $\beta_{g}$ are the power output dispatch and
the participation factor of generating unit $g$, respectively (see, e.g.,
Bienstock et al. (2014), Hou and Roald (2020), Peña-Ordieres et al. (2021)).
The minimum and maximum capacity of generator $g$ is denoted by
$\underline{p}_{g}$ and $\overline{p}_{g}$, respectively.

  4. 4.

_Power balance_ : Given the affine control policy of the previous point, the
power balance equation takes the following form:

$\displaystyle\sum_{g\in\mathcal{G}}\tilde{p}_{g}-\sum_{n\in\mathcal{N}}\tilde{d}_{n}=\sum_{g\in\mathcal{G}}\left(p_{g}-\beta_{g}\Omega\right)-\sum_{n\in\mathcal{N}}\left(d_{n}-\omega_{n}\right)=0.$

Hence, to ensure the power balance for _any_ realization of the forecast
errors $(\omega_{n})_{n\in\mathcal{N}}$, it must hold:

$\displaystyle\sum_{g\in\mathcal{G}}p_{g}-\sum_{n\in\mathcal{N}}d_{n}$
$\displaystyle=$ $\displaystyle 0$
$\displaystyle\sum_{g\in\mathcal{G}}\beta_{g}$ $\displaystyle=$ $\displaystyle
1.$

  5. 5.

_Power flows_ : Line flows are modeled using the well-known approximation
based on the power transfer distribution factors (PTDFs), $B_{ln}$,
$l\in\mathcal{L}$, $n\in\mathcal{N}$, which sets a linear relation between the
power flow through line $l$ and the power injected at node $n$. The maximum
capacity of line $l$ is denoted by $\overline{f}_{l}$.

  6. 6.

_Power production cost_ : The cost function of each generating unit is assumed
to be quadratic and, as a result, the total power production cost is given by

$\displaystyle\sum_{g\in\mathcal{G}}C_{2,g}\,\left(p_{g}-\Omega\beta_{g}\right)^{2}+C_{1,g}\,\left(p_{g}-\Omega\beta_{g}\right)+C_{0,g},$

where $C_{2,g}$, $C_{1,g}$, $C_{0,g}$ are the coefficients defining the
quadratic cost function of generating unit $g$. On the assumption that
$\omega_{n}$, for each $n\in\mathcal{N}$, is a random variable with zero mean,
we have (see, for instance, Peña-Ordieres et al. (2021))

$\displaystyle\mathbb{E}\left[\sum_{g\in\mathcal{G}}C_{2,g}\,\left(p_{g}-\Omega\beta_{g}\right)^{2}+C_{1,g}\,\left(p_{g}-\Omega\beta_{g}\right)+C_{0,g}\right]=\sum_{g\in\mathcal{G}}C_{2,g}\,p_{g}^{2}+C_{1,g}\,p_{g}+C_{0,g}+\mathbb{V}\left(\Omega\right)\,C_{2,g}\,\beta_{g}^{2},$

where $\mathbb{V}(\Omega)$ denotes the variance of the random variable
$\Omega$.

With the above ingredients, the JCC-OPF problem that we tackle in this paper
is formulated as follows:

$\displaystyle\min_{p_{g},\beta_{g}\geq 0,\forall g\in\mathcal{G}}\quad$
$\displaystyle\sum_{g\in\mathcal{G}}C_{2,g}\,p_{g}^{2}+C_{1,g}\,p_{g}+C_{0,g}+\mathbb{V}\left(\Omega\right)\,C_{2,g}\,\beta_{g}^{2}$
(8a) s.t. $\displaystyle\sum_{g\in\mathcal{G}}\beta_{g}=1$ (8b)
$\displaystyle\sum_{g\in\mathcal{G}}p_{g}-\sum_{n\in\mathcal{N}}d_{n}=0$ (8c)
$\displaystyle\underline{p}_{g}\leq p_{g}\leq\overline{p}_{g},\quad\forall
g\in\mathcal{G}$ (8d)
$\displaystyle-\overline{f}_{l}\leq\sum_{n\in\mathcal{N}}B_{ln}\left(\sum_{g\in\mathcal{G}_{n}}p_{g}-d_{n}\right)\leq\overline{f}_{l},\quad\forall
l\in\mathcal{L}$ (8e)
$\displaystyle\mathbb{P}\left(\begin{array}[]{l}\underline{p}_{g}\leq
p_{g}-\Omega\beta_{g}\leq\overline{p}_{g},\quad\forall g\in\mathcal{G}\\\
-\overline{f}_{l}\leq\displaystyle\sum_{n\in\mathcal{N}}B_{ln}\left(\displaystyle\sum_{g\in\mathcal{G}_{n}}\left(p_{g}-\Omega\beta_{g}\right)+\omega_{n}-d_{n}\right)\leq\overline{f}_{l},\quad\forall
l\in\mathcal{L}\end{array}\right)\geq 1-\epsilon.$ (8h)

The objective (8a) is the minimization of the expected total generation cost.
The equality constraints (8b) and (8c) enforce the power balance in the
system, whereas constraints (8d) and (8e) ensure a feasible power dispatch
which corresponds to an error-free scenario, i.e., to a realization of the
net-load forecast errors such that $\omega_{n}=0$ $\forall n\in\mathcal{N}$.
Finally, expression (8h) constitutes the joint chance-constraint system by
which the decision-maker states that the OPF solution must be feasible with a
probability greater than or equal to $1-\epsilon$. Accordingly, parameter
$\epsilon$ is the maximum allowed probability of constraint violation set by
the user. Formulation (8) is quite standard and has been used before by
Bienstock et al. (2014), Hou and Roald (2020), Peña-Ordieres et al. (2021),
among others.

Problem (8) can be written in the form of (1), i.e., as a CCP with linear JCC
and random RHS and LHS. To see this, define the vector of continuous decision
variables $x$ in (1) as $x:=(p_{g},\beta_{g})_{g\in\mathcal{G}}$, and group
all these variables by means of the set $\mathcal{I}$ with elements $i$
running from 1 to $|\mathcal{I}|=2|\mathcal{G}|$. In this way, we have that
$x\in\mathbb{R}^{|\mathcal{I}|}_{+}$, the set
$X\subseteq\mathbb{R}^{|\mathcal{I}|}$ represents the polyhedron defined by
the deterministic constraints (8b)–(8e), and
$f:\mathbb{R}^{|\mathcal{I}|}\longrightarrow\mathbb{R}$ is the convex function
providing the expected total generation cost (8a). Likewise, if we collect all
the constraints involved in the joint chance-constraint system (8h) into the
set $\mathcal{J}$ (hence $|\mathcal{J}|=2|\mathcal{G}|+2|\mathcal{L}|$), this
system can be represented by way of a technology matrix with random rows
$a_{j}(\omega)\in\mathbb{R}^{|\mathcal{I}|}$, $j\in\mathcal{J}$, and random
RHS $b_{j}(\omega)\in\mathbb{R}$, $j\in\mathcal{J}$, where the uncertainty is
again represented through the random vector $\omega$ taking values in
$\mathbb{R}^{|\mathcal{N}|}$.

As discussed in Section 2, the CCP (8) can be easily reformulated into a MIP
using SAA. Thus, we assume that $\omega$ has a finite discrete support defined
by a collection of points
$\\{\omega_{s}\in\mathbb{R}^{|\mathcal{N}|},s\in\mathcal{S}\\}$ and respective
probability masses $\mathbb{P}(\omega=\omega_{s})=\frac{1}{|\mathcal{S}|}$,
$\forall s\in\mathcal{S}=\\{1,\dots,|\mathcal{S}|\\}$. Accordingly,
$\omega_{ns}$ and $\Omega_{s}$ are realizations of the respective random
variables under scenario $s$. We define
$p:=\left\lfloor\epsilon|\mathcal{S}|\right\rfloor$, the vector $y$ of binary
variables $y_{s}$, $\forall s\in\mathcal{S}$, and the large enough constants
$M^{1}_{gs},M^{2}_{gs},M^{3}_{ls},M^{4}_{ls}$. Thus, the MIP reformulation of
problem (8) writes as follows:

$\displaystyle\min_{p_{g},\beta_{g}\geq 0,y_{s}}\quad$
$\displaystyle\sum_{g\in\mathcal{G}}C_{2,g}\,p_{g}^{2}+C_{1,g}\,p_{g}+C_{0,g}+\widehat{\mathbb{V}}\left(\Omega\right)\,C_{2,g}\,\beta_{g}^{2}$
(9a) s.t. $\displaystyle\eqref{eq:OPF_balance1}-\eqref{eq:OPF_flow-det}$ (9b)
$\displaystyle
p_{g}-\Omega_{s}\beta_{g}\geq\underline{p}_{g}-y_{s}M^{1}_{gs},\quad\forall
g\in\mathcal{G},s\in\mathcal{S}$ (9c) $\displaystyle
p_{g}-\Omega_{s}\beta_{g}\leq\overline{p}_{g}+y_{s}M^{2}_{gs},\quad\forall
g\in\mathcal{G},s\in\mathcal{S}$ (9d)
$\displaystyle\sum_{n\in\mathcal{N}}B_{ln}\left(\sum_{g\in\mathcal{G}_{n}}\left(p_{g}-\Omega_{s}\beta_{g}\right)-d_{n}+\omega_{ns}\right)\geq-\overline{f}_{l}-y_{s}M^{3}_{ls},\quad\forall
l\in\mathcal{L},s\in\mathcal{S}$ (9e)
$\displaystyle\sum_{n\in\mathcal{N}}B_{ln}\left(\sum_{g\in\mathcal{G}_{n}}\left(p_{g}-\Omega_{s}\beta_{g}\right)-d_{n}+\omega_{ns}\right)\leq\overline{f}_{l}+y_{s}M^{4}_{ls},\quad\forall
l\in\mathcal{L},s\in\mathcal{S}$ (9f)
$\displaystyle\sum_{s\in\mathcal{S}}y_{s}\leq p$ (9g) $\displaystyle
y_{s}\in\\{0,1\\},\quad\forall s\in\mathcal{S}.$ (9h)

Constraints (9c)-(9f) represent the sample-based reformulation of the joint
chance constraint (8h). For a given scenario $s\in\mathcal{S}$, inequalities
(9c)-(9f) guarantee that all the original constraints are satisfied when
$y_{s}=0$. If $y_{s}=1$, some of the original constraints can be violated for
scenario $s$. Finally, the inequality (9g) ensures that the probability of the
JCC is met and the binary character of variables $y_{s}$ is declared in (9h).

###  3.2 Numerical Experiments

This section discusses a series of numerical experiments with which we
evaluate the different approaches presented in Section 2 to solve the SAA-
based MIP reformulation of the JCC-OPF. In particular, we compare the
performance of approaches T, TS and TS+V using five standard power systems
widely employed in the technical literature on the topic, namely the IEEE-
RTS-24, IEEE-57, IEEE-RTS-73, IEEE-118, and IEEE-300 test systems. All data
pertaining to these systems are publicly available in the repository Power
Grid Lib (2022) under version 21 and their main features are listed in Table
1. All optimization problems have been solved using GUROBI 9.1.2 (Gurobi
Optimization, LLC (2022)) on a Linux-based server with CPUs clocking at 2.6
GHz, 6 threads and 32 GB of RAM. In all cases, the optimality GAP has been set
to $10^{-9}\%$ and the time limit to 10 hours. All data used and codes
implemented throughout this numerical experiments can be found in the
repository OASYS (2023).

Table 1: Short Description of Test Power Systems | IEEE-RTS-24 | IEEE-57 | IEEE-RTS-73 | IEEE-118 | IEEE-300
---|---|---|---|---|---
# Nodes | 24 | 57 | 73 | 118 | 300
# Generators | 32 | 4 | 96 | 19 | 57
# Lines | 38 | 41 | 120 | 186 | 411

Similarly to Peña-Ordieres et al. (2021), we assume that the error of net
loads is normally distributed, i.e., $\omega\sim N(\mathbf{0},\Sigma)$, where
$\mathbf{0}$ and $\Sigma$ represent, respectively, the zero vector and the
covariance matrix. We also assume that the standard deviation of $\omega_{n}$
at node $n$ is proportional to the net nodal demand $d_{n}$ according to a
parameter $\zeta$ between 0 and 1. Thus, this parameter controls the magnitude
of net demand fluctuations. Under these assumptions, the procedure to model
uncertainty proposed in Peña-Ordieres et al. (2021) runs as follows. First, we
compute the positive definite matrix $C=\widehat{C}\widehat{C}^{\top}$ where
each element of matrix $\widehat{C}$ is a sample randomly drawn from a uniform
distribution with support in $[-1,1]$. Then, to obtain a positive definite
matrix $\Sigma$ in which the diagonal elements are equal to $(\zeta
d_{n})^{2}$, we define each of its entries ($\sigma_{nn^{\prime}}$) as
follows:

$\displaystyle\sigma_{nn^{\prime}}=\zeta^{2}\frac{c_{nn^{\prime}}}{\sqrt{c_{nn}c_{n^{\prime}n^{\prime}}}}d_{n}d_{n^{\prime}},\quad\forall
n,n^{\prime}\in\mathcal{N},$

where $c_{nn^{\prime}}$ denotes the element of matrix $C$ in row $n$ and
column $n^{\prime}$. To avoid generating infeasible instances of the JCC-OPF
problem, the parameter $\zeta$ has been set to 0.15 for the four smallest
systems and to 0.05 for the IEEE-300 system. To characterize the net demand
uncertainty, we consider 1000 scenarios and a tolerable probability of
violation of the joint chance constraint of 5% (i.e., $\epsilon=0.05$ and
$p=50$). Finally, each solution strategy is run for ten different sets of
randomly generated samples. Accordingly, in this section we provide tables
with figures averaged over these ten instances. Also, the generated samples
can be downloaded from the repository OASYS (2023).

Table 2 includes the results of solving the mixed-integer quadratic
optimization model (9) if the large constants $M^{1}_{gs}$, $M^{2}_{gs}$,
$M^{3}_{ls}$ and $M^{4}_{ls}$ are set to a high enough arbitrary value,
specifically $10^{4}$. Despite being remarkably computationally expensive,
this approach has been used in the technical literature (e.g., Zhang et al.
(2015)), and thus we refer to it as _benchmark approach_ (BN). Table 2
includes the number of constraints in the model (#CON), the linear relaxation
gap (LR-GAP) calculated using the optimal solution of each instance, the
optimality gap given by the difference between the best lower bound and the
best integer solution found by the MIP solver (MIP-GAP), the number of
instances solved to global optimality in less than 10 hours (#OPT) and the
solution time in seconds (Time). As expected, the computational time needed to
solve the OPF with the Big-M model (9) increases significantly with the size
of the instances. While the 10 instances from systems IEEE-RTS-24, IEEE-57 and
IEEE-73 are solved in less than 10 hours, none of the instances for systems
IEEE-118 and IEEE-300 are solved to global optimality within that time limit,
and the average MIP-GAP after 10 hours amounts to 0.29% and 0.27%,
respectively. Interestingly, despite the fact that the LR-GAP is relatively
low for the IEEE-RTS-73 system, the average computational time required to
solve this case is particularly high compared to the two smaller systems.

Table 2: Benchmark approach (BN): Results | IEEE-RTS-24 | IEEE-57 | IEEE-RTS-73 | IEEE-118 | IEEE-300
---|---|---|---|---|---
#CON | 140143 | 168171 | 432435 | 410413 | 936939
LR-GAP | 1.756% | 0.623% | 0.061% | 0.956% | 1.114%
MIP-GAP (#OPT) | 0.00% (10) | 0.00% (10) | 0.00% (10) | 0.29% (0) | 0.27% (0)
Time (s) | 1121.3 | 103.2 | 11161.2 | 36000.0 | 36000.0

As discussed in the technical literature, a proper tuning of the Big-Ms makes
model (9) tighter and generally easier to solve (Qiu et al. (2014)) by the MIP
routine. Therefore, we evaluate the computational performance of the
“Iterative Coefficient Strengthening” Algorithm and provide the corresponding
results in Table 3. In particular, T(1), T(2) and T(3) represent the results
obtained by solving model (9) with the Big-M values provided by Algorithm 1
with $\kappa=1$, 2 and 3, respectively. Table 3 includes the average values of
LR-GAP and MIP-GAP, the number of instances solved to optimality in less than
10 hours (#OPT) and the speedup factor with respect to the benchmark approach.
To determine this factor, we have considered that the total computational time
of approach T is given as the sum of the time required to run Algorithm 1
$\kappa$ times to determine the Big-Ms plus the time it takes to solve problem
(9).

Table 3: Coefficient tightening approach (T): Results | IEEE-RTS-24 | IEEE-57 | IEEE-RTS-73 | IEEE-118 | IEEE-300
---|---|---|---|---|---
LR-GAP |  T(1) | 1.755% | 0.510% | 0.061% | 0.711% | 0.472%
T(2) | 1.662% | 0.330% | 0.055% | 0.522% | 0.324%
T(3) | 1.386% | 0.255% | 0.029% | 0.434% | 0.264%
MIP-GAP (#OPT) |  T(1) | 0.00% (10) | 0.00% (10) | 0.00% (10) | 0.27% (0) | 0.16% (0)
T(2) | 0.00% (10) | 0.00% (10) | 0.03% (2) | 0.16% (0) | 0.09% (0)
T(3) | 0.00% (10) | 0.00% (10) | 0.00% (10) | 0.12% (0) | 0.07% (0)
Speedup factor |  T(1) | 0.22x | 0.07x | 0.61x | 1.00x | 1.00x
T(2) | 0.17x | 0.13x | 0.31x | 1.00x | 1.00x
T(3) | 0.48x | 0.24x | 0.66x | 1.00x | 1.00x

Since reducing the Big-M values makes model (9) tighter, the results in Table
3 show lower values of LR-GAP with respect to BN. Furthermore, this effect
grows with the number of iterations since Algorithm 1 ensures that the Big-Ms
never increase between iterations. Although decreasing the values of the Big-
Ms leads to tighter MIPs for all the test systems, the numerical results in
Table 3 clearly indicate that computational savings are not guaranteed in all
cases. Indeed, while the ten instances are solved by BN in less than 10 hours
for the IEEE-RTS-73 system, T(2) only provides the optimal solution for two
instances. On top of that, the speedup factors for the three smaller systems
are always lower than 1, which means that the computational times actually
increase in these cases. On the contrary, the average MIP-GAP of the two
largest systems is significantly decreased with respect to BN. Therefore, we
conclude that, due to the heuristics implemented in current commercial MIP
solvers, the computational advantages that one could expect a priori from
“Iterative Coefficient Strengthening” are not always guaranteed and are
contingent on the structure and data of the problem.

Next, in Table 4 we provide the computational results related to the TS
method, in which Algorithm 1 is extended to remove redundant constraints from
model (9). Here, #CON is provided as the percentage of the number of
constraints of the reference model BN (indicated in Table 2) that are retained
by TS in each iteration. Table 4 also includes the average MIP-GAP, the number
of instances solved to optimality and the average speedup factor in relation
to BN.

Table 4: Tightening and Screening (TS): Results | IEEE-RTS-24 | IEEE-57 | IEEE-RTS-73 | IEEE-118 | IEEE-300
---|---|---|---|---|---
#CON |  TS(1) | 23.9% | 2.8% | 26.0% | 8.9% | 12.7%
TS(2) | 23.3% | 2.2% | 23.8% | 6.3% | 9.1%
TS(3) | 23.1% | 2.1% | 23.0% | 5.6% | 8.0%
MIP-GAP (#OPT) |  TS(1) | 0.00% (10) | 0.00% (10) | 0.00% (10) | 0.15% (0) | 0.08% (0)
TS(2) | 0.00% (10) | 0.00% (10) | 0.00% (10) | 0.03% (2) | 0.04% (0)
TS(3) | 0.00% (10) | 0.00% (10) | 0.00% (10) | 0.01% (6) | 0.01% (4)
Speedup factor |  TS(1) | 1.5x | 1.8x | 4.7x | 1.0x | 1.0x
TS(2) | 1.5x | 3.8x | 2.6x | 1.1x | 1.0x
TS(3) | 3.8x | 4.4x | 15.1x | 1.4x | 1.2x

Table 4 shows that the upgraded Algorithm 1 screens out a large percentage of
the constraints in model (9), only retaining between 2% and 26% of them. The
results in this table demonstrate that combining the tightening of the Big-Ms
and the elimination of superfluous constraints leads to significant
computational savings. For instance, the speedup factor for the three smallest
systems ranges now between 1.5x and 15.1x. In addition, TS(3) is able to solve
six and four instances to optimality for systems IEEE-118 and IEEE-300,
respectively, in less than 10 hours, and the average MIP-GAP is reduced to
0.01% in these two largest power systems. Therefore, the computational
performance of model (9) has been drastically improved by combining the
strengthening of the Big-Ms (making model (9) tighter) and the removal of
redundant constraints (making model (9) more compact). Equally important, the
elimination of superfluous constraints notably reduces the MIP solver’s need
for RAM memory, which decreases from around 100 GB in T to 32 GB in TS.

We continue the numerical experiments by evaluating the impact of including
the valid inequalities derived in Section 2.2 in model (9). The so obtained
results are collated in Table 5. In what follows, this approach is called BN+V
for short. The results in Table 5 comprise, in order and following the
previous notation, the average number of constraints, the average values of
LR-GAP and MIP-GAP, the number of instances solved to optimality, and the
average speedup factor with respect to the BN approach of Table 2.

Table 5: Tightening by valid inequalities (BN+V): Results | IEEE-RTS-24 | IEEE-57 | IEEE-RTS-73 | IEEE-118 | IEEE-300
---|---|---|---|---|---
#CON | 101.0% | 101.2% | 101.2% | 101.6% | 101.5%
LR-GAP | 0.3374% | 0.2038% | 0.0001% | 0.4784% | 0.3192%
MIP-GAP (#OPT) | 0.00% (10) | 0.00% (10) | 0.00% (10) | 0.03% (1) | 0.08% (0)
Speedup factor | 46.4x | 13.8x | 65.2x | 1.1x | 1.0x

As can be seen, our valid inequalities only increase the total number of
constraints by 1.0-1.6%. However, the LR-GAP is significantly reduced compared
to that obtained by BN. This effect is particularly noticeable for the IEEE-
RTS-73 system with an average value of the linear relaxation gap equal to
0.0001%, meaning that the linear relaxation of problem (9) with the proposed
valid inequalities is very tight and its solution very close to the actual
solution of that problem. Furthermore, the inclusion of the valid inequalities
lead to average speedup factors that range between 13.8x and 65.2x for the
three smallest systems. For the two largest systems, the time limit is reached
in most instances, but the average MIP-GAP is reduced to 0.03% and 0.08%,
respectively.

Finally, we present similar simulation results for the setup in which the
valid inequalities are also used to boost the tightening and screening power
of Algorithm 1 (that is, the valid inequalities are also included in problem
(3)), leading to method TS+V. Table 6 provides the average number of
constraints, the average values of LR-GAP and MIP-GAP, and the average speedup
factor of TS+V with respect to BN in Table 2. Since increasing the number of
iterations of Algorithm 1 barely affects the performance of TS+V, all the data
shown in Table 6 correspond to $\kappa=1$.

Table 6: Tightening and screening with valid inequalities (TS+V): Results | IEEE-RTS-24 | IEEE-57 | IEEE-RTS-73 | IEEE-118 | IEEE-300
---|---|---|---|---|---
#CON | 3.49% | 1.61% | 3.84% | 2.68% | 3.30%
LR-GAP | 0.3365% | 0.1519% | 0.0001% | 0.2821% | 0.1603%
MIP-GAP (#OPT) | 0.00% (10) | 0.00% (10) | 0.00% (10) | 0.00% (10) | 0.00% (10)
Speedup factor | 706.8x | 35.3x | 1470.0x | 23.1x | 8.5x

The comparison of the results in Tables 4, 5 and 6 yields the following
observations. First, including the valid inequalities in Algorithm 1
strengthens the Big-Ms even further, which in turn remarkably reduces the
linear relaxation gap and increases the number of constraints identified as
redundant in model (9). Indeed, the total number of constraints eventually
retained by TS+V ranges between 1.61% and 3.84% if compared with BN. Second,
TS+V can solve the 10 instances to global optimality in less than 10 hours for
the five test systems considered in these numerical experiments. In fact, the
optimal solutions obtained by TS+V are the ones we use to compute the values
of the linear relaxation gap throughout these simulations. Third, TS+V is able
to achieve speedup factors between 8.5x and 1470.0x depending on the test
system. All in all, TS+V features the best computational performance in terms
of resolution time and MIP-GAP among all the methods tested so far.

Next, in Table 7 the results of TS+V are contrasted with those provided by
state-of-the-art approximations available in the literature. In particular, we
consider the CVaR-based and ALSO-X conservative approximations, both described
in Jiang and Xie (2022). Table 7 provides the average cost increase in
percentage with respect to the optimal cost and the speedup factor with regard
to BN for the different methodologies compared. As expected, the CVaR-based
approximation leads to conservative results and involves average cost
increases that range between 0.49% and 2.88%. Interestingly enough, TS+V
computes the optimal solution and involves a higher speedup factor than the
CVaR-based approach for two of the five systems. Compared with the two ALSO-X
approximations, TS+V obtains the global optimal solution in all cases with
speedup factors that are still tantamount to those of these approximate
methods.

Table 7: Comparison of the proposed TS+V approach and existing approximate methods | IEEE-RTS-24 | IEEE-57 | IEEE-RTS-73 | IEEE-118 | IEEE-300
---|---|---|---|---|---
Average cost increase | TS+V | 0.00% | 0.00% | 0.00% | 0.00% | 0.00%
CVaR | 2.88% | 0.53% | 1.71% | 0.57% | 0.49%
ALSO-X | 0.80% | 0.08% | 0.41% | 0.08% | 0.05%
ALSO-X+ | 0.53% | 0.07% | 0.12% | 0.05% | 0.04%
Speedup factor | TS+V | 706.8x | 35.3x | 1470.0x | 23.1x | 8.5x
CVaR | 387.1x | 117.0x | 265.1x | 4045.5x | 779.7x
ALSO-X | 18.0x | 1.2x | 21.3x | 148.6x | 11.13x
ALSO-X+ | 7.8x | 0.7x | 6.6x | 49.3x | 3.7x

To conclude this study, we compare our solution approach with the one proposed
in Luedtke (2014). This author presents an exact approach to solve CCPs based
on the addition of quantile cuts in a so-called lazy fashion to a relaxation
of the problem that only includes the deterministic constraints. In this way,
at each node of the branching tree we seek for a violated constraint from a
scenario $s$ with $y_{s}=0$. Using the coefficients of this violated
constraint, a valid inequality is derived by means of the resolution of
$|\mathcal{S}|$ linear separation subproblems. These strengthened cuts contain
binary variables associated to several scenarios and are obtained by applying
the star inequalities of Atamtürk et al. (2000).

We have implemented the procedure including, at each iteration, both the
initial or base violated inequality and the strengthened one, since this
approach delivered the best results. The comparison between our approach and
Luedtke’s (which we name BCD for short, from _branch-and-cut decomposition_)
is shown in Table 8. As can be seen, a similar number of constraints are
generated in the BCD approach (that includes the deterministic constraints and
the cuts generated dynamically). Furthermore, our procedure clearly
outperforms the BCD approach in terms of computational time. One possible
explanation, already pointed out in the paper, is that the generation of the
cuts requires solving 1000 subproblems at each iteration. In the particular
application presented in Luedtke (2014), however, this potential bottleneck is
overcome by leveraging the specific structure of the subproblem to develop a
closed-form expression of its solution, thus avoiding the need to optimize
over all the scenario sets. With that said, we remark that the times provided
in Table 8 under the acronym BCD-P correspond to a _parallel_ implementation
of the BCD algorithm, in which only the time required to solve the most time-
consuming subproblem out of the 1000 to be solved at each iteration is added
to the final solution time reported. For completeness and because the
_parallel_ implementation of BCD using off-the-shelf optimization solvers is
not trivial at all, we also report in Table 8 the speedup factor that is
attained by a _serial_ implementation of BCD, termed BCD-S in the table. As
can be seen, this implementation is, however, far from being competitive in
all cases.

Table 8: Comparison of the proposed TS+V approach and the BCD approach | IEEE-RTS-24 | IEEE-57 | IEEE-RTS-73 | IEEE-118 | IEEE-300
---|---|---|---|---|---
#CONS | TS+V | 3.49% | 1.61% | 3.84% | 2.68% | 3.30%
BCD | 11.37% | 2.69% | 3.58% | 2.46% | 1.53%
Speedup factor | TS+V | 706.8x | 35.3x | 1470.0x | 23.1x | 8.5x
BCD-P | 10.12x | 2.05x | 69.00x | 9.24x | 4.74x
BCD-S | 0.12x | 0.04x | 0.30x | 1.45x | 0.42x

##  4 Conclusions and future research

In this paper, we propose a novel exact resolution technique for a MIP SAA-
based reformulation of joint chance-constrained problems in the form of (1).
Our methodology includes a screening method to eliminate superfluous
constraints based on an iterative procedure to repeatedly tighten the Big-Ms
present in the MIP. These procedures are combined with the addition of
inequalities that are valid provided that the technology matrix in the chance-
constraint system exhibits the structure presented in (4). Said inequalities
strengthen the linear relaxation of the MIP SAA-based reformulation and allow
for additional screening of constraints. The resultant model is thus compact
and tight.

We have applied our methodology to solve the joint chance-constrained DC
Optimal Power Flow. In the case study, we show that, in comparison with the
benchmark model, our methodology provides remarkable results in terms of the
linear relaxation bounds, the RAM memory needed to solve the instances, and
the total computational resolution time. Specifically, our method TS+V solves
to optimality all of the instances generated for the IEEE-RTS-118 and the
IEEE-RTS-300 test systems, the majority of which are not solved within 10
hours of computational time using the benchmark approach. Furthermore, the
average number of constraints eliminated from all instances with TS+V always
exceeds a 95% of them, and the lower bound is markedly increased by the
inclusion of the valid inequalities, showing the outstanding results of the
combination of the methods developed.

The comparison of our results with those provided by existing approximate and
exact methods shows that our approach is computationally very competitive for
small and medium-sized instances, always providing the best results in terms
of cost. For the large instances addressed, while outperformed by the
approximate methods in terms of computational time (as expected), our exact
solution strategy not only provides a certificate of optimality, but also
returns the optimal solution within the set time limit. A promising future
research line consists in the development of a generalized set of valid
inequalities that combine variables from pairs or subgroups of constraints
$j\in\mathcal{J}$ in the chance-constrained system (1c).

## Acknowledgements

This work was supported in part by the European Research Council (ERC) under
the EU Horizon 2020 research and innovation program (grant agreement No.
755705), in part by the Spanish Ministry of Science and Innovation
(AEI/10.13039/501100011033) through project PID2020-115460GB-I00, and in part
by the Junta de Andalucía (JA) and the European Regional Development Fund
(FEDER) through the research project P20_00153. Á. Porras is also financially
supported by the Spanish Ministry of Science, Innovation and Universities
through the University Teacher Training Program with fellowship number
FPU19/03053. Finally, the authors thankfully acknowledge the computer
resources, technical expertise, and assistance provided by the SCBI
(Supercomputing and Bioinformatics) center of the University of Málaga.

## References

  * Abdi and Fukasawa (2016) Abdi A, Fukasawa R (2016) On the mixing set with a knapsack constraint. _Mathematical Programming_ 157(1):191–217.
  * Agarwal et al. (1998) Agarwal PK, Aronov B, Chan TM, Sharir M (1998) On levels in arrangements of lines, segments, planes, and triangles. _Discrete & Computational Geometry_ 19(3):315–331.
  * Agarwal and Matoušek (1995) Agarwal PK, Matoušek J (1995) Dynamic half-space range reporting and its applications. _Algorithmica_ 13(4):325–345.
  * Ahmed et al. (2017) Ahmed S, Luedtke J, Song Y, Xie W (2017) Nonanticipative duality, relaxations, and formulations for chance-constrained stochastic programs. _Mathematical Programming_ 162(1):51–81.
  * Ahmed and Xie (2018) Ahmed S, Xie W (2018) Relaxations and approximations of chance constraints under finite distributions. _Mathematical Programming_ 170(1):43–65.
  * Atamtürk et al. (2000) Atamtürk A, Nemhauser GL, Savelsbergh MWP (2000) Conflict graphs in solving integer programming problems. _European Journal of Operational Research_ 121(1):40–55, ISSN 0377-2217, URL http://dx.doi.org/10.1016/S0377-2217(99)00015-6.
  * Ben-Tal and Nemirovski (2000) Ben-Tal A, Nemirovski A (2000) Robust solutions of linear programming problems contaminated with uncertain data. _Mathematical programming_ 88(3):411–424.
  * Bienstock et al. (2014) Bienstock D, Chertkov M, Harnett S (2014) Chance-constrained optimal power flow: Risk-aware network control under uncertainty. _SIAM Review_ 56(3):461–495.
  * Borobia (1994) Borobia A (1994) Mirror property for nonsingular mixed configurations of lines and points in $\mathbb{R}$3\. _Discrete & Computational Geometry_ 11(3):311–320.
  * Chan et al. (1997) Chan TM, Snoeyink J, Yap CK (1997) Primal dividing and dual pruning: Output-sensitive construction of four-dimensional polytopes and three-dimensional Voronoi diagrams. _Discrete & Computational Geometry_ 18(4):433–454.
  * Chazelle and Edelsbrunner (1987) Chazelle B, Edelsbrunner H (1987) An improved algorithm for constructing kth-order Voronoi diagrams. _IEEE Transactions on Computers_ 100(11):1349–1354.
  * Chazelle and Preparata (1986) Chazelle B, Preparata FP (1986) Halfspace range search: An algorithmic application of k-sets. _Discrete & Computational Geometry_ 1(1):83–93.
  * Cheema et al. (2014) Cheema MA, Shen Z, Lin X, Zhang W (2014) A unified framework for efficiently processing ranking related queries. _EDBT_ , 427–438.
  * Chen et al. (2021) Chen G, Zhang H, Hui H, Song Y (2021) Time-efficient joint chance-constrained optimal power flow with a learning-based robust approximation. _arXiv preprint arXiv:2112.09827_ .
  * Clarkson (1987) Clarkson KL (1987) New applications of random sampling in computational geometry. _Discrete & Computational Geometry_ 2(2):195–222.
  * Conforti et al. (2014) Conforti M, Cornuéjols G, Zambelli G (2014) _Integer Programming_ , volume 271 (Springer).
  * Daníelsson et al. (2008) Daníelsson J, Jorgensen BN, de Vries CG, Yang X (2008) Optimal portfolio allocation under the probabilistic VaR constraint and incentives for financial innovation. _Annals of Finance_ 4(3):345–367.
  * Das et al. (2007) Das G, Gunopulos D, Koudas N, Sarkas N (2007) Ad-hoc top-k query answering for data streams. _Proceedings of the 33rd International Conference on Very Large Data Bases_ , 183–194.
  * Dentcheva et al. (2000) Dentcheva D, Prékopa A, Ruszczynski A (2000) Concavity and efficient points of discrete distributions in probabilistic programming. _Mathematical Programming_ 89(1):55–77.
  * Dey (1998) Dey TK (1998) Improved bounds for planar k-sets and related problems. _Discrete & Computational Geometry_ 19(3):373–382.
  * Edelsbrunner and Welzl (1985) Edelsbrunner H, Welzl E (1985) On the number of line separations of a finite set in the plane. _Journal of Combinatorial Theory, Series A_ 38(1):15–29.
  * Edelsbrunner and Welzl (1986) Edelsbrunner H, Welzl E (1986) Constructing belts in two-dimensional arrangements with applications. _SIAM Journal on Computing_ 15(1).
  * Elçi et al. (2018) Elçi Ö, Noyan N, Bülbül K (2018) Chance-constrained stochastic programming under variable reliability levels with an application to humanitarian relief network design. _Computers & Operations Research_ 96:91–107.
  * Esteban-Pérez and Morales (2023) Esteban-Pérez A, Morales JM (2023) Distributionally robust optimal power flow with contextual information. _European Journal of Operational Research_ 306(3):1047–1058.
  * Frank et al. (2012) Frank S, Steponavice I, Rebennack S (2012) Optimal power flow: A bibliographic survey I. _Energy Systems_ 3(3):221–258.
  * Günlük and Pochet (2001) Günlük O, Pochet Y (2001) Mixing mixed-integer inequalities. _Mathematical Programming_ 90(3):429–457.
  * Gurobi Optimization, LLC (2022) Gurobi Optimization, LLC (2022) Gurobi Optimizer Reference Manual. URL https://www.gurobi.com.
  * Halperin and Sharir (2017) Halperin D, Sharir M (2017) Arrangements. _Handbook of Discrete and Computational Geometry_ , 723–762 (Chapman and Hall/CRC), ISBN 978-1-315-11960-1, URL http://dx.doi.org/10.1201/9781315119601-28.
  * Henrion (2007) Henrion R (2007) Structural properties of linear probabilistic constraints. _Optimization_ 56(4):425–440.
  * Henrion and Strugarek (2011) Henrion R, Strugarek C (2011) Convexity of chance constraints with dependent random variables: The use of copulae. _Stochastic Optimization Methods in Finance and Energy_ , 427–439 (Springer).
  * Hong et al. (2011) Hong LJ, Yang Y, Zhang L (2011) Sequential convex approximations to joint chance constrained programs: A Monte Carlo approach. _Operations Research_ 59(3):617–630.
  * Hou and Roald (2020) Hou AM, Roald LA (2020) Chance constraint tuning for optimal power flow. _2020 International Conference on Probabilistic Methods Applied to Power Systems (PMAPS)_ , 1–6, URL http://dx.doi.org/10.1109/PMAPS47429.2020.9183552.
  * Jiang and Xie (2022) Jiang N, Xie W (2022) ALSO-X and ALSO-X+: Better convex approximations for chance constrained programs. _Operations Research_ 70(6):3581–3600.
  * Küçükyavuz and Jiang (2022) Küçükyavuz S, Jiang R (2022) Chance-constrained optimization under limited distributional information: A review of reformulations based on sampling and distributional robustness. _EURO Journal on Computational Optimization_ 10:100030, ISSN 2192-4406, URL http://dx.doi.org/https://doi.org/10.1016/j.ejco.2022.100030.
  * Lagoa et al. (2005) Lagoa CM, Li X, Sznaier M (2005) Probabilistically constrained linear programs and risk-adjusted controller design. _SIAM Journal on Optimization_ 15(3):938–951.
  * Lejeune and Dehghanian (2020) Lejeune MA, Dehghanian P (2020) Optimal power flow models with probabilistic guarantees: A boolean approach. _IEEE Transactions on Power Systems_ 35(6):4932–4935.
  * Luedtke (2014) Luedtke J (2014) A branch-and-cut decomposition algorithm for solving chance-constrained mathematical programs with finite support. _Mathematical Programming_ 146(1-2):219–244, ISSN 0025-5610, 1436-4646, URL http://dx.doi.org/10.1007/s10107-013-0684-6.
  * Luedtke and Ahmed (2008) Luedtke J, Ahmed S (2008) A sample approximation approach for optimization with probabilistic constraints. _SIAM Journal on Optimization_ 19(2):674–699.
  * Luedtke et al. (2010) Luedtke J, Ahmed S, Nemhauser GL (2010) An integer programming approach for linear programs with probabilistic constraints. _Mathematical programming_ 122(2):247–272.
  * Miller and Wagner (1965) Miller BL, Wagner HM (1965) Chance constrained programming with joint constraints. _Operations Research_ 13(6):930–945.
  * Nair and Miller-Hooks (2011) Nair R, Miller-Hooks E (2011) Fleet management for vehicle sharing operations. _Transportation Science_ 45(4):524–540.
  * Najjarbashi and Lim (2020) Najjarbashi A, Lim GJ (2020) A decomposition algorithm for the two-stage chance-constrained operating room scheduling problem. _IEEE Access_ 8:80160–80172.
  * Natarajan et al. (2008) Natarajan K, Pachamanova D, Sim M (2008) Incorporating asymmetric distributional information in robust value-at-risk optimization. _Management Science_ 54(3):573–585.
  * Nemirovski and Shapiro (2006) Nemirovski A, Shapiro A (2006) Scenario approximations of chance constraints. _Probabilistic and randomized methods for design under uncertainty_ 3–47.
  * Nemirovski and Shapiro (2007) Nemirovski A, Shapiro A (2007) Convex approximations of chance constrained programs. _SIAM Journal on Optimization_ 17(4):969–996.
  * Nivasch (2008) Nivasch G (2008) An improved, simple construction of many halving edges. _Contemporary Mathematics_ 453:299–306.
  * OASYS (2023) OASYS (2023) Data and Code for a Tight and Compact Model of the SAA-based Joint Chance-constrained OPF. _GitHub repository (https://github.com/groupoasys/TC_SAA_JCC-OPF)_ URL https://github.com/groupoasys/TC_SAA_JCC-OPF.
  * Peña-Ordieres et al. (2020) Peña-Ordieres A, Luedtke JR, Wächter A (2020) Solving chance-constrained problems via a smooth sample-based nonlinear approximation. _SIAM Journal on Optimization_ 30(3):2221–2250, URL http://dx.doi.org/10.1137/19M1261985.
  * Peña-Ordieres et al. (2021) Peña-Ordieres A, Molzahn DK, Roald LA, Wächter A (2021) DC optimal power flow with joint chance constraints. _IEEE Transactions on Power Systems_ 36(1):147–158, URL http://dx.doi.org/10.1109/TPWRS.2020.3004023.
  * Power Grid Lib (2022) Power Grid Lib (2022) GitHub repository, available at: (https://github.com/power-grid-lib/pglib-opf).
  * Prékopa (2003) Prékopa A (2003) Probabilistic programming. _Handbooks in Operations Research and Management Science_ 10:267–351.
  * Preparata and Shamos (2012) Preparata FP, Shamos MI (2012) _Computational Geometry: An Introduction_ (Springer Science & Business Media).
  * Qiu et al. (2014) Qiu F, Ahmed S, Dey SS, Wolsey LA (2014) Covering linear programming with violations. _INFORMS Journal on Computing_ 26(3):531–546, URL http://dx.doi.org/10.1287/ijoc.2013.0582.
  * Ramos (1999) Ramos EA (1999) On range reporting, ray shooting and k-level construction. _Proceedings of the Fifteenth Annual Symposium on Computational Geometry_ , 390–399.
  * Roald and Andersson (2018) Roald L, Andersson G (2018) Chance-constrained AC optimal power flow: Reformulations and efficient algorithms. _IEEE Transactions on Power Systems_ 33(3):2906–2918, URL http://dx.doi.org/10.1109/TPWRS.2017.2745410.
  * Rockafellar and Uryasev (2000) Rockafellar RT, Uryasev S (2000) Optimization of conditional value-at-risk. _The Journal of Risk_ 2(3):21–41, ISSN 14651211, URL http://dx.doi.org/10.21314/JOR.2000.038.
  * Roos and Widmayer (1994) Roos T, Widmayer P (1994) K-Violation linear programming. _Information Processing Letters_ 52(2):109–114, ISSN 0020-0190, URL http://dx.doi.org/10.1016/0020-0190(94)00134-0.
  * Shapiro (2003) Shapiro A (2003) Monte Carlo sampling methods. _Handbooks in operations research and management science_ 10:353–425.
  * Sharir (2011) Sharir M (2011) An improved bound for k-sets in four dimensions. _Combinatorics, Probability and Computing_ 20(1):119–129.
  * Sharir et al. (2000) Sharir M, Smorodinsky S, Tardos G (2000) An improved bound for k-sets in three dimensions. _Proceedings of the Sixteenth Annual Symposium on Computational Geometry_ , 43–49.
  * Shaw et al. (2016) Shaw K, Irfan M, Shankar R, Yadav SS (2016) Low carbon chance constrained supply chain network design problem: A Benders decomposition based approach. _Computers & Industrial Engineering_ 98:483–497.
  * Shen et al. (2012) Shen Z, Cheema MA, Lin X (2012) Loyalty-based selection: Retrieving objects that persistently satisfy criteria. _Proceedings of the 21st ACM International Conference on Information and Knowledge Management_ , 2189–2193.
  * Song et al. (2014) Song Y, Luedtke JR, Küçükyavuz S (2014) Chance-constrained binary packing problems. _INFORMS Journal on Computing_ 26(4):735–747.
  * Sun et al. (2014) Sun H, Xu H, Wang Y (2014) Asymptotic analysis of sample average approximation for stochastic optimization problems with joint chance constraints via conditional value at risk and difference of convex functions. _Journal of Optimization Theory and Applications_ 161(1):257–284, ISSN 0022-3239, 1573-2878, URL http://dx.doi.org/10.1007/s10957-012-0127-1.
  * Taleizadeh et al. (2012) Taleizadeh AA, Niaki STA, Makui A (2012) Multiproduct multiple-buyer single-vendor supply chain problem with stochastic demand, variable lead-time, and multi-chance constraint. _Expert Systems with Applications_ 39(5):5338–5348.
  * Tanner et al. (2008) Tanner MW, Sattenspiel L, Ntaimo L (2008) Finding optimal vaccination strategies under parameter uncertainty using stochastic programming. _Mathematical Biosciences_ 215(2):144–151, ISSN 0025-5564, URL http://dx.doi.org/10.1016/j.mbs.2008.07.006.
  * Tayur et al. (1995) Tayur SR, Thomas RR, Natraj NR (1995) An algebraic geometry algorithm for scheduling in presence of setups and correlated demands. _Mathematical Programming_ 69(1):369–401.
  * Tóth et al. (2017) Tóth CD, O’Rourke J, Goodman JE (2017) _Handbook of Discrete and Computational Geometry_ (CRC press).
  * Tóth (2000) Tóth G (2000) Point sets with many k-sets. _Proceedings of the Sixteenth Annual Symposium on Computational Geometry_ , 37–42.
  * Van Ackooij et al. (2011) Van Ackooij W, Zorgati R, Henrion R, Möller A (2011) Chance constrained programming and its applications to energy management. _Stochastic Optimization-Seeing the Optimal for the Uncertain_ 291–320.
  * Vrakopoulou et al. (2013) Vrakopoulou M, Margellos K, Lygeros J, Andersson G (2013) Probabilistic Guarantees for the N-1 Security of Systems with Wind Power Generation. _Reliability and Risk Evaluation of Wind Integrated Power Systems_ , 59–73 (Springer).
  * Xie and Ahmed (2016) Xie W, Ahmed S (2016) On the quantile cut closure of chance-constrained problems. _International Conference on Integer Programming and Combinatorial Optimization_ , 398–409 (Springer).
  * Xie and Ahmed (2018) Xie W, Ahmed S (2018) On quantile cuts and their closure for chance constrained optimization problems. _Mathematical Programming_ 172(1):621–646.
  * Yu et al. (2012) Yu A, Agarwal PK, Yang J (2012) Processing a large number of continuous preference top-k queries. _Proceedings of the 2012 ACM SIGMOD International Conference on Management of Data_ , 397–408.
  * Zhang et al. (2015) Zhang Y, Shen S, Mathieu JL (2015) Data-driven optimization approaches for optimal power flow with uncertain reserves from load control. _2015 American Control Conference (ACC)_ , 3013–3018 (IEEE).

##  Appendix Appendix A Relationship with _quantile cuts_

Consider problem (1), with feasible region $\mathcal{FR}$, assume $X$ is a
non-empty and compact set given by deterministic constraints and
$a_{j}(\omega)=a^{0}_{j}+\Omega_{j}(\omega)\hat{a}_{j}$, $j\in\mathcal{J}$.
For a given scenario $s\in\mathcal{S}$, define $\mathcal{X}_{s}$ as the set of
points that satisfy the constraints for that scenario, i.e.,

$\mathcal{X}_{s}:=\left\\{x\in\mathbb{R}^{|\mathcal{I}|}:\Omega_{js}\hat{a}_{j}^{\top}x\leq
b_{js}-x^{\top}a_{j}^{0},\quad\forall j\in\mathcal{J}\right\\}.$

To characterize the feasible region $\mathcal{FR}$, let
$\mathcal{S}^{F}:=\left\\{S^{\prime}\subseteq\mathcal{S}:|S^{\prime}|\geq|\mathcal{S}|-p\right\\}$
be the collection of all feasible subsets of scenarios. Then $\mathcal{FR}$
can be represented as:

$\mathcal{FR}:=\bigcup_{S^{\prime}\in\mathcal{S}^{F}}\left[X\bigcap_{s\in
S^{\prime}}\mathcal{X}_{s}\right].$

In Xie and Ahmed (2018) it is proved that $\mathcal{FR}$ can also be written
in a conjunctive normal form. To do so, we introduce the definition of a
_minimal partial covering subset_.

######  Definition A.1.

A set $S^{\prime\prime}\subseteq\mathcal{S}$ is a _partial covering subset_ if
it intersects with all the feasible scenario subsets in $\mathcal{S}^{F}$,
i.e., for any $S^{\prime}\in\mathcal{S}^{F}$ we have $S^{\prime}\cap
S^{\prime\prime}\neq\emptyset$. Also, a set $S^{\prime\prime}$ is a _minimal_
partial covering subset if there does not exist another partial covering
subset $\hat{S}\subseteq\mathcal{S}$ such that $\hat{S}\subsetneq S^{\prime}$.
Let $\mathcal{S}^{C}$ denote the collection of all minimal partial covering
subsets. When the samples are independent and identically distributed, like in
our case,
$\mathcal{S}^{C}:=\left\\{S^{\prime\prime}\subseteq\mathcal{S}:|S^{\prime\prime}|=p+1\right\\}$.

######  Proposition A.2 ((Xie and Ahmed 2018, Proposition 3)).

$\mathcal{FR}=\bigcap_{S^{\prime\prime}\in\mathcal{S}^{C}}\left[\bigcup_{s\in
S^{\prime\prime}}\left(X\bigcap\mathcal{X}_{s}\right)\right].$

Consider now the relaxation of problem (1) associated to a given row
$j\in\mathcal{J}$:

$\displaystyle\min_{x}\quad$ $\displaystyle f(x)$ (10a) s.t.
$\displaystyle\mathbb{P}\left\\{a_{j}(\omega)^{\top}x\leq
b_{j}(\omega)\right\\}\geq 1-\epsilon$ (10b) $\displaystyle x\in X$ (10c)

and its feasible region $\mathcal{FR}_{j}$. Clearly, it holds that
$\mathcal{FR}\subseteq\bigcap_{j\in\mathcal{J}}\mathcal{FR}_{j}$. This is
equivalent to considering the problem without joint constraints. In a similar
manner, for a given scenario $s\in\mathcal{S}$ and row $j\in\mathcal{J}$,
define $\mathcal{X}_{js}$ as the set of points that satisfy the scenario for
that row:

$\mathcal{X}_{js}:=\left\\{x\in\mathbb{R}^{|\mathcal{I}|}:\Omega_{js}\hat{a}_{j}^{\top}x\leq
b_{js}-x^{\top}a_{j}^{0}\right\\},$

with $\mathcal{X}_{s}=\bigcap_{j\in\mathcal{J}}\mathcal{X}_{js}$. Making use
of the previous results, we have

$\mathcal{FR}_{j}:=\bigcup_{S^{\prime}\in\mathcal{S}^{F}}\left[X\bigcap_{s\in
S^{\prime}}\mathcal{X}_{js}\right]=\bigcap_{S^{\prime\prime}\in\mathcal{S}^{C}}\left[\bigcup_{s\in
S^{\prime\prime}}\left(X\bigcap\mathcal{X}_{js}\right)\right].$

Using this characterization of a relaxation of the feasible region, we are
able to prove the following:

######  Proposition A.3.

Let $\mathcal{U}_{j}:=\left\\{x\in
X:U_{j}^{p+1}(\hat{a}_{j}^{\top}x)\leq-x^{\top}a_{j}^{0}\right\\}$. Then for
$j\in\mathcal{J}$, $\mathcal{U}_{j}=\mathcal{FR}_{j}$.

###### Proof.

$[\mathcal{U}_{j}\subseteq\mathcal{FR}_{j}]$ Assume $x\in\mathcal{U}_{j}$.
Then $x\in X$ satisfies the chance constraint for (at least) $|\mathcal{S}|-p$
scenarios. But then, for any minimal partial covering
$S^{\prime\prime}\in\mathcal{S}^{C}$, since $|S^{\prime\prime}|=p+1$, we must
have $x\in\bigcup_{s\in S^{\prime\prime}}(X\bigcap\mathcal{X}_{js})$. Given
that this occurs for any $S^{\prime\prime}$, then $x$ belongs to the
intersection of $S^{\prime\prime}$, and hence $x\in\mathcal{FR}_{j}$.

$[\mathcal{FR}_{j}\subseteq\mathcal{U}_{j}]$ This proof is straightforward
considering the validity of the proof of Proposition 2.3, but assuming
$x\in\mathcal{FR}_{j}$ instead of $x\in\mathcal{FR}$. Thus, if
$x\in\mathcal{FR}_{j}$, then $x$ satisfies the single chance constraint $j$ by
definition. Since
$\mathcal{U}_{j}^{p+1}(\hat{a}_{j}^{\top}x)\leq-x^{\top}a_{j}^{0}$ is a valid
inequality for problem (10) associated to $j$, $x$ must satisfy it as well. ∎

Now, we are going to establish the relationship between our valid inequalities
and the quantile closure described in Xie and Ahmed (2018). Let us first
introduce the necessary definitions to relate the quantile closure with
$\mathrm{conv}(\mathcal{FR})$.

######  Definition A.4.

Given $\alpha\in\mathbb{R}^{|\mathcal{I}|}$, let
$\\{\gamma^{\alpha}_{s}(X)\\}$ be the optimal values of

$\gamma^{\alpha}_{s}(X):=\min\left\\{\alpha^{\top}x:x\in
X\cap\mathcal{X}_{s}\right\\}\quad\forall s\in\mathcal{S}.$

The $(1-\epsilon)$– _quantile_ (or simply _quantile_) of
$\\{\gamma^{\alpha}_{s}(X)\\}_{s\in S}$ is denoted by $\gamma^{\alpha}_{q}(X)$
and is given by

$\gamma^{\alpha}_{q}(X):=\min_{S^{\prime}\in\mathcal{S}^{F}}\max_{s\in
S^{\prime}}\gamma^{\alpha}_{s}(X)$

and the associated _quantile cut_ is

$\alpha^{\top}x\geq\gamma^{\alpha}_{q}(X).$

Note that the definition of quantile cuts depends on $X$, a set of
deterministic constraints that we may change if we apply cuts successively.

######  Definition A.5.

The _first quantile closure_ of $X$ is defined as

$X^{1}:=\bigcap_{\alpha\in\mathbb{R}^{|\mathcal{I}|}}\left\\{x\in\mathbb{R}^{|\mathcal{I}|}:\alpha^{\top}x\geq\gamma_{q}^{\alpha}(X)\right\\}.$

Inductively, we define the _$r$ th round quantile closure_ $X^{r}$ as:

$X^{r}:=\bigcap_{\alpha\in\mathbb{R}^{|\mathcal{I}|}}\left\\{x\in\mathbb{R}^{|\mathcal{I}|}:\alpha^{\top}x\geq\gamma_{q}^{\alpha}(X^{r-1})\right\\}\quad
r\geq 2.$

######  Theorem A.6 ((Xie and Ahmed 2018, Theorem 1)).

$X^{1}=\bigcap_{S^{\prime\prime}\in\mathcal{S}^{C}}\mathrm{conv}\left[\bigcup_{s\in
S^{\prime\prime}}\left(X\bigcap\mathcal{X}_{s}\right)\right],$ (11)

and for each $r\geq 2$, $r\in\mathbb{Z}_{++}$,

$X^{r}=\bigcap_{S^{\prime\prime}\in\mathcal{S}^{C}}\mathrm{conv}\left[\bigcup_{s\in
S^{\prime\prime}}\left(X^{r-1}\bigcap\mathcal{X}_{s}\right)\right],$ (12)

where $X^{0}=X$.

######  Theorem A.7 ((Xie and Ahmed 2018, Theorem 2)).

The set sequence $\\{X^{r}\\}$ converges to $\mathrm{conv}(\mathcal{FR})$ with
respect to the Hausdorff distance; i.e.,
$\widebar{X}:=\lim_{r\rightarrow\infty}X^{r}=\mathrm{conv}(\mathcal{FR})$.

If, like previously, we consider these definitions applied to the relaxation
of problem (1) that results from considering a single row $j\in\mathcal{J}$,
then we can establish the following results:

######  Definition A.8.

Given $\alpha\in\mathbb{R}^{|\mathcal{I}|}$, let
$\\{\gamma^{\alpha}_{js}(X)\\}$ be the optimal values of

$\gamma^{\alpha}_{js}(X)=\min\left\\{\alpha^{\top}x:x\in
X\cap\mathcal{X}_{js}\right\\}\quad\forall s\in\mathcal{S}.$

The $(1-\epsilon)$– _quantile_ (or simply _quantile_) of
$\\{\gamma^{\alpha}_{js}(X)\\}_{s\in S}$ is denoted by
$\gamma^{\alpha}_{jq}(X)$ and is given by

$\gamma^{\alpha}_{jq}(X):=\min_{S^{\prime}\in\mathcal{S}^{F}}\max_{s\in
S^{\prime}}\gamma^{\alpha}_{js}(X),$

and the associated _quantile cut_ is

$\alpha^{\top}x\geq\gamma^{\alpha}_{jq}(X).$

######  Definition A.9.

The _first quantile closure_ of $X$ associated to $j\in\mathcal{J}$ is defined
as

$X^{1}_{j}:=\bigcap_{\alpha\in\mathbb{R}^{|\mathcal{I}|}}\left\\{x\in\mathbb{R}^{|\mathcal{I}|}:\alpha^{\top}x\geq\gamma_{jq}^{\alpha}(X)\right\\}.$

Inductively, we define the _$r$ th round quantile closure_ $X^{r}_{j}$
associated to $j\in\mathcal{J}$ as:

$X^{r}_{j}:=\bigcap_{\alpha\in\mathbb{R}^{|\mathcal{I}|}}\left\\{x\in\mathbb{R}^{|\mathcal{I}|}:\alpha^{\top}x\geq\gamma_{jq}^{\alpha}(X^{r-1})\right\\}\quad
r\geq 2.$

######  Theorem A.10 ((Xie and Ahmed 2018, Theorem 1)).

$X^{1}_{j}=\bigcap_{S^{\prime\prime}\in\mathcal{S}^{C}}\mathrm{conv}\left[\bigcup_{s\in
S^{\prime\prime}}\left(X\bigcap\mathcal{X}_{js}\right)\right],$ (13)

and for each $r\geq 2$, $r\in\mathbb{Z}_{++}$,

$X^{r}_{j}=\bigcap_{S^{\prime\prime}\in\mathcal{S}^{C}}\mathrm{conv}\left[\bigcup_{s\in
S^{\prime\prime}}\left(X^{r-1}_{j}\bigcap\mathcal{X}_{js}\right)\right],$ (14)

where $X^{0}_{j}=X$.

######  Theorem A.11 ((Xie and Ahmed 2018, Theorem 2)).

The set sequence $\\{X^{r}_{j}\\}$ converges to
$\mathrm{conv}(\mathcal{FR}_{j})$ with respect to the Hausdorff distance;
i.e.,
$\widebar{X}_{j}:=\lim_{r\rightarrow\infty}X^{r}_{j}=\mathrm{conv}(\mathcal{FR}_{j})$.

######  Corollary A.12.

$\mathcal{U}_{j}\subseteq\mathrm{conv}(\mathcal{U}_{j})=\widebar{X}_{j}$ (15)

This result is particularly relevant because, as it is proven in the same
work, the separation over the first quantile closure is, in general, NP-hard.
Therefore, Corollary A.12 gives a way of deriving
$\mathrm{conv}(\mathcal{FR}_{j})$ through the explicit set of valid
inequalities (6).

##  Appendix Appendix B Complexity and extension to the multidimensional case

Suppose there exist matrices $A_{j}\in\mathbb{R}^{|\mathcal{I}|\times
q_{j}},q_{j}\leq|\mathcal{I}|$ (this includes the identity),
$j\in\mathcal{J}$, such that $a_{j}(\omega)$ in the chance constraint can be
written as $a_{j}(\omega)=a^{0}_{j}+A_{j}\Omega_{j}(\omega)$,
$j\in\mathcal{J}$, with $a_{j}^{0}\in\mathbb{R}^{|\mathcal{I}|}$ and
$\Omega_{j}(\omega)$ being a $q_{j}$-dimensional vector of real-valued
functions whose domains include the support of $\omega$. Note that the case
described in Subsection 2.2 corresponds to $q_{j}=1,\forall j\in\mathcal{J}$.

Now, define $z_{j}:=A_{j}^{\top}x$, with $z_{j}^{d}$ and $z_{j}^{u}$ being
vectors containing the component-wise lower and upper bounds on $z$ induced by
the feasibility set (4b). Then, we can define the sets of hyperplanes
$\Pi_{j}$ for $j\in\mathcal{J}$

$\Pi_{j}:=\left\\{f_{js}(z_{j})=\Omega_{js}^{\top}z_{j}-b_{js},\quad\text{for
}z_{j}\in B(z_{j}^{d},z_{j}^{u}),\forall s\in\mathcal{S}\right\\},$ (16)

where
$B(z_{j}^{d},z_{j}^{u}):=\bigtimes\limits_{i=1}^{q_{j}}[z_{ji}^{d},z_{ji}^{u}]$.

The $k$-lower envelope in higher dimension is defined as follows:

######  Definition B.1.

Given a set of $q$-dimensional hyperplanes $\Pi$, the lower (resp. upper)
score of a point $t$ is the number of hyperplanes that lie strictly below
(resp. strictly above) $t$. The $k$-lower (resp. $k$-upper) envelope is the
closure of the set of points in $\Pi$ that have lower (resp. upper) score
equal to $k-1$. The $k$-lower envelope is also known as the $k$-level.

There is an extensive literature in the field of geometry dealing with the
development of $k$-envelopes, since they find applications in seemingly
unrelated problems like designing data structures (Chazelle and Preparata
(1986), Clarkson (1987), Agarwal and Matoušek (1995)), processing ranking
queries (Yu et al. (2012), Das et al. (2007), Shen et al. (2012), Cheema et
al. (2014)) or solving hyperplane partitioning problems (Borobia (1994), Tóth
(2000)). By geometric transformations, the time bounds found for $k$-envelopes
also apply to closely related structures (and vice versa) like _$k$ -order
Voronoi diagrams_ (Chan et al. (1997), Chazelle and Edelsbrunner (1987), Ramos
(1999)), or $k$-belts and $k$-sets:

######  Definition B.2.

The _$k$ -belt_ of a set of hyperplanes $\Pi$, for $0\leq
k\leq\left\lceil\frac{|\Pi|}{2}\right\rceil$, is the closure of the set of
points in the space that have lower and upper score greater than or equal to
$k$.

######  Definition B.3.

Given a fixed set $T$ of points in $\mathbb{R}^{d}$, a _$k$ -set_ ($k\leq|T|$)
of $T$ is any subset of $T$ of size $k$ which can be separated from the
complement of $T$ by a hyperplane. The $k$-set problem seeks for the maximum
possible number of $k$-sets.

Figure 3: In grey, the 4-belt of a set $\mathcal{L}_{j}$ of 8 lines in the
plane

Figure 3 shows (in light grey) the $4$-belt of a set of 8 lines. Obviously,
the boundary of the $k$-belt is the union of the $k$-upper and $k$-lower
envelopes. Less obvious is the fact that $k$-envelopes and $k$-sets are
closely related. In fact, using duality one can map a set of hyperplanes in
$\mathbb{R}^{d}$ (the primal space) into a set of points in $\mathbb{R}^{d}$
(the dual space). For instance, in two dimensions it suffices to map a point
$(\tilde{z},\tilde{t})$ to a line $L:\tilde{z}z-\tilde{t}$ and vice versa.
Hence, the worst-case complexity of the $k$-levels is within a constant factor
of the worst-case number of $k$-sets in point sets (Dey (1998), Yu et al.
(2012)).

In Edelsbrunner and Welzl (1986) it is proved that, if the intersections of
lines are computed dynamically in their manner, then Algorithm 2 constructs
the boundary of the $k$-belt of $\mathcal{L}$ in
$O(n_{k}\log^{2}|\mathcal{L}|)$ time and $O(n_{k}+|\mathcal{L}|)$ space, where
$n_{k}$ denotes the maximum number of edges bounding the $k$-belt of any set
$\mathcal{L}$ with $|\mathcal{L}|$ lines. Upper bounds on $n_{k}$ are given in
Edelsbrunner and Welzl (1985) by transforming the geometric problem into a
combinatorial one. In particular, the upper bound established depends on $k$
and $|\mathcal{L}|$: $n_{k}=O(|\mathcal{L}|k^{\frac{1}{2}})$. Dey (1998), Tóth
(2000), Nivasch (2008) have subsequently improved the upper bound of the
complexity on the $k$-set problem in the two-dimensional case. As for the
development of efficient algorithms to compute $k$-envelopes or $k$-sets in
higher dimensional spaces, they constitute a major challenge mainly because
the computational complexity of $k$-envelopes increases exponentially with the
increase in dimensionality (Agarwal and Matoušek (1995)). However, there exist
algorithms for the cases $q=3$ (Sharir et al. (2000), Ramos (1999), Tóth
(2000)), $q=4$ (Sharir (2011)) and $q\geq 5$ (Agarwal et al. (1998), Halperin
and Sharir (2017)).

As for the convexification procedure, the convex hull problem is a classical
problem with immense literature. The generalization of the Jarvis March to the
multidimensional case is known as the _graph transversal method_ (and also as
the _gift-wrapping_ algorithm), whereas the Graham Scan extension can be found
under the name of _incremental method_. For small dimensions ($q=2,3$), we
encounter parallel algorithms in the literature that use the _divide-and-
conquer_ paradigm to obtain better running times than the ones proposed here.
However, we chose sequential ones for the ease of computation and explanation,
and mainly because, given that the sets studied in our computational
experiments are quite small, the difference in performance in the
preprocessing steps is negligible. Surveys on algorithms and on the complexity
of the convex hull problem can be found, for instance, in Tóth et al. (2017),
Preparata and Shamos (2012).