# Learning Decision-Focused Uncertainty Sets
in Robust Optimization

Irina Wang, Cole Becker, Bart Van Parys, and Bartolomeo Stellato

###### Abstract

We propose a data-driven technique to automatically learn the uncertainty sets
in robust optimization. Our method reshapes the uncertainty sets by minimizing
the expected performance across a family of problems subject to guaranteeing
constraint satisfaction. Our approach is very flexible and can learn a wide
variety of uncertainty sets while preserving tractability. We solve the
constrained learning problem using a stochastic augmented Lagrangian method
that relies on differentiating the solutions of the robust optimization
problems with respect to the parameters of the uncertainty set. Due to the
nonsmooth and nonconvex nature of the augmented Lagrangian function, we apply
the nonsmooth conservative implicit function theorem to establish convergence
to a critical point, which is a feasible solution of the constrained problem
under mild assumptions. Using empirical process theory, we show finite-sample
probabilistic guarantees of constraint satisfaction for the resulting
solutions. Numerical experiments show that our method outperforms traditional
approaches in robust and distributionally robust optimization in terms of out-
of-sample performance and constraint satisfaction guarantees.

##  1 Introduction

Over the past years, robust optimization (RO) has become a widely adopted
efficient tool for decision-making under uncertainty. A robust approach first
defines an uncertainty set where the uncertain vector lives, and subsequently
hedges its decision against any worst-case realization in this set (Ben-Tal
and Nemirovski,, 2008; Ben-Tal et al.,, 2009; Bertsimas et al.,, 2011;
Bertsimas and den Hertog,, 2022). To attain proper solution quality and
practical robustness the choice of uncertainty set is paramount.

Traditionally, uncertainty sets are shaped and sized to guarantee they contain
an appropriately large fraction of the probability mass of the unknown
distribution generating the samples (i.e., high coverage). The main appeal of
this approach is that the uncertainty set is independent of the subsequent
optimization problem considered. In this vein, many approaches to designing
uncertainty sets assume structural information of the underlying
distributions, and rely on these a priori assumptions to build guarantees of
constraint satisfaction (Ben-Tal and Nemirovski,, 2000; Bertsimas et al.,,
2011; Bandi and Bertsimas,, 2012). In particular, a large amount of existing
literature studies how to bound the probability of constraint satisfaction
based on the size and shape of the uncertainty sets (Li et al.,, 2012; Guzman
et al.,, 2017; Bertsimas et al.,, 2021), where the aforementioned structural
information is exploited, but the optimization problem itself is not
leveraged. By decoupling the design of the uncertainty set from the subsequent
optimization problem, the resulting theoretical guarantees can often be quite
conservative (Roos and den Hertog,, 2020; Gao and Kleywegt,, 2023).

The trend of constructing uncertainty sets separately from the subsequent
optimization problem has continued in the current era of data-driven RO, where
the explosion in the availability of data has lead to new paradigms where
uncertainty sets are designed directly from data (Bertsimas et al.,, 2018).
Whether by combining a priori assumptions with the confidence regions of
statistical hypothesis tests on the data (Bertsimas et al.,, 2018), or by
approximating a high-probability-region using quantile estimation (Hong et
al.,, 2021), significant attention has been given to the probability
distribution generating the data and less so to the subsequent optimization
problem. More recently in distributionally robust optimization (DRO), the
constructed ambiguity set of distributions has historically been required to
contain the true data-generating distribution with a certain level of
confidence (Esfahani and Kuhn,, 2018). Again, the intuitive appeal of this
approach is that the design of the ambiguity set can be carried out in
isolation of the subsequent optimization problem.

Unfortunately, decoupling the uncertainty set construction from the subsequent
optimization problem in a two-step approach may lead to very conservative
decisions. In DRO with Wasserstein ambiguity sets, the two-step approach
suffers from the curse of dimensionality as their radius must scale with the
dimension of the uncertainty (Weed and Bach,, 2019). On the other hand, one
can overcome the curse of dimensionality by exploiting Lipschitz conditions on
the objective function of the subsequent optimization problem when designing
the ambiguity set (Gao,, 2023). The situation in DRO with $f$-divergence
ambiguity sets is even more striking. For instance, Kullback-Leibler ambiguity
sets around an empirical distribution can structurally not contain any
continuous distribution however large its radius is chosen. Consequently, if
the data-generating distribution is continuous a two-step approach is not
possible here. Nevertheless, the Kullback-Leibler can be shaped to ensure that
the solution has optimal statistical properties when one is willing to
exploiting the fact that the subsequent optimization problem has a continuous
loss function (Van Parys et al.,, 2021). In both cases, significant reduction
in conservatism can be achieved by designing the uncertainty set not only
based on the distributions generating the data, but also by exploiting
structural properties of the subsequent optimization problem.

The drive for even better performance in RO, then, leads us to a new
direction: a method of constructing uncertainty sets that deviates from the
conventional decoupled approach and explicitly takes into account the
optimization problem as well. This idea has been around for multiple years; as
stated by (Ben-Tal et al.,, 2009, page 34), the requirement of containing a
high probability mass is only a sufficient (and often restrictive and not
needed) condition to guarantee high probability of constraint satisfaction.
However, standardized methodology to exploit the notion have seldom been
developed. The customary restriction, in essence, allows the uncertainty set
developed to hold for all classes of optimization problems, which brings about
high levels of conservatism. However, the uncertainty set need only provide
high solution quality for the desired optimization problems, and in turn need
only include worst-case realizations of uncertainty for these instances.

In this work, we propose a new approach to calibrate uncertainty sets for a
particular family of parametric optimization instances. This is a common
scenario in many applications, including inventory management, where we solve
similar optimization problems with varying initial inventory levels (i.e., the
parameters), while satisfying uncertain demand with high probability. Our
method minimizes an expected objective over a family of parametric problems,
while ensuring data-driven probabilistic guarantees enforced as a conditional
value-at-risk constraint (Rockafellar and Uryasev,, 2002). We construct the
uncertainty set, which is the family-specific region of worst-case
realizations, by learning its shape and size automatically. Specifically, we
develop a stochastic augmented Lagrangian approach, relying on nonsmooth
implicit differentiation techniques (Agrawal et al., 2019b, ; Agrawal et al.,
2019a, ; Bolte et al.,, 2021; Bertrand et al.,, 2022) to compute the
derivative of the solution of robust optimization problem with respect to the
shape and size of the uncertainty sets. By incorporating the problem objective
and the data-driven constraint satisfaction into the learning problem, we
control the trade-off between performance and robustness. Our approach offers
significant performance improvements over traditional techniques constructing
uncertainty sets with high probability mass requirements, and is more flexible
than data-driven RO methods that only tune the radius of the uncertainty sets
(often via grid search), and not its shape.

###  1.1 Our contributions

In this work, we present a new approach to automatically learn the uncertainty
sets in RO. We obtain high-quality solutions across a family of parametric
optimization problems while ensuring probabilistic guarantees of constraint
satisfaction. By directly taking into account the quality of the RO solution,
the sets found are decision-focused and eliminate the unnecessary conservatism
common in traditional data-driven RO methods.

  * •

We minimize the expectation of the objective function over a family of
parametric RO problems, while ensuring an appropriate probability level of
constraint satisfaction. The latter is formulated as constraint on the
conditional-value-of-risk ($\mathop{\bf CVaR}$) of the problem family’s
uncertain constraints.

  * •

The aforementioned objective and constraints rely on the solutions of the
parametric RO problems. The overall formulation is thus a constrained bi-level
problem, with the inner level represented by the parametric RO problems.

  * •

We develop a stochastic augmented Lagrangian algorithm to solve the
constrained learning problem to find the optimal parameters of the uncertainty
set. Due to the potential nonsmooth, nonconvex, and therefore
nondifferentiable nature of the constructed problem, we show path-
differentiability of the Lagrangian, and apply nonsmooth implicit
differentiation. We show convergence to a critical point under step-size and
boundedness conditions, and establish the critical point as a feasible
solution of the constrained problem under mild assumptions.

  * •

Using empirical process theory, we give finite sample probabilistic guarantees
of constraint satisfaction for the resulting solution.

  * •

We benchmark our method on examples in portfolio optimization and multi-stage
inventory management, outperforming data-driven robust and distributionally
robust approaches in terms of out-of-sample objective value and probabilistic
guarantees. The code is available at
https://github.com/stellatogrp/lropt_experiments.

###  1.2 Related work

##### Data-driven robust optimization.

Data-driven robust optimization has gained wide popularity, by constructing
uncertainty sets through data-driven approximations of the unknown data-
generating distribution. Using hypothesis testing, Bertsimas et al., (2018)
pair a priori assumptions on the distribution with different statistical
tests, and obtain various uncertainty sets with different shapes,
computational properties, and modeling power. Similar approaches construct
data-driven sets using quantile estimation (Hong et al.,, 2021), type-$\infty$
Wasserstein distance (Bertsimas et al.,, 2022), and deep learning clustering
techniques (Goerigk and Kurtz,, 2023). All these methods rely on a two-step
procedure which separates the construction of the uncertainty set from the
resulting robust optimization problem, and often enforces a high coverage
requirement on the selected uncertainty sets. As noted by Ben-Tal et al.,
(2009), however, this requirement may be unnecessary. Recently, Gupta, (2019)
proves that while convex problems require DRO ambiguity sets with high
coverage to guarantee constraint satisfaction, traditional DRO ambiguity sets
for concave problems are significantly larger than the asymptotically optimal
sets. In order to find non coverage-based and potentially smaller uncertainty
sets that offer improved performance for targeted robust optimization
problems, we turn to ideas in differentiable optimization. Recent developments
in differentiable convex optimization have enabled a popular end-to-end
approach of training a predictive model to minimize loss on a downstream
optimization task (Donti et al.,, 2017; Demirović et al.,, 2019; Wilder et
al., 2019a, ; Wilder et al., 2019b, ; Cameron et al.,, 2022; Elmachtoub and
Grigas,, 2022), for which we also point to the survey by Kotary et al.,
(2021). Our approach adapts this idea, training instead the shape of the
uncertainty set on its performance on parametric robust optimization problems.
While developing this work, we became aware of a concurrent work with a
similar goal of fitting decision-focused uncertainty sets. In their work,
Schuurmans and Patrinos, (2023) construct DRO ambiguity sets which constrain
only the worst-case distribution along the direction in which the expected
cost of an approximate solution increases most rapidly. Our work, however,
considers a broader class of optimization problems and uncertainty sets, and
can directly tackle joint chance constraints using hyperparameter tuning.

##### Contexual optimization.

Contextual optimization and end-to-end approaches, beginning from the Smart
Predict, then Optimize framework (Elmachtoub and Grigas,, 2022) have been
shown to improve performance as opposed to traditional two-stage approaches,
where the model is trained separately from the optimization problem (Costa and
Iyengar,, 2023). Costa and Iyengar, (2023) propose an end-to-end
distributionally robust approach for portfolio construction, where they
integrate a _prediction layer_ that predicts asset returns with data on
financial features and historical returns, and use these predictions in a
_decision layer_ that solves the robust optimization problem. Concurrent to
our work, Chenreddy and Delage, (2024) propose an end-to-end training
algorithm to produce contextual uncertainty sets that achieve a certain
conditional coverage. This is to decrease the conservatism of their previous
two-stage approach, conditional robust optimization, where they learn a
contextually adapted uncertainty set known to contain with high probability
the realized uncertainty (Chenreddy et al.,, 2022). We, however, adopt a
different notion for context parameters, denoting them family parameters that
parameterize different instances of the downstream optimization problem,
instead of the uncertain parameter. Our learned uncertainty set therefore does
not learn a mapping from the context parameters to the uncertainty set (as in
contextual and conditional RO (Sadana et al.,, 2023)), as it may be
computationally expensive and hard to interpret. We do not enforce coverage
constraints, which may introduce unnecessary conservatism, and instead freely
shape the uncertainty sets based on their performance on the entire
parametralized family of robust problems. As pointed out in the recent survey
by Sadana et al., (2023), current work in the end-to-end framework seldom
focus on uncertainty in the downstream optimization problem, especially in the
constraints. Our approach thus fills this gap, where the process of
constructing the uncertainty set is linked in an end-to-end fashion to the
solution of the inner robust optimization problem. Our formulation is
flexible, as it allows for a wide variety of uncertain constraints and
uncertainty sets.

##### Modeling languages for robust optimization

Several open-source packages facilitate building models for RO, including as
ROME (Goh and Sim,, 2009) in Matlab, RSOME (Xiong and Chen,, 2021) in Python,
JuMPeR (Dunning et al.,, 2017) in Julia, and ROmodel (Wiebe and Misener,,
2022) for the Pyomo modeling language in Python. These packages mostly support
RO problems with affine uncertain constraints, which are solved via
reformulations and/or cutting plane procedures. Schiele et al., (2023)
recently introduced a Python extension of cvxpy (Diamond and Boyd,, 2016) for
solving the larger class of convex-concave saddle point problems, by
automating a conic dualization procedure reformulating min-max type problems
into min-min problems which can then be solved efficiently (Juditsky and
Nemirovski,, 2022). In the paper, they describe a set of rules termed
disciplined saddle programming (DSP) detailing when a saddle point
optimization problem can be equivalently formulated as convex optimization
ones, and efficiently be solved. In this work, we build a Python extension to
cvxpy tailored to formulating and solving RO problems. While some RO problems
can be described in terms of the DSP framework, our library augments DSP
capabilities in RO by supporting max-of-concave uncertain constraint and
structure to intuitively describe well known uncertainty sets for ease of use.
Finally, our method automatically learns the uncertainty set through
differentiable optimization via a tight integration with the works by Agrawal
et al., 2019a ; Agrawal et al., 2019b .

##### Differentiable convex optimization.

There has been extensive work on embedding differentiable optimization
problems as layers within deep learning architectures (Amos and Kolter,,
2017). Agrawal et al., 2019a  propose an approach to differentiating through
disciplined convex programs, and implement their methodology in cvxpy.
Specifically, they implement differentiable layers for disciplined convex
programs in PyTorch and TensorFlow 2.0, in the package cvxpylayers. They
convert the convex programs into conic programs, then implicitly differentiate
the residual map of the homogeneous self-dual embedding associated with the
conic program. These developments have fueled the recent advances in end-to-
end approaches, and is also the tool of choice in this work.

##### Nonsmooth implicit differentiation.

Traditional techniques in differentiable convex optimization rely on the
implicit function theorem, which requires derivatives to exist almost
everywhere. Recently, a nonsmooth variant has been established, using
conservative Jacobians, provided that a nonsmooth form of the classical
invertibility condition is fulfilled (Bolte et al.,, 2021). As a key
ingredient, definable and locally Lipschitz continuous functions are shown to
admit conservative Jacobians, which follow the chain rule, and are thus
compatible with backpropagation techniques in stochastic gradient descent
(Davis et al.,, 2018; Bolte et al.,, 2021). These techniques can be used to
obtain convergence for nonsmooth problems.

##### Automatic hyperparameter selection.

Implicit differentiation techniques can greatly accelerate tuning the
regularization parameters in regression and classification problems (Bertrand
et al.,, 2020, 2022; Blanchet et al.,, 2019). However, most literature focuses
on tuning hyperparameters of the objective function in regression and
classification tasks, such as hyper-parameters for Lasso regularization
(Bertrand et al.,, 2022). In this work, we tune various hyperparameters of the
uncertainty sets in RO using implicit differentiation techniques.

##### Risk measures.

Value at risk ($\mathbf{VaR}$) and conditional value at risk ($\mathop{\bf
CVaR}$) are common risk measures for quantifying costs encountered in the
tails of distributions. A $\mathbf{VaR}_{\eta}$ constraint is a chance
constraint that measures the maximum cost incurred with respect to a given
confidence level $\eta$, but while intuitive, the function is often
intractable, and does not quantify costs beyond the found threshold value.
Conditional value at risk, on the other hand, measures the expected cost given
that the cost exceeds $\mathbf{VaR}$, and is known to have a tractable
formulation (Uryasev and Rockafellar,, 2001; Rockafellar and Uryasev,, 2002).
Since its development, $\mathop{\bf CVaR}$ optimization has been a common
method for controlling risk in portfolio optimization and option hedging
problems, and has also been explored in risk-constrained Markov decision
problems, where the risk is presented as a constraint on the $\mathop{\bf
CVaR}$ of the cumulative cost (Chow et al.,, 2015). The use of $\mathop{\bf
CVaR}$ constraints to substitute chance constraints has also been explored in
reliability engineering, where the $\mathop{\bf CVaR}$ constraint has been
given a different name: buffered failure probability (Rockafellar and Royset,,
2010). A related concept for controlling risk is the buffered probability of
exceedance (bPOE) (Mafusalov and Uryasev,, 2018), which is equivalent to 1
minus the inverse of $\mathop{\bf CVaR}$, and has a level set equivalence with
$\mathop{\bf CVaR}$. In our work, we constrain risk, which in RO is the
probability of constraint violation, by constraining the $\mathop{\bf CVaR}$
of the uncertain constraint functions.

###  1.3 Layout of the paper

In Section 2, we introduce the notion of the reshaping parameters with a
motivating example. In Section 3, we describe the data-driven problem and the
algorithm for learning reshaping parameters, and prove its convergence. In
Section 4, we discuss the probabilistic guarantees implied by the
aforementioned formulations. In Section 5, we present various numerical
examples problems. For completeness, we give example convex reformulations of
robust optimization problems for common uncertainty sets, with our reshaped
parameters, in the appendices.

##  2 The robust optimization problem

We consider a family of parametric robust optimization (RO) problems, indexed
by the family parameter $y\in{\mbox{\bf R}}^{p}$, where each instance takes
the form:

$x(\theta,y)\in\begin{array}[t]{ll}\mathop{\rm argmin}&f(x)\\\ \mbox{subject
to}&g(x,u,y)\leq 0\quad\forall u\in{\mathcal{U}}(\theta),\end{array}$ (1)

where $x\in{\mbox{\bf R}}^{n}$ is the optimization variable, $u\in{\mbox{\bf
R}}^{m}$ is the uncertain parameter, $f:{\mbox{\bf R}}^{n}\to{\mbox{\bf R}}$
is the convex objective function, and $g:{\mbox{\bf R}}^{n}\times{\mbox{\bf
R}}^{m}\times{\mbox{\bf R}}^{p}\to{\mbox{\bf R}}$ is the uncertain constraint
which we assume to be the maximum of functions concave in $u$ and convex in
$x$ and $y$, e.g., $g(x,u,y)=\max_{l=1,\dots,L}g_{l}(x,u,y)$. This structure
allows us to express $L$ joint uncertain constraints using the maximum of each
component $l$. Problem (1) guarantees that, for a given $y$, the solution
$x(\theta,y)$ satisfies the constraints for any value of the uncertain
parameter $u$ in ${\mathcal{U}}(\theta)\subseteq{\mbox{\bf R}}^{m}$, a convex
uncertainty set that depends on the set parameter $\theta\in{\mbox{\bf
R}}^{q}$. For example, $\theta$ can represent the center, size, and shape of
${\mathcal{U}}(\theta)$. In this setting, the family parameter $y$ represents
information available when we make our decisions (e.g., current allocations in
portfolio optimization). The uncertain parameter $u$, however, is not known
(e.g., future returns in portfolio optimization). Importantly, we do not
assume that $y$ affects $u$, which is common in several practical applications
(e.g., when our portfolio allocations do no affect future market returns).

Here, the set ${\mathcal{U}}(\theta)$ does not depend on individual $y$ as
this would require either learning a different uncertainty set for every
instance of the family, which may be costly, or learning a mapping from $y$ to
the uncertainty set (as in conditional RO (Sadana et al.,, 2023)), which may
be hard to interpret. Instead, we construct a fixed set that works for all
realizations of $y$ encountered in target application.

Our goal is to find a set parameter $\theta$ in a decision-focused manner,
i.e., such that the solutions $x(\theta,y)$ perform well in terms of objective
value and constraint satisfaction guarantees across the family of problems.
Specifically, we choose the uncertainty set ${\mathcal{U}}(\theta)$ such that
the solution $x(\theta,y)$ from (1) implies that the uncertain constraint is
satisfied with probability $1-\eta$ with $\eta>0$, across all instances of
problem (1):

${\mathbf{P}}_{(u,y)}(g(x(\theta,y),u,y)\leq 0)\geq 1-\eta.$ (2)

Here, ${\mathbf{P}}_{y}$ is the distribution of $y$ with domain
$\mathcal{D}_{y}\subseteq{\mbox{\bf R}}^{p}$, ${\mathbf{P}}_{u}$ is the
distribution of $u$ with domain $\mathcal{D}_{u}\subseteq{\mbox{\bf R}}^{m}$,
and ${\mathbf{P}}_{(u,y)}$ is their joint distribution, all of which are
unknown. We note that this is an aggregate probabilistic guarantee for the
parametric family of problems, and not conditional guarantee which holds for
all instances $y$. Whereas such conditional guarantee (Chenreddy and Delage,,
2024) are perhaps more desirable they require stringent assumptions on the
distribution ${\mathbf{P}}_{y}$ which are typically hard to verify in
practice. We can formulate constraint (2) in terms of the corresponding value
at risk, i.e.,

$\mathop{\bf
VaR}(g(x(\theta,y),u,y),\eta)=\inf\\{\gamma\mid{\mathbf{P}}_{(u,y)}(g(x(\theta,y),u,y)\leq\gamma)\geq
1-\eta\\}\leq 0.$

Unfortunately, the value at risk function often leads to intractable
formulations (Uryasev and Rockafellar,, 2001). To get a tractable
approximation, we can adopt the conditional value at risk (Uryasev and
Rockafellar,, 2001; Rockafellar and Uryasev,, 2002), defined as

$\mathop{\bf
CVaR}(g(x(\theta,y),u,y),\eta)=\inf_{\alpha}\left\\{\frac{{\mbox{\bf
E}}_{(u,y)}(g(x(\theta,y),u,y)-\alpha)_{+}}{\eta}+\alpha\right\\},$ (3)

where $(a)_{+}=\max\\{a,0\\}$. It is well known (Uryasev and Rockafellar,,
2001) that, for any $x\in{\mbox{\bf R}}^{n}$, the relationship between these
probabilistic guarantees of constraint satisfaction is

$\mathop{\bf CVaR}(g(x,u,y),\eta)\leq 0~{}~{}\Longrightarrow~{}~{}\mathop{\bf
VaR}(g(x,u,y),\eta)\leq
0~{}~{}\Longleftrightarrow~{}~{}{\mathbf{P}}_{(u,y)}(g(x,u,y)\leq 0)\geq
1-\eta.$ (4)

Therefore, if the solution $x(\theta,y)$ of (1) satisfies $\mathop{\bf
CVaR}(g(x(\theta,y),u,y),\eta)\leq 0$, we have the desired probabilistic
guarantee (2). This constraint, together with the goal of minimizing the
expected objective value, are constructed into an evaluation metric for
updating $\theta$; see Figure 1. Typically, by constructing uncertainty sets
that contain at least $1-\eta$ probability mass, i.e.,
${\mathbf{P}}_{u}(u\in{\mathcal{U}}(\theta))\geq 1-\eta$, we can ensure that
any feasible solution of (1) implies a probabilistic guarantee (2). However,
this sufficient condition may lead to overly conservative solutions (Ben-Tal
et al.,, 2009, page 34), as it works for any RO problem and any constraint
using the uncertainty set ${\mathcal{U}}(\theta)$. In this work, instead, we
avoid such conservatism by taking into account the cost of the robust
solutions $x(\theta,y)$ while constructing the uncertainty sets.

Figure 1: Our high-level procedure. We evaluate the performance of the
decision $x(\theta,y)$ in terms of the expected objective value and constraint
satisfaction. The evalution metric is defined in Section 3.2.

###  2.1 Motivating example

We consider the portfolio management problem where we select a portfolio
$x\in{\mbox{\bf R}}^{n}$ of stocks to maximize $u^{T}x$, where $u\in{\mbox{\bf
R}}^{n}$ are the uncertain asset returns. For visualization purposes, we
choose $n=2$. The trading costs are proportional to the deviations from our
previous holdings $x^{\text{prev}}\in{\mbox{\bf R}}^{n}$, i.e., $\lambda^{\rm
trade}\|x-x^{\text{prev}}\|_{1}$ with $\lambda^{\rm trade}>0$. We seek to
minimize the objective

$-u^{T}x+\lambda^{\rm trade}\|x-x^{\text{prev}}\|_{1}.$

To account for uncertainty in the returns $u$, we introduce a new variable $t$
to write the objective in epigraph form, obtaining the RO problem

$\begin{array}[]{ll}\mbox{minimize}&t+\lambda^{\rm
trade}\|x-x^{\text{prev}}\|_{1}\\\ \mbox{subject to}&-u^{T}x\leq t\quad\forall
u\in{\mathcal{U}}(\theta)\\\ &\mathbf{1}^{T}x=1,\quad x\geq 0.\end{array}$ (5)

We would like to construct uncertainty sets without knowing the distributions
of $y$ and $u$.

##### Problem setup.

The uncertain parameter $u$ is distributed as a normal distribution, with
parameters

$\mu=\left[\begin{array}[]{l}0.3\\\
0.3\end{array}\right],\quad\Sigma=\left[\begin{array}[]{ll}\phantom{-}0.5&-0.3\\\
-0.3&\phantom{-}0.4\end{array}\right].$

We have access to $500$ samples. We parametrize the problem (5) by
$y=x^{\text{prev}}$, the previous holdings. We generate 10 values of $y$ on
the standard simplex using a Dirichlet distribution, with parameter value
$(2.5,1)$, corresponding to the 2 different stocks. We set $\lambda^{\rm
trade}=0.2$. We generate data for $u$ and $y$ as detailed above, for 5
repetitions of the experiment. Below, we show the exact values of the
uncertainty sets obtained for one such repetition.

Figure 2: Top (bottom) row: mean-variance (reshaped) uncertainty set
corresponding to the empirical probability of constraint violation
$\hat{\eta}$. White dots: training data points. Green dot: training data mean.
Red ellipsoid: mean-variance uncertainty set. Purple lines: lines
corresponding to a constraint function $-u^{T}x(\theta,y)-t(\theta,y)=0$; for
each plot, there is a line for each data driven solution $x(\theta,y)$ and
$t(\theta,y)$. Middle row: out-of-sample trade-off curves for the objective
value vs. empirical probability of constraint violation, averaged across all
problem instances and all repetitions of the experiment. The shaded regions
represent the 0.3 to 0.7 quantiles of the repetitions. Each marker corresponds
to an $\rho$ value. The uncertainty sets depicted above and below achieves the
probabilities of constraint violation $\hat{\eta}$ given by the red dotted
lines, for one repetition of the experiment.

##### Mean-Variance uncertainty set.

Standard methods in RO methods construct uncertainty sets based on the
empirical mean $\hat{\mu}$ and covariance $\hat{\Sigma}$ of the uncertainty,

${\mathcal{U}}{(\theta)}=\\{u=\hat{\mu}+\hat{\Sigma}^{1/2}z\mid\|z\|_{2}\leq\rho\\}=\\{b^{\rm
mv}+A^{\rm mv}z\mid\|z\|_{2}\leq\rho\\},$ (6)

where $\theta=(A^{\rm mv},b^{\rm mv},\rho)$. The parameter $\rho$ represents
the size of the uncertainty set, and, for one set of data, we have

$A^{\rm
mv}=\hat{\Sigma}^{1/2}=\left[\begin{array}[]{ll}\phantom{-}0.69&-0.24\\\
-0.24&\phantom{-}0.59\end{array}\right],\quad b^{\rm
mv}=\hat{\mu}=\left[\begin{array}[]{l}0.28\\\ 0.33\end{array}\right].$

##### Decision-focused (reshaped) uncertainty set.

Consider now the same set as in (6) with different shape, i.e.,
${\mathcal{U}}(\theta)=\\{u=b^{\rm re}+A^{\rm re}z\mid\|z\|_{2}\leq\rho\\}$,
where $\theta=(A^{\rm re},b^{\rm re},\rho)$ and

$A^{\rm re}=\left[\begin{array}[]{ll}0.42&\phantom{-}0.83\\\
0.78&-0.08\end{array}\right],\quad b^{\rm re}=\left[\begin{array}[]{l}0.30\\\
0.41\end{array}\right].$

##### Comparison.

In Figure 2, we compare the two uncertainty sets where we vary $\rho$ while
$A^{\rm st},b^{\rm st}$ and $A^{\rm re},b^{\rm re}$ are fixed. For each value
of $\rho$ we compare the average out-of-sample objective value at the
optimizers $x(\theta,y)$ against the out-of-sample empirical probability of
constraint violation. The average is taken across problem instances $y$’s,
testing datasets of $N=4500$ uncertainty realizations, and over all
repetitions of the experiment. We note that the reshaped uncertainty set
achieves data-driven solutions that give, on average, better trade-offs
between the out-of-sample objective value and empirical probability of
constraint violation. While the mean-variance uncertainty set conforms to the
shape of data, it might be conservative because it does not take into account
the structure of the optimization problem. In contrast, the reshaped
uncertainty set gives better (lower) out-of-sample costs for the same
probability of constraint satisfaction, even though it does not conform as
well to the shape of the data, and covers a lower percentage of the data
points. For some target values of empirical probability of constraint
violation, $\hat{\eta}$, we plot the constraint values evaluated at
$x(\theta,y)$ and $t(\theta,y)$ for all $y$, i.e. the contour lines
$-u^{T}x(\theta,y)-t(\theta,y)=0$, for both the mean-variance and reshaped
sets. We note that across the family of $y$’s, the reshaped uncertainty set
identifies reasonable worst-case realizations of the uncertainty, denoted by
the intersections between the lines and the uncertainty set, despite the set’s
less conventional orientation.

In the next sections, we introduce an automatic technique find such sets with
better performance, by tuning their shape and size such that the problem
objective is minimized while guaranteeing constraint satisfaction across the
family of parametric RO problems (1).

##  3 Learning the uncertainty set

###  3.1 Constrained learning formulation

Consider a RO problem (1). We have a dataset of $N$ independent samples of the
uncertain parameter $u$,
$U^{N}=\\{u^{i}\\}_{i=1}^{N}\subseteq\mathcal{D}_{u}$, governed by the product
distribution ${\mathbf{P}}^{N}_{u}$, and a dataset of $J$ independent samples
of the family parameter $y$,
$Y^{J}=\\{y^{j}\\}_{j=1}^{J}\subseteq\mathcal{D}_{y}$, governed by the product
distribution ${\mathbf{P}}_{y}^{J}$. We define the combined samples
$(u^{i},y^{j})$ for $i=1,\dots,N$ and $j=1,\dots,J$, and the combined dataset
$U^{N}\times Y^{J},$ with product distribution defined as
${\mathbf{P}}^{N\times J}$. Our main requirement is to identify a set
parameter $\theta$ for which the solutions $x(\theta,y)$ to problem (1) with
uncertainty set $\mathcal{U}(\theta)$, satisfy the following finite-sample
probabilistic guarantee

${\mathbf{P}}^{N\times J}\left\\{{\mathbf{P}}_{(u,y)}(g(x(\theta,y),u,y)\leq
0)\geq 1-\eta\right\\}\geq 1-\beta,$ (7)

where $\beta>0$ is the target confidence level across draws of the dataset
$U^{N}\times Y^{J}$. Instead of working with the probability
${\mathbf{P}}_{(u,y)}$, which often leads to intractable formulations (Uryasev
and Rockafellar,, 2001) as discussed in Section 2, we adopt a tractable
approximation based on the conditional value at risk (Uryasev and
Rockafellar,, 2001; Rockafellar and Uryasev,, 2002) defined in (3), obtaining
the constraint

$\mathop{\bf
CVaR}(g(x(\theta,y),u,y),\eta)\leq\kappa\quad\iff\quad\inf_{\alpha}\left\\{\frac{{\mbox{\bf
E}}_{(u,y)}(g(x(\theta,y),u,y)-\alpha)_{+}}{\eta}+\alpha\right\\}\leq\kappa,$

where $\kappa>0$ is a threshold value that we define in Section 4 to obtain
the finite sample guarantees (7). Our goal is to find the set parameter
$\theta$ that minimizes the average objective of the RO problem across
realizations of the family parameter $y$, which can be written as

$\begin{array}[]{ll}\mbox{minimize}&{\mbox{\bf
E}}_{y}\left(f(x(\theta,y))\right)\\\ \mbox{subject to}&{\mbox{\bf
E}}_{(u,y)}\left((g(x(\theta,y),u,y)-\alpha)_{+}/\eta+\alpha\right)\leq\kappa,\\\
\end{array}$ (8)

where $x(\theta,y)$ are solutions of the inner RO problems, defined as in (1).
We given in Appendix A examples of parameterized uncertainty sets and the
corresponding convex reformulations of (1), which we use to solve for
$x(\theta,y)$. The expectation constraint corresponds to the $\mathop{\bf
CVaR}$ constraint where $\alpha$ is shifted to become a minimization variable.
When the functions $f$ and $g$ are convex in $x$, this shift in $\alpha$ has
no impact on the optimal solution (Sarykalin et al.,, 2008). With this
formulation, we tune $\alpha$ and the set parameter $\theta$ to improve the
quality of the solutions $x(\theta,y)$ of the inner RO problem (1). Solving
the training problem (8) is challenging because of the stochasticity in the
objective and the constraints, and the nonconvex bi-level structure introduced
by the inner RO problem. We address these challenges using a stochastic
augmented Lagrangian algorithm.

###  3.2 Training the uncertainty set using a Stochastic augmented Lagrangian
algorithm

As the distributions ${\mathbf{P}}_{u}$ and ${\mathbf{P}}_{y}$ are unknown, we
are unable to obtain the true expectation in (8). As such, we work with
sampled versions of the expectations. We define a condition on the sampled
$\mathop{\bf CVaR}$,

$\widehat{\mathop{\bf
CVaR}}(\theta)=\inf_{\alpha}\frac{1}{NJ}\sum_{j=1}^{J}\sum_{i=1}^{N}\left(\frac{(g(x(\theta,y^{j}),u^{i},y^{j})-\alpha)_{+}}{\eta}+\alpha\right)\leq\kappa,$
(9)

and require this condition to hold for feasible $\theta$. To find such a
$\theta$, we introduce finite-sample approximations of the objective and
constraint functions,

$F(\theta)=\frac{1}{J}\sum_{j=1}^{J}f(x(\theta,y^{j}))\quad\text{and}\quad
H(\alpha,\theta)=\frac{1}{NJ}\sum_{j=1}^{J}\sum_{i=1}^{N}(g(x(\theta,y^{j}),u^{i},y^{j})-\alpha)_{+}/\eta+\alpha,$

evaluated over the dataset $U^{N}\times Y^{J}$. We use these approximations to
define the sampled version of problem (8),

$\begin{array}[]{ll}\mbox{minimize}&F(\theta)\\\ \mbox{subject
to}&H(\alpha,\theta)+s=\kappa\\\ &s\geq 0,\\\ \end{array}$ (10)

where $s$ is the non-negative slack variable in ${\mbox{\bf R}}^{d}$ to
represent the inequality constraint as an equality. We can now write the
augmented Lagrangian function as

$\displaystyle L(\alpha,s,\theta,\lambda,\mu)=F(\theta)$
$\displaystyle+\lambda(H(\alpha,\theta)+s-\kappa)+\frac{\mu}{2}\|H(\alpha,\theta)+s-\kappa\|^{2},$
(11)

where $\lambda$ is the Lagrange multiplier in R, and $\mu$ is the penalty
parameter in R. This function $L$ is evaluated over the entire dataset
$U^{N}\times Y^{J}$. We denote the estimated augmented Lagrangian over a
subset (batch) of data $B=(U\times Y)\subseteq(U^{N}\times Y^{J})$ as
$\hat{L}$. Using the stochastic augmented Lagrangian approach outlined in
Algorithms 1 and 2, and under mild assumptions, we can find a feasible
$\theta$ satisfying condition (9). We also define a slightly weaker optimality
condition: an _$\epsilon$ -stationary point_ for problem (10). Given
$\epsilon>0$, a point $(\alpha,s,\theta)\in{\mbox{\bf R}}^{q+2}$ is called an
$\epsilon$-stationary point of (10) if there is a scalar $\lambda\in{\mbox{\bf
R}}$ such that

$\|H(\alpha,\theta)+s-\kappa\|_{2}\leq\epsilon,\quad\mathop{\bf
dist{}}(0,J_{F}(\theta)+\lambda
J_{H+s-\kappa}(\alpha,s,\theta)))\leq\epsilon,$ (12)

where $J_{F}(\theta)$ is the conservative Jacobian of the objective function
$F$ in the variable $\theta$, $J_{H+s-\kappa}(\alpha,s,\theta)$ denotes the
conservative Jacobian (Definition 3.1 below) of the constraint function in
variables $(\alpha,s,\theta)$, and the distance between a vector $x$ and set
$\mathcal{X}$ is defined as $\mathop{\bf
dist{}}{(x,\mathcal{X})}=\min_{y\in\mathcal{X}}\|x-y\|_{2}$. In addition,
given $\lambda$ and $\mu$, a point $(\alpha,s,\theta)\in{\mbox{\bf R}}^{q+2}$
is called a stationary point of the augmented Lagrangian function (11), if
$0\in J_{L}(\alpha,s,\theta)$. While useful in our convergence analysis, this
stationarity condition is hard to evaluate. That’s why in practical
applications, instead of using it to terminate Algorithm 2, we set the number
of iterations $t_{\text{max}}$ to a reasonably large value.

Algorithm 1 Safeguarded stochastic augmented Lagrangian algorithm to solve
(10)

1:given $\theta^{0}=(A^{\text{init}},b^{\text{init}})$, $\alpha^{0}$, $s^{0}$,
$\lambda^{0},\mu^{0},\epsilon\geq 0,\kappa$, $\sigma>1$, $\tau<1$,
$H_{\text{best}}=\infty$

2:for $k=1,\dots,k_{\rm max}$ do

3: $(\alpha^{k},s^{k},\theta^{k})\leftarrow$ minimize
$L(\alpha,s,\theta,\lambda^{k-1},\mu^{k-1})$ over $(\alpha,s,\theta)$ via
Algorithm 2

4: if $\|H(\alpha^{k},\theta^{k})+s^{k}-\kappa\|_{2}\leq\epsilon$ then

5: return $\theta^{k}$ $\triangleright$ $\epsilon$-stationary point found

6: else if
$\|H(\alpha^{k},\theta^{k})+s^{k}-\kappa\|_{2}\leq\tau\|H_{\text{best}}\|_{2}$
then

7:
$\lambda^{k}\leftarrow\Pi_{[{\lambda^{\rm{min}}},\lambda^{\rm{max}}]}(\lambda^{k-1}+\mu^{k-1}(H(\alpha^{k},\theta^{k})+s^{k}-\kappa)$)
$\triangleright$ update Lagrange multipliers

8: $H_{\text{best}}\leftarrow\|H(\alpha^{k},\theta^{k})+s^{k}-\kappa\|_{2}$

9: $\mu^{k}\leftarrow\mu^{k-1}$

10: else

11: $\mu^{k}\leftarrow\sigma\mu^{k-1}$ $\triangleright$ update penalty

12: $\lambda^{k}\leftarrow\lambda^{k-1}$

13:return $\theta^{k_{\rm max}}$

Algorithm 2 Stochastic gradient descent for the inner loop

1:given $\theta^{0}$, $\alpha^{0}$, $s^{0}$,
$\lambda,\mu,\\{\delta^{t}\\}_{t\in\mathbf{N}},U^{N}\times Y^{J}$

2:for $t=1,\dots,t_{\rm max}$ do

3: sample data $(U^{t}\times Y^{t})\subseteq(U^{N}\times Y^{J})$

4: $x^{t}(\theta^{t-1},y)\leftarrow$ solve inner RO problem (1) for all $y\in
Y^{t}$

5: $J_{\hat{L}}\leftarrow J_{\hat{L}}(\alpha^{t-1},\theta^{t-1},s^{t-1})$
$\triangleright$ compute Jacobian of $\hat{L}$ over $U^{t}\times Y^{t}$

6:
$\alpha^{t}\leftarrow\alpha^{t-1}-\delta^{t}v_{\alpha}^{t},\quad\text{with}\quad
v_{\alpha}^{t}\in J_{\hat{L}}(\alpha)$

7: $s^{t}\leftarrow(s^{t-1}-\delta^{t}v_{s}^{t})_{+},\quad\text{with}\quad
v_{s}^{t}\in J_{\hat{L}}(s)$

8:
$\theta^{t}\leftarrow\theta^{t-1}-\delta^{t}v_{\theta}^{t},\quad\text{with}\quad
v_{\theta}^{t}\in J_{\hat{L}}(\theta)$

9:return $(\alpha^{t},s^{t},\theta^{t})$

####  3.2.1 Path differentiability and conservative Jacobians

In order to apply Algorithms 1 and 2, we need to establish the existence of
$J_{L}$ and $J_{\hat{L}}$, the conservative Jacobians of the augmented
Lagrangian functions $L$ and $\hat{L}$. Conservative Jacobians, introduced in
(Bolte and Pauwels,, 2019), are a generalized form of Jacobians for functions
that are almost everywhere differentiable.

######  Definition 3.1 (Conservative Jacobian (Bolte et al.,, 2021)).

Given a locally Lipschitz continuous function $F:{\mbox{\bf
R}}^{m}\rightarrow{\mbox{\bf R}}^{m}$, we call $J_{F}:{\mbox{\bf
R}}^{n}\rightrightarrows{\mbox{\bf R}}^{n\times m}$ its conservative Jacobian
if $J_{F}$ has a closed graph, is locally bounded, and is nonempty with

$\frac{\text{d}}{\text{d}t}F(\gamma(t))\in
J_{F}(\gamma(t))\dot{\gamma}(t)~{}~{}\text{a.e.},$

for any absolutely continuous curve $\gamma$ in ${\mbox{\bf R}}^{n}$. A
locally Lipschitz function $F$ that admits a conservative Jacobian $J_{F}$ is
called path differentiable.

To obtain the conservative Jacobians $J_{L}$ and $J_{\hat{L}}$, we note the
following assumptions.

######  Assumption 3.1 (Definability).

All functions $f$ and $g_{l}$ in the robust optimization problem (1) are
_definable_ (van den Dries and Miller,, 1999; Coste,, 2000; Bolte et al.,,
2021).

The class of definable functions is very wide and includes semi-algebraic
functions; for an overview, see the works by Bolte and Pauwels, (2019) and
Bolte et al., (2021). For a foundational work, see van den Dries and Miller,
(1999). Assuming the definability of the objective and constraint functions,
it follows that $F$ and $H$, and subsequently $L$, as compositions of
definable pieces, are also definable (Bolte et al.,, 2021).

######  Assumption 3.2 (Local Lipschitzness).

The Lagrangian function $L(\alpha,s,\theta,\lambda,\mu)$ is locally Lipschitz
in variables $\alpha,s,\theta,\lambda,\mu$. A function $F:{\mbox{\bf
R}}^{n}\rightarrow{\mbox{\bf R}}^{m}$ is locally Lipschitz if, for each
$x\in{\mbox{\bf R}}^{n}$, there exists a neighborhood $\mathcal{V}$ of $x$
such that $F$ is Lipschitz on $\mathcal{V}$.

We also require the augmented Lagrangian to be locally Lipschitz continuous,
as definable locally Lipschitz functions are shown to be path differentiable,
meaning they admit conservative Jacobians and follow the chain rule (Bolte et
al.,, 2021). Admittance of the chain rule is imperative, as $L$ depends on
$\theta$ indirectly through the solution $x(\theta,y)$ of the inner robust
problem (1). In order to establish the conservative Jacobians of $L$ in the
variable $\theta$, we thus need to show the path differentiability of the
inner problem and compute $J_{x}(\theta)$, for which we require a uniqueness
assumption.

######  Assumption 3.3 (Uniqueness).

The robust optimization problem (1) has a convex conic reformulation with a
unique solution.

We can then make use of the results of Agrawal et al., 2019b  and Bolte et
al., (2021), which enables us to find conservative Jacobians of optimal
solutions for conic problems by differentiating through the KKT optimality
conditions. We state the theorem below, with the proof delayed to Appendix C.

######  Theorem 3.1 (Path differentiability of the RO solution).

Suppose Assumption 3.3 holds. Let $x(\theta,y)$ be the unique solution, for a
given input pair $(\theta,y)$, to the convex conic reformulation of the inner
optimization problem (1). This solution is path differentiable with nonempty
conservative Jacobian $J_{x}(\theta)\neq\emptyset$.

We can then evaluate the conservative Jacobians of the augmented Lagrangian
function.

######  Proposition 3.1 (Conservative Jacobians of the augmented Lagrangian
(11)).

Suppose Assumptions 3.1 and 3.2 hold, and for each given input pair
$(\theta,y)$, the solution $x(\theta,y)$ of the optimization problem (1) is
path differentiable with conservative Jacobian $J_{x}(\theta)$. Then, the
Lagrangian function $L(\alpha,s,\theta,\lambda,\mu)$, from equation (11), and
its batched version $\hat{L}(\alpha,s,\theta,\lambda,\mu)$, are path
differentiable with conservative Jacobians $J_{L}(\alpha,s,\theta)$ and
$J_{\hat{L}}(\alpha,s,\theta)$.

The specific forms of the conservative Jacobians established in Proposition
3.1 are given in Appendix B.

####  3.2.2 Convergence analysis

To show convergence of the augmented Lagrangian algorithm, given in Algorithm
1, we require the convergence of the inner loop, Algorithm 2. The convergence
of Algorithm 2 relies on the path-differentiability of the Lagrangian
function, as well as the boundedness of the iterates and the Lagrangian
multiplier $\lambda$.

######  Assumption 3.4 (Boundedness).

There exists a constant $\alpha^{\rm max}>0$ and compact set
$\mathcal{D}_{\alpha}\subset{\mbox{\bf R}}$, such that, for
$\alpha^{0}\in\mathcal{D}_{\alpha}$, $\|\alpha^{t}\|\leq\alpha^{\rm max}$
almost surely. Similarly, there exist constants $\theta^{\rm max},s^{\rm
max},x^{\rm max},y^{\rm max}>0$, compact sets
$\mathcal{D}_{\theta}\subset{\mbox{\bf R}}^{q}$,
$\mathcal{D}_{s}\subset{\mbox{\bf R}}$, and $\mathcal{X}\subset{\mbox{\bf
R}}^{n}$, such that, for $\theta^{0}\in\mathcal{D}_{\theta}$,
$s^{0}\in\mathcal{D}_{s}$, and $y\in\mathcal{D}_{y}$, we have
$\|\theta^{t}\|\leq\theta^{\rm max}$, $\|s^{t}\|\leq s^{\rm max}$, and
$x(\theta^{t},y)\in\mathcal{X}=\\{x\in{\mbox{\bf R}}^{n}\mid\|x\|\leq x^{\rm
max}\\}$ almost surely. In addition, the optimal Lagrangian multiplier
$\lambda^{\star}\in[\lambda^{\rm{min}}$, $\lambda^{\rm{max}}]$.

In addition, we require that the step sizes are nonnegative, square summable,
but not summable.

######  Assumption 3.5 (Step sizes).

The step sizes $\\{\delta^{t}\\}_{t\in\mathbf{N}}$ are chosen such that

$\displaystyle\sum_{t=1}^{\infty}\delta^{t}=\infty,\quad\sum_{t=1}^{\infty}(\delta^{t})^{2}<\infty,\quad\delta^{t}=o(1/\log(t)).$

Under these assumptions, by Theorem 6.2 of Davis et al., (2018) and Theorem 3
of Bolte et al., (2021), the sequence $\\{\alpha^{t},s^{t},\theta^{t}\\}$ in
Algorithm 2 converges to a stationary point of the Lagrangian function (11)
with arbitrary $\lambda$ and $\mu$.

We note that $\lambda^{k}$ and $\mu^{k}$ in Algorithm 1 are potentially
unbounded, in which case they do not converge. The stationary point
$(\alpha^{k},s^{k},\theta^{k})$ of the augmented Lagrangian function for any
$\lambda^{k}$ and $\mu^{k}$ is, therefore, not guaranteed to be a feasible
solution of the $\mathop{\bf CVaR}$-constrained problem (10). However, with
the use of the safeguards $[\lambda^{\rm{min}}$, $\lambda^{\rm{max}}]$, we can
guarantee that the limit point of the augmented Lagrangian algorithm is at
least a stationary point of constraint violation (Birgin and Martínez,, 2009,
Theorem 6.3). More specifically, below we show that the algorithm finds a
feasible point of condition (9) under mild assumptions, and otherwise
converges to a stationary point of the constraint violation function
$(1/2)\|H(\alpha,\theta)+s-\kappa\|^{2}$. The proof is delayed to Appendix D.

######  Theorem 3.2.

Suppose Assumptions 3.1 to 3.5 hold. If the sequence of penalty parameters
$\\{\mu^{k}\\}$ is bounded, then the sequence of updates $\\{\theta^{k}\\}$ in
Algorithm 1 converges almost surely to a feasible point of the
$\widehat{\mathop{\bf CVaR}}$ condition (9) as $k\to\infty$. Otherwise,
$\\{\alpha^{k},s^{k},\theta^{k}\\}$ converges to a stationary point of the
constraint violation function $(1/2)\|H(\alpha,\theta)+s-\kappa\|^{2}$.

Note that a feasible point $\theta^{\star}$ of the $\widehat{\mathop{\bf
CVaR}}$ condition (9) also implies a feasible point of the contrained problem
(10). In practice, with a finite number of iterations and some $\epsilon>0$,
we obtain an $\epsilon$-stationary solution of (10) satisfying condition (12).

##  4 Probabilistic guarantees of constraint satisfaction

In the previous section, we have proved convergence to a feasible solution of
the sampled problem (10), which does not exactly correspond to the
probabilistic guarantees of constraint satisfaction (7). To arrive at our
desired guarantees, where the inner probability is with respect to the unknown
joint distribution ${\mathbf{P}}_{(u,y)}$, we need to analyze the
generalization behavior of the robust solutions on out-of-sample data. We
begin with a boundedness assumption on the function inside the expectation
from the $\mathop{\bf CVaR}$ definition, which appears in (8).

######  Assumption 4.1.

The function $g(x(\theta,y),u,y)-\alpha)_{+}/\eta+\alpha$ maps to values in
$(-c^{\rm min},c^{\rm max})$, where $0\leq c^{\rm min},c^{\rm max}<\infty$,
for all $\theta,\alpha,u,y$ in its support.

Let $C=c^{\rm min}+c^{\rm max}$ and define the normalized function

$\psi_{z}(u,y)=\frac{1}{C}\left(c^{\rm
min}+\frac{(g(x,u,y)-\alpha)_{+}}{\eta}+\alpha\right),$ (13)

where $z=(x,\alpha)\in\mathcal{Z}\subset{\mbox{\bf R}}^{n+1}$. Note that when
Assumption 3.4 holds, we have
$\mathcal{Z}=\mathcal{X}\times\mathcal{D}_{\alpha}$ being a compact set. The
function $\psi_{z}$ then maps to values in $(0,1)$, and we consider the entire
function class

$\Psi=\\{\psi_{z}\\}_{z\in\mathcal{Z}},$

which we use to acquire the finite sample probabilistic guarantee in Theorem
4.1.

######  Theorem 4.1 (Finite sample probabilistic guarantee).

Suppose Assumptions 3.4 and 4.1 and hold. For the function class $\Psi$
defined above, suppose the following condition on the covering number $N$
holds with constants $V$ and $K$, relative to the $L_{2}(Q)$ norm
$\|\psi\|_{Q,2}=(\int|\psi|^{2}dQ)^{1/2}$ with probability measure $Q$,

$\sup_{Q}N(\delta,\Psi,L_{2}(Q))\leq\left(\frac{K}{\delta}\right)^{V},$ (14)

for every $0<\delta<K$. Then, when Algorithm 1 converges to an
$\epsilon$-stationary point of (10), we have the finite sample probabilistic
guarantee

${\mathbf{P}}^{N\times J}\left({\mathbf{P}}_{(u,y)}(g(x,u,y)\leq 0)\geq
1-\eta\quad\forall x\in\mathcal{X}\right)\geq 1-\beta,$ (15)

where $\beta=(D\sqrt{N}C\tau/\sqrt{V})^{V}\exp(-2NC^{2}\tau^{2})$, $D$ is a
constant that depends only on $K$, and $\tau\leq(c^{\min}+\kappa+\epsilon)/C$.

###### Proof.

By the definition of the function class $\Psi$ and Assumptions 3.4 and 4.1, we
can apply Theorem 2.14.9 by van der Vaart and Wellner, (1996) to get

${\mathbf{P}}^{N\times
J}\left(\sup_{\psi_{z}\in\Psi}\left\|\frac{1}{NJ}\sum_{j=1}^{J}\sum_{i=1}^{N}\psi_{z}(u^{i},y^{j})-{\mbox{\bf
E}}_{(u,y)}[\psi_{z}(u,y)]\right\|\leq\tau\right)\geq
1-\left(\frac{D\sqrt{N}C\tau}{\sqrt{V}}\right)^{V}\exp(-2NC^{2}\tau^{2}).$

By the convergence of Algorithm 1 to an $\epsilon$-stationary point of (10),
we have that the empirical $\mathop{\bf CVaR}$ will be close to $\kappa$,
which, given equation (13), means

$\frac{1}{NJ}\sum_{j=1}^{J}\sum_{i=1}^{N}\psi_{z}(u^{i},y^{j})\leq(1/C)(c^{\min}+\kappa+\epsilon).$

Therefore, we have

${\mathbf{P}}^{N\times J}\left(\sup_{\psi_{z}\in\Psi}{\mbox{\bf
E}}_{(u,y)}[\psi_{z}(u,y)]\leq(1/C)(c^{\min}+\kappa+\epsilon)+\tau\right)\geq
1-\beta,$

where $\beta=(D\sqrt{N}C\tau/\sqrt{V})^{V}\exp(-2NC^{2}\tau^{2})$. If $\kappa$
is set such that $\kappa\leq-C\tau-\epsilon-c^{\min}$, we then have

${\mathbf{P}}^{N\times J}\left(\mathop{\bf CVaR}(g(x,u,y),\eta)\leq
0\quad\forall x\in\mathcal{X}\right)\geq 1-\beta,$

which implies

${\mathbf{P}}^{N\times J}\left({\mathbf{P}}_{(u,y)}(g(x,u,y)\leq 0)\geq
1-\eta\quad\forall x\in\mathcal{X}\right)\geq 1-\beta.$ (16)

∎

As this probabilistic guarantee holds for all $x$ within the bounded set
$\mathcal{X}$, it also holds for our data-driven solutions $x(\theta,y)$.

######  Corollary 4.1.1.

Suppose Assumption 3.4 holds. For all $\theta^{0}\in\mathcal{D}_{\theta}$,
$s^{0}\in\mathcal{D}_{s}$, and $y\in\mathcal{D}_{y}$, Theorem 4.1 implies the
probabilistic guarantee (7) for the data-driven solution
$x(\theta^{\star},y)$, i.e., ${\mathbf{P}}^{N\times
J}\left\\{{\mathbf{P}}_{(u,y)}(g(x(\theta^{\star},y),u,y)\leq 0)\geq
1-\eta\right\\}\geq 1-\beta.$

We note that the covering number requirement holds when $\Psi$ is a class of
Lipschitz functions, which is implied when the constraint function $g$ is
Lipschitz. For instance, when $g$ is affine or a maximum of affine functions
in variable $x$, it is Lipschitz in $x$ for all $u\in\mathcal{D}_{u}$ and
$y\in\mathcal{D}_{y}$, i.e.,

$|g(x,u,y)-g(x^{\prime},u,y)|\leq\bar{L}\|x-x^{\prime}\|_{2}\quad\forall
x,x^{\prime}\in\mathcal{X},~{}u\in\mathcal{D}_{u},~{}y\in\mathcal{D}_{y},$

with some Lipschitz constant $\bar{L}$. Then, $\psi_{z}\in\Psi$ is Lipschitz
in $z=(x,\alpha)$ with Lipschitz constant $L\leq(1/C)((\bar{L}+1)/\eta+1)$. We
can now derive the following covering number condition.

##### Covering number for $L$-Lipschitz functions.

When $\psi_{z}$ is $L$-Lipschitz in $z$ for all $u$ and $y$, we can apply
standard bounds on the covering number of Lipschitz losses (van der Vaart and
Wellner,, 1996; Pollard,, 1890; Namkoong and Duchi,, 2017). In this case, the
covering number condition (14) reduces to

$\sup_{Q}N(\delta,\Psi,L_{2}(Q))\leq
N(\delta/L,\mathcal{Z},\|\cdot\|_{2})\leq\left(1+\frac{(x^{\rm
max}+\alpha^{\rm
max})L}{\delta}\right)^{n+1}\leq\left(\frac{K}{\delta}\right)^{n+1},$

for every $0<\delta<K$, where $K=\delta+(x^{\rm max}+\alpha^{\rm max})L$.
Then, $\beta=(D\sqrt{N}C\tau/\sqrt{n+1})^{n+1}\exp(-2NC^{2}\tau^{2})$.

##  5 Examples

We run the experiments on the Princeton Institute for Computational Science
and Engineering (PICSciE) facility with 20 parallel 2.4 GHz Skylake cores. We
solve the inner optimization problems (1) with the Clarabel optimizer with
default settings. The learning procedure relies on the Python package
cvxpylayers (Agrawal et al., 2019a, ), and the code to replicate the
experiments can be accessed at

https://github.com/stellatogrp/lropt_experiments/.

For each experiment, we compare our method with a RO and a DRO benchmark,
given below.

Learning for robust optimization (LRO) (our method).

We apply the training method given in Section 3.2, for a variety of $\eta$
values, to achieve various empirical probabilities of constraint violation
$\hat{\eta}$. The uncertainty set can be RO or DRO based. Default training
parameters are given in Table 2.

LRO \+ Fine Tuning (our method).

Upon completion of the LRO training procedure, we also fine tune the size of
the uncertainty set using grid search to precisely achieve a target in-sample
empirical probability of constraint violation.

Mean-Variance RO.

We perform mean-variance RO, where the shape of the uncertainty set is
determined by the empirical mean and variance of the data (Ben-Tal and
Nemirovski,, 2000). We use grid search over a the set size to achieve various
empirical probabilities of constraint violation $\hat{\eta}$.

Wasserstein DRO.

We perform Wasserstein-1 DRO (Esfahani and Kuhn,, 2018), where the uncertainty
set is constructed around individual data points. We again use grid search
over the set size to achieve various empirical probabilities of constraint
violation $\hat{\eta}$.

For each setting, we repeat the experiment 20 times, each with a new
independent dataset. For each method, we evaluate the performance metrics in
Table 1.

Table 1: Performance metrics in numerical examples, averaged over 20 runs.

Name Description Obj. Out-of-sample objective value estimate for problem (8)
$\hat{\eta}$ Out-of-sample empirical probability of constraint violation
$\hat{\beta}$ Out-of-sample estimated confidence, given as
${\mathbf{P}}(\hat{\eta}>\eta)$, where the probability is over the testing
data ${\rm time}$ Average solution time of the final robust problem after
training

Table 2: Default hyperparameter values for Algorithm 1. We set the initial uncertainty set parameter values $(A^{0},b^{0})=(\hat{\Sigma}^{1/2},\hat{\mu})$, depending on the empirical covariance and mean of the training data. $\theta^{0}$ | $\alpha^{0}$ | $s^{0}$ | $\lambda^{0}$ | $\mu^{0}$ | $\kappa$ | $\sigma$ | $\tau$ | $k_{\text{max}}$ | $t_{\text{max}}$ | $\eta$ | $|B|$
---|---|---|---|---|---|---|---|---|---|---|---
$(A^{0},b^{0})$ | 0 | 0 | 1 | 1 | $-0.01$ | 1.015 | 0.95 | 1000 | 100 | $[0,0.3]$ | 30

###  5.1 Portfolio optimization

##### Problem description.

We consider larger instances of the portfolio optimization problem from
Section 2.1, where we select a portfolio $x\in{\mbox{\bf R}}^{n}$ of stocks to
maximize returns, and refrain from deviating too far from the previous
holdings $x^{\text{prev}}\in{\mbox{\bf R}}^{n}$. As we know these previous
holding values, we parametralize the problem with the family parameter
$y=x^{\text{prev}}$, and solve the RO problem

$\begin{array}[]{ll}\mbox{minimize}&t+\lambda^{\text{trade}}\|x-x^{\text{prev}}\|_{1}\\\
\mbox{subject to}&-u^{T}x\leq t\quad\forall u\in{\mathcal{U}}(\theta)\\\
&\mathbf{1}^{T}x=1,\quad x\geq 0,\end{array}$

where the uncertain parameter $u$ stands for the unknown returns.

##### Data.

We consider the case where $n=5$ and $n=10$ stocks. We let
$\lambda^{\text{trade}}=0.2$. The uncertain returns $u$ is sampled from a
multivariate normal distribution $N(\mu,\Sigma)$, where $\mu$ is generated
uniformly on $[0.5,1]^{n}$, and $\Sigma$ is generated as $0.1FF^{T}$, where
$F$ is generated on $[0,1]^{n\times 2}$. The nominal previous holdings
$\bar{x}^{\text{prev}}$ is generated from a Dirichlet distribution with
parameter values equal to $\mu$. We consider 10 instances of
$x^{\text{prev}}$, distributed normally around $\bar{x}^{\text{prev}}$, with a
variance of $0.05$. Each sample is then projected to be non-negative, and
normalized to sum to 1. The total training sample size is $N=500$, with a
testing sample size of $20000$.

##### Uncertainty set.

We learn both a RO and a DRO uncertainty set. For RO, the set is ellipsoidal
with with norm $p=2$, and for DRO, we train the shape of a mean robust
optimization (MRO) set with $K=5$ and power $p=1$. This set is closely related
to the Wasserstein type-1 DRO set. See Appendix A for details on the MRO
uncertainty set. By choosing $K\ll N$, we speed up the training procedure with
minimal performance degradation compared to training with Wasserstein DRO
(Wang et al.,, 2022). For both $n=5$ and $n=10$, we learn uncertainty sets for
$\eta\in\eta^{\text{list}}=(0.01,0.03,0.05,0.08,0.10,0.15,0.2,0.3)$.

##### Results.

In Table 3, we compare the performance of our LRO methods with mean-variance
RO and Wasserstein DRO, and note the desirable trade-off given by our learned
sets.

Table 3: Portfolio optimization comparisons averaged across 20 runs.
LRO${}_{\text{RO}}$: learning the shape and size of the RO uncertainty set
with $\eta=0.3$. LRO-T${}_{\text{RO},0.03}$: Best value given by
LRO${}_{\text{RO}}$ with $\eta\in\eta^{\text{list}}$, tuned for an in-sample
empirical probability of constraint violation ($\hat{\eta}^{\text{in}}$) below
0.03. LRO${}_{\text{MRO}}$: learning the shape and size of the MRO uncertainty
set with $\eta=0.3$. LRO-T${}_{\text{DRO},0.03}$: Best value given by
LRO${}_{\text{MRO}}$ with $\eta\in\eta^{\text{list}}$, tuned for an in-sample
empirical probability of constraint violation ($\hat{\eta}^{\text{in}}$) below
0.03. Mean-Var RO0.03: mean-variance RO where the size is tuned for an
$\hat{\eta}^{\text{in}}$ below $0.03$. W-DRO0.03: Wasserstein DRO where the
size is tuned for an $\hat{\eta}^{\text{in}}$ below $0.03$. Obj.: the average
out-of-sample objective value ${\mbox{\bf
E}}[-u^{T}x+0.2\|x-x^{\text{prev}}\|_{1}]$. ${\rm time}$: average solve time
after training.

$n=5$ (dimension of $u$) Method LRO${}_{\text{RO}}$ LRO-T${}_{\text{RO},0.03}$
LRO${}_{\text{MRO}}$ LRO-T${}_{\text{MRO},0.03}$ MV-RO0.03 W-DRO0.03 [head to
column names]table_vals2.csv $n=10$ (dimension of $u$) [head to column
names]table_vals2.csv

In Figure 3, for the 20 runs of the experiment, and for $m\in\\{5,10\\}$, we
show the trade-off between the average objective value and empirical
probability of constraint violation. The trade-offs for the LRO sets are
obtained by applying the training algorithm for the $\eta$ values provided
above, fine-tuning the size of the uncertainty set, and reporting the best
objective values for a range of $\hat{\eta}$ values. We observe lower
objective values for the LRO sets across the relevant levels of constraint
violation. Note that in general, at the same out-of-sample objective value,
the mean-variance RO and nominal Wasserstein DRO sets suffer from a much
higher empirical probability of constraint violation. Therefore, the decision-
focused uncertainty sets are more effective than uncertainty sets based on the
distribution alone.

Comparing the LRO-RO and LRO-MRO sets, we note that the RO-based set has a
higher variance in performance. However, as seen for $n=10$, it could
outperform the MRO-based set. In terms of computation time, while the RO
problems are less expensive to solve, the MRO framework (Wang et al.,, 2022)
already presents a significant computational time reduction compared to
Wasserstein DRO, especially when $K\ll N$.

In Figure 4, for $n=10$ stocks, a single $y$, and a single run of the
experiment, we show the distribution of the optimal portfolio weights for all
four uncertainty sets. The mean of the uncertain returns is
$\mu=(0.5,0.55,0.57,0.59,0.65,0.67,0.7,0.71,0.77,0.86)$, sorted in ascending
order (bottom to top in Figure 4), and the corresponding diagonals of $\Sigma$
are $(0.64,0.08,0.25,0.03,0.06,0.25,0.11,0.36,0.13,0.02)$. The covariances
between the entries are omitted here for brevity. To mitigate the risks, when
$\hat{\eta}$ is low the distribution of portfolio weights is either skewed
towards the relatively lower variance stocks, as is the case for the the RO
sets, or is close to an even distribution, as is the case for the DRO sets.
Intuitively, one would invest more in stocks with higher returns and higher
risks when the proxy of risk tolerance, $\hat{\eta}$, is high. Indeed, we
observe that the relatively high-risk, high-return stock, the third-to-last
stock (depicted in orange), receives larger investments as $\hat{\eta}$
increases. Note, however, that the reshaped sets, LRO-RO and LRO-MRO, are able
to identify this profitable stock and increase its portfolio weight even at
relatively low values of $\hat{\eta}$. By identifying a more compact region of
uncertainty targeted for the problem family, they are able reduce conservatism
in the portfolio selection process and take more profitable positions.

Figure 3: Portfolio. Average trade-offs between out-of-sample objective values and probability of constraint violation across 20 experiments. The shaded region represents the 0.3 to 0.7 quantiles. The LRO results are obtained through training procedures with $\eta$ values in the range $[0,0.3]$.  |
---|---
Figure 4: Portfolio. Optimal portfolio weights chosen by the four methods for
various $\hat{\eta}$. The colors represent the true mean $\mu$ of the returns.

###  5.2 Inventory management

##### Problem description.

We consider a two-stage adjustable RO problem on inventory management (Iancu
and Trichakis,, 2014, Section 5.2.1). There exist a retail network consisting
of a single warehouse and $n$ different retail points, indexed by
$i=1,\dots,n$. For simplicity, only a single product is sold. There are a
total of $C$ units of the product available, and each retail point starts with
$0$ inventory and is capable of stocking at most $c_{i}$ units. The
transportation costs for distributing the items at the $i$th retail point is
$t_{i}$ currency units per unit of inventory, and the operating cost is
$h_{i}$ per unit. The revenues per unit is $r_{i}$, which we take as our
family parameter, ${\it i.e.},~{}y=r$, as these values are known but may be
subject to minor changes across problem instances. Customer demand is
uncertain, with a nominal demand plus an additional amount dependent on $m$
market factors. The formula for the demand at retail point $i$ is

$a_{i}=\bar{a}_{i}+q_{i}^{T}u,$

where $\bar{a}_{i}\in{\mbox{\bf R}}$ denotes the nominal demand, and
$q_{i}\in{\mbox{\bf R}}^{m}$ denotes the exposure of demand $a_{i}$ to the
market factors. We denote with $u$ the unknown market factors, which belong to
an uncertainty set ${\mathcal{U}}(\theta)$. The decision variable is the
stocking decisions for all retail points, denoted $s\in{\mbox{\bf R}}^{n}$. In
addition, the problem involves a vector $w(u)$ of realized sales at each
point, which is dependent on stocking decision $s$ and the unknown factors $u$
affecting the demand. We have

$w_{i}(u)=\min\\{s_{i},a_{i}\\},$

for each retail point. The manager makes stocking decisions to maximize worst-
case profits, which is equivalent to minimizing worst-case loss,
$-r^{T}w(u)+(t+h)^{T}s$, $\forall u\in{\mathcal{U}}(\theta)$. For each problem
instance, the RO formulation is given as

$\begin{array}[]{ll}\mbox{minimize}&\tau\\\ \mbox{subject
to}&-r^{T}w(u)+(t+h)^{T}s\leq\tau,\quad\forall u\in{\mathcal{U}}(\theta)\\\
&w_{i}(u)\leq s_{i},\quad i=1,\dots,n,\quad\forall
u\in{\mathcal{U}}(\theta)\\\ &w_{i}(u)\leq\bar{a}_{i}+q_{i}^{T}u,\quad
i=1,\dots,n,\quad\forall u\in{\mathcal{U}}(\theta)\\\ &\mathbf{1}^{T}s=C,\quad
0\leq s\leq c,\\\ \end{array}$ (17)

where variables $s\in{\mbox{\bf R}}^{n}$, $w\in{\mbox{\bf R}}^{n}$, and
$\tau\in{\mbox{\bf R}}$. To get a tractable approximation, we replace the
adjustable decision $w_{i}(u)$ with their affine adjustable robust
counterparts $z^{0}_{i}+v_{i}^{T}u$, where $z^{0}_{i}\in{\mbox{\bf R}}$ and
$v_{i}\in{\mbox{\bf R}}^{f}$ are auxiliary introduced as part of the linear
decision rules (Ben-Tal et al.,, 2009). The new formulation is

$\begin{array}[]{ll}\mbox{minimize}&\tau\\\ \mbox{subject
to}&-r^{T}z^{0}-r^{T}Vu+(t+h)^{T}s\leq\tau,\quad\forall
u\in{\mathcal{U}}(\theta)\\\ &z^{0}_{i}+v_{i}^{T}u\leq s_{i},\quad
i=1,\dots,n,\quad\forall u\in{\mathcal{U}}(\theta)\\\
&z^{0}_{i}+v_{i}^{T}u\leq\bar{a}_{i}+q_{i}^{T}u,\quad i=1,\dots,n,\quad\forall
u\in{\mathcal{U}}(\theta)\\\ &\mathbf{1}^{T}s=C,\quad 0\leq s\leq c,\\\
\end{array}$ (18)

where $V=[v_{1}\cdots v_{n}]^{T}\in{\mbox{\bf R}}^{n\times f}$. We treat the
multiple uncertain constraints as a single joint chance constraint represented
as the maximum of concave functions.

##### Data.

As in the work by Iancu and Trichakis, (2014), we let $n=10$. We consider the
cases where the number of uncertain market factors is $m=4$ and $m=8$. We set
$C=200$, and all other problem data are independently sampled from uniformly
distributed random variables. Inventory capacities $c$ are between 30 and 50.
Transportation and operation costs are between 1 and 3. Sales revenues $r$ are
between 2 and 4, and nominal demand values are between 10 and 20. We consider
10 instances for $r$, perturbing each by random sample of a normal
distribution with mean 0 and standard deviation $0.05$. The market factor
exposure coefficients $q_{i}$ are between $-2$ and $2$. We generate uncertain
market factors $u$ from 3 normal distributions, to simulate 3 uncertain market
environments. For distribution $j$, the mean is $\mu_{i}=i\gamma_{j}$ for all
$i=1,\dots,m$, with $\gamma=(0.3,0.6,0.9)$, and the variance is
$\sigma_{i}=0.02^{2}+(0.1i)^{2}$ for all distributions. The total training
sample size is $N=100$, split evenly across the 3 distributions. The testing
sample size is $20000$.

##### Uncertainty set.

We train the shape and size of an ellipsoidal uncertainty set, with norm
$p=2$, and for
$\eta\in\eta^{\text{list}}=(0.01,0.03,0.05,0.08,0.10,0.15,0.2,0.3)$.

##### Results.

In Table 4, we compare the performance of our LRO method with mean-variance RO
and Wasserstein DRO, and note the desirable trade-off given by our learned
reshaped set.

Table 4: Average inventory management performance over 20 runs. LRO: learning
the shape and size of the uncertainty set with $\eta=0.3$. LRO-T0.03: Best
value given by LRO with $\eta\in\eta^{\text{list}}$, tuned for an in-sample
empirical probability of constraint violation ($\hat{\eta}^{\text{in}}$) below
0.03. LRO-T0.05: Best value given by LRO with $\eta\in\eta^{\text{list}}$,
tuned for an $\hat{\eta}^{\text{in}}$ below 0.05. Mean-Var RO0.03: mean-
variance RO where the size is tuned for an $\hat{\eta}^{\text{in}}$ below
$0.03$. Mean-Var RO0.05: mean-variance RO where the size is tuned for an
$\hat{\eta}^{\text{in}}$ below $0.05$. W-DRO0.03: Wasserstein DRO where the
size is tuned for an $\hat{\eta}^{\text{in}}$ below $0.03$. W-DRO0.05:
Wasserstein DRO where the size is tuned for an $\hat{\eta}^{\text{in}}$ below
$0.05$. Obj.: the average out-of-sample objective value ${\mbox{\bf
E}}[-r^{T}w(u)+(t+h)^{T}s]$. ${\rm time}$: average solve time after training.

$m=4$ (dimension of $u$) Method LRO LRO-T0.03 LRO-T0.05 MV-RO0.03 MV-RO0.05
W-DRO0.03 W-DRO0.05 [head to column names]table_vals_inv.csv $m=8$ (dimension
of $u$) [head to column names]table_vals_inv.csv

In Figure 5, we show the trade-off between the average objective value and
empirical probability of constraint violation for the 20 runs of the
experiment. We observe an improved objective value for the LRO set across
relevant levels of constraint violation. In Figure 5, for $m=4$, we note that
the LRO set with a constraint violation $\hat{\eta}$ of around 0.028 achieves
an objective value that the mean-variance RO and Wasserstein DRO sets only
achieve at signficantly higher $\hat{\eta}$’s, of around 0.12 and 0.07,
respectively. A similar comparison holds for $m=8$. It is apparent that our
decision-focused manner of learning uncertainty sets identifies more effective
regions of uncertainty and therefore reduces the conservatism of the data-
driven solutions.

Figure 5: Inventory. Average trade-offs between out-of-sample objective values
and probability of constraint violation across 20 experiments. The shaded
region represents the 0.3 to 0.7 quantiles. The LROPT results are obtained
through training procedures with $\eta$ values in the range $[0,0.3]$.

##  6 Conclusion

We develop a data-driven method to learn the shape and size uncertainty sets
in robust optimization in a decision-focused manner, i.e., using the structure
and the solution of the robust counterpart. By formulating the learning
problem as a stochastic optimization problem with $\mathop{\bf CVaR}$
constraints, we choose the shape and size of the uncertainty sets to optimize
the data-driven performance while enforcing desired probabilistic guarantees.
We show convergence of our algorithm using nonsmooth implicit differentiation,
and provide finite sample probabilistic guarantees on unseen data. Through
numerical examples, we demonstrate the performance gain achieved by our method
compared to traditional manual tuning methods.

##  7 Acknowledgments

Irina Wang and Bartolomeo Stellato are supported by the NSF CAREER Award ECCS
2239771. The simulations presented in this article were performed on
computational resources managed and supported by Princeton Research Computing,
a consortium of groups including the Princeton Institute for Computational
Science and Engineering (PICSciE) and the Office of Information Technology’s
High Performance Computing Center and Visualization Laboratory at Princeton
University.

## References

  * (1) Agrawal, A., Amos, B., Barratt, S., Boyd, S., Diamond, S., and Kolter, Z. (2019a).  Differentiable convex optimization layers.  In Advances in Neural Information Processing Systems.
  * (2) Agrawal, A., Barratt, S., Boyd, S., Busseti, E., and Moursi, W. (2019b).  Differentiating through a cone program.  Journal of Applied and Numerical Optimization, 1(2):107–115.
  * Amos and Kolter, (2017) Amos, B. and Kolter, J. Z. (2017).  Optnet: differentiable optimization as a layer in neural networks.  In Proceedings of the 34th International Conference on Machine Learning - Volume 70, ICML’17, page 136–145. JMLR.org.
  * Bandi and Bertsimas, (2012) Bandi, C. and Bertsimas, D. (2012).  Tractable stochastic analysis in high dimensions via robust optimization.  Mathematical Programming, 134(1):23–70.
  * Bank et al., (1984) Bank, B., Guddat, J., Klatte, D., Kummer, B., and Tammer, K. (1984).  Non-linear parametric optimization.  SIAM Review, 26(4):594–595.
  * Bauschke and Combettes, (2017) Bauschke, H. H. and Combettes, P. L. (2017).  Convex Analysis and Monotone Operator Theory in Hilbert Spaces.  Springer Publishing Company, Incorporated, 2nd edition.
  * Ben-Tal et al., (2009) Ben-Tal, A., El Ghaoui, L., and Nemirovski, A. (2009).  Robust Optimization.  Princeton University Press.
  * Ben-Tal and Nemirovski, (2000) Ben-Tal, A. and Nemirovski, A. (2000).  Robust solutions of Linear Programming problems contaminated with uncertain data.  Mathematical Programming, 88(3):411–424.
  * Ben-Tal and Nemirovski, (2008) Ben-Tal, A. and Nemirovski, A. (2008).  Selected topics in robust convex optimization.  Math. Program., 112:125–158.
  * Bertrand et al., (2020) Bertrand, Q., Klopfenstein, Q., Blondel, M., Vaiter, S., Gramfort, A., and Salmon, J. (2020).  Implicit differentiation of Lasso-type models for hyperparameter optimization.  In Proceedings of the 37th International Conference on Machine Learning, volume 119 of Proceedings of Machine Learning Research, pages 810–821. PMLR.
  * Bertrand et al., (2022) Bertrand, Q., Klopfenstein, Q., Massias, M., Blondel, M., Vaiter, S., Gramfort, A., and Salmon, J. (2022).  Implicit differentiation for fast hyperparameter selection in non-smooth convex learning.  The Journal of Machine Learning Research, 23(1):6680–6722.
  * Bertsimas et al., (2011) Bertsimas, D., Brown, D. B., and Caramanis, C. (2011).  Theory and applications of robust optimization.  SIAM Review, 53(3):464–501.
  * Bertsimas and den Hertog, (2022) Bertsimas, D. and den Hertog, D. (2022).  Robust and Adaptive Optimization.  Dynamic Ideas.
  * Bertsimas et al., (2021) Bertsimas, D., den Hertog, D., and Pauphilet, J. (2021).  Probabilistic Guarantees in Robust Optimization.  SIAM Journal on Optimization, 31(4):2893–2920.
  * Bertsimas et al., (2018) Bertsimas, D., Gupta, V., and Kallus, N. (2018).  Data-driven robust optimization.  Mathematical Programming, 167(2):235–292.
  * Bertsimas et al., (2022) Bertsimas, D., Shtern, S., and Sturt, B. (2022).  Technical note—two-stage sample robust optimization.  Operations Research, 70(1):624–640.
  * Birgin and Martínez, (2009) Birgin, E. G. and Martínez, J. M. (2009).  Practical augmented Lagrangian methodsPractical Augmented Lagrangian Methods, pages 3013–3023.  Springer US, Boston, MA.
  * Blanchet et al., (2019) Blanchet, J., Kang, Y., Murthy, K., and Zhang, F. (2019).  Data-driven optimal transport cost selection for distributionally robust optimization.  2019 Winter Simulation Conference (WSC).
  * Bolte et al., (2021) Bolte, J., Le, T., Pauwels, E., and Silveti-Falls, T. (2021).  Nonsmooth implicit differentiation for machine-learning and optimization.  In Advances in Neural Information Processing Systems, volume 34, pages 13537–13549.
  * Bolte and Pauwels, (2019) Bolte, J. and Pauwels, E. (2019).  Conservative set valued fields, automatic differentiation, stochastic gradient method and deep learning.  Mathematical Programming, 188.
  * Busseti et al., (2019) Busseti, E., Moursi, W., and Boyd, S. (2019).  Solution refinement at regular points of conic problems.  Computational Optimization and Applications, 74(3):627–643.
  * Cameron et al., (2022) Cameron, C., Hartford, J., Lundy, T., and Leyton-Brown, K. (2022).  The perils of learning before optimizing.  Proceedings of the AAAI Conference on Artificial Intelligence, 36:3708–3715.
  * Chenreddy et al., (2022) Chenreddy, A., Bandi, N., and Delage, E. (2022).  Data-driven conditional robust optimization.  In Advances in Neural Information Processing Systems, volume 35, pages 9525–9537. Curran Associates, Inc.
  * Chenreddy and Delage, (2024) Chenreddy, A. and Delage, E. (2024).  End-to-end conditional robust optimization.  arXiv preprint arXiv:2403.04670.
  * Chow et al., (2015) Chow, Y., Ghavamzadeh, M., Janson, L., and Pavone, M. (2015).  Risk-constrained reinforcement learning with percentile risk criteria.  Journal of Machine Learning Research, 18.
  * Clarke, (1975) Clarke, F. H. (1975).  Generalized gradients and applications.  Transactions of the American Mathematical Society, 205:247–262.
  * Costa and Iyengar, (2023) Costa, G. and Iyengar, G. N. (2023).  Distributionally robust end-to-end portfolio construction.  Quantitative Finance, 23(10):1465–1482.
  * Coste, (2000) Coste, M. (2000).  An Introduction to O-minimal Geometry.  Dottorato di ricerca in matematica / Università di Pisa, Dipartimento di Matematica. Istituti editoriali e poligrafici internazionali.
  * Davis et al., (2018) Davis, D., Drusvyatskiy, D., Kakade, S. M., and Lee, J. (2018).  Stochastic subgradient method converges on tame functions.  Foundations of Computational Mathematics, 20:119–154.
  * Demirović et al., (2019) Demirović, E., Stuckey, P., Bailey, J., Chan, J., Leckie, C., Ramamohanarao, K., and Guns, T. (2019).  An investigation into prediction + optimisation for the knapsack problem.  In Rousseau, L.-M. and Stergiou, K., editors, Integration of Constraint Programming, Artificial Intelligence, and Operations Research, Lecture Notes in Computer Science, pages 241–257. Springer.
  * Diamond and Boyd, (2016) Diamond, S. and Boyd, S. (2016).  CVXPY: A Python-embedded modeling language for convex optimization.  Journal of Machine Learning Research, 17(83):1–5.
  * Donti et al., (2017) Donti, P. L., Amos, B., and Kolter, J. Z. (2017).  Task-based end-to-end model learning in stochastic optimization.  In Proceedings of the 31st International Conference on Neural Information Processing Systems, NIPS’17, page 5490–5500, Red Hook, NY, USA. Curran Associates Inc.
  * Dunning et al., (2017) Dunning, I., Huchette, J., and Lubin, M. (2017).  JuMP: A Modeling Language for Mathematical Optimization.  SIAM Review, 59(2):295–320.
  * Elmachtoub and Grigas, (2022) Elmachtoub, A. N. and Grigas, P. (2022).  Smart “predict, then optimize”.  Manage. Sci., 68(1):9–26.
  * Esfahani and Kuhn, (2018) Esfahani, P. M. and Kuhn, D. (2018).  Data-driven distributionally robust optimization using the Wasserstein metric: performance guarantees and tractable reformulations.  Mathematical Programming, 171:115–166.
  * Gao, (2023) Gao, R. (2023).  Finite-sample guarantees for wasserstein distributionally robust optimization: Breaking the curse of dimensionality.  Operations Research, 71(6):2291–2306.
  * Gao and Kleywegt, (2023) Gao, R. and Kleywegt, A. (2023).  Distributionally robust stochastic optimization with wasserstein distance.
  * Goerigk and Kurtz, (2023) Goerigk, M. and Kurtz, J. (2023).  Data-driven robust optimization using deep neural networks.  Computers and Operations Research, 151:106087.
  * Goh and Sim, (2009) Goh, J. and Sim, M. (2009).  Robust optimization made easy with rome.  Operations Research, 59:973–985.
  * Gupta, (2019) Gupta, V. (2019).  Near-optimal bayesian ambiguity sets for distributionally robust optimization.  Management Science, 65(9):4242–4260.
  * Guzman et al., (2017) Guzman, Y. A., Matthews, L. R., and Floudas, C. A. (2017).  New a priori and a posteriori probabilistic bounds for robust counterpart optimization: Iii. exact and near-exact a posteriori expressions for known probability distributions.  Computers & Chemical Engineering, 103:116–143.
  * Hong et al., (2021) Hong, L. J., Huang, Z., and Lam, H. (2021).  Learning-based robust optimization: Procedures and statistical guarantees.  Management Science, 67(6):3447–3467.
  * Iancu and Trichakis, (2014) Iancu, D. A. and Trichakis, N. (2014).  Pareto efficiency in robust optimization.  Management Science, 60(1):130–147.
  * Juditsky and Nemirovski, (2022) Juditsky, A. and Nemirovski, A. (2022).  On well-structured convex–concave saddle point problems and variational inequalities with monotone operators.  Optimization Methods and Software, 37(5):1567–1602.
  * Kanzow et al., (2009) Kanzow, C., Ferenczi, I., and Fukushima, M. (2009).  On the local convergence of semismooth newton methods for linear and nonlinear second-order cone programs without strict complementarity.  SIAM Journal on Optimization, 20(1):297–320.
  * Kotary et al., (2021) Kotary, J., Fioretto, F., Van Hentenryck, P., and Wilder, B. (2021).  End-to-end constrained optimization learning: A survey.  International Joint Conference on Artificial Intelligence, pages 4475–4482.
  * Li et al., (2012) Li, Z., Tang, Q., and Floudas, C. A. (2012).  A comparative theoretical and computational study on robust counterpart optimization: Ii. probabilistic guarantees on constraint satisfaction.  Industrial & Engineering Chemistry Research, 51(19):6769–6788.
  * Mafusalov and Uryasev, (2018) Mafusalov, A. and Uryasev, S. (2018).  Buffered probability of exceedance: Mathematical properties and optimization.  SIAM Journal on Optimization, 28(2):1077–1103.
  * Namkoong and Duchi, (2017) Namkoong, H. and Duchi, J. C. (2017).  Variance-based regularization with convex objectives.  In Advances in Neural Information Processing Systems, volume 30. Curran Associates, Inc.
  * Pollard, (1890) Pollard, D. (1890).  Convergence of Stochastic Processes.  Springer New York, NY.
  * Rockafellar and Royset, (2010) Rockafellar, R. and Royset, J. (2010).  On buffered failure probability in design and optimization of structures.  Reliability Engineering & System Safety, 95(5):499–510.
  * Rockafellar and Uryasev, (2002) Rockafellar, R. T. and Uryasev, S. (2002).  Conditional value-at-risk for general loss distributions.  Journal of banking & finance, 26(7):1443–1471.
  * Rockafellar and Wets, (1998) Rockafellar, R. T. and Wets, R. J. (1998).  Variational analysis.  Grundlehren der mathematischen Wissenschaften.
  * Roos and den Hertog, (2020) Roos, E. and den Hertog, D. (2020).  Reducing conservatism in robust optimization.  INFORMS Journal on Computing, 32(4):1109–1127.
  * Sadana et al., (2023) Sadana, U., Chenreddy, A., Delage, E., Forel, A., Frejinger, E., and Vidal, T. (2023).  A survey of contextual optimization methods for decision making under uncertainty.  arXiv preprint arXiv:2306.10374.
  * Sarykalin et al., (2008) Sarykalin, S., Serraino, G., and Uryasev, S. (2008).  Value-at-Risk vs Conditional Value-at-Risk in Risk Management and Optimization, chapter 13, pages 270–294.  INFORMS.
  * Schiele et al., (2023) Schiele, P., Luxenberg, E., and Boyd, S. (2023).  Disciplined saddle programming.  arXiv preprint arXiv:2301.13427.
  * Schuurmans and Patrinos, (2023) Schuurmans, M. and Patrinos, P. (2023).  Distributionally robust optimization using cost-aware ambiguity sets.  IEEE Control Systems Letters, PP:1–1.
  * Uryasev and Rockafellar, (2001) Uryasev, S. and Rockafellar, R. T. (2001).  Conditional value-at-risk: optimization approach.  In Stochastic optimization: algorithms and applications, pages 411–435. Springer.
  * van den Dries and Miller, (1999) van den Dries, L. and Miller, C. (1999).  Geometric categories, o-minimal structures and control.  In Vaandrager, F. W. and van Schuppen, J. H., editors, Hybrid Systems: Computation and Control, pages 4–4, Berlin, Heidelberg. Springer Berlin Heidelberg.
  * van der Vaart and Wellner, (1996) van der Vaart, A. and Wellner, J. (1996).  Weak Convergence and Empirical Processes: With Applications to Statistics.  Springer New York, NY.
  * Van Parys et al., (2021) Van Parys, B. P., Esfahani, P. M., and Kuhn, D. (2021).  From data to decisions: Distributionally robust optimization is optimal.  Management Science, 67(6):3387–3402.
  * Wang et al., (2022) Wang, I., Becker, C., Van Parys, B., and Stellato, B. (2022).  Mean robust optimization.  arXiv preprint arXiv:2207.10820.
  * Weed and Bach, (2019) Weed, J. and Bach, F. (2019).  Sharp asymptotic and finite-sample rates of convergence of empirical measures in Wasserstein distance.  Bernoulli, 25(4A):2620 – 2648.
  * Wiebe and Misener, (2022) Wiebe, J. and Misener, R. (2022).  Romodel: modeling robust optimization problems in pyomo.  Optimization and Engineering, 23(4):1873–1894.
  * (66) Wilder, B., Dilkina, B., and Tambe, M. (2019a).  Melding the data-decisions pipeline: Decision-focused learning for combinatorial optimization.  Proceedings of the AAAI Conference on Artificial Intelligence, 33(01):1658–1665.
  * (67) Wilder, B., Ewing, E., Dilkina, B., and Tambe, M. (2019b).  End to end learning and optimization on graphs.  In Wallach, H., Larochelle, H., Beygelzimer, A., d'Alché-Buc, F., Fox, E., and Garnett, R., editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc.
  * Xiong and Chen, (2021) Xiong, P. and Chen, Z. (2021).  Rsome in python: An open-source package for robust stochastic optimization made easy.
  * Ye et al., (1994) Ye, Y., Todd, M. J., and Mizuno, S. (1994).  An $o(\sqrt{n}l)$-iteration homogeneous and self-dual linear programming algorithm.  Mathematics of Operations Research, 19(1):53–67.
  * Zhen et al., (2021) Zhen, J., Kuhn, D., and Wiesemann, W. (2021).  A unified theory of robust and distributionally robust optimization via the primal-worst-equals-dual-best principle.

##  Appendix A Uncertainty sets and general reformulations

We derive convex reformulations of (1) for a general class of constraints $g$
and common uncertainty sets. We express $u$ as an affine transformation of an
uncertain parameter $z$, i.e., $u=Ax+b$, where $A\in{\mbox{\bf
R}}^{m\times\hat{m}}$, and $\hat{m}\leq m$, and $b\in{\mbox{\bf R}}^{m}$. The
proofs are in Appendix E. For notational simplicity, we suppress the
dependency of $g$ on $y$ here and for all following sections on
reformulations. We consider two general types of functions $g$: the maximum of
concave functions (left), where each $g_{l}$ is concave in $u$, and its
special case, the maximum of affine functions (right), where $P_{l}(x)$ and
$a_{l}(x)$ are convex functions in $x$:

$g(x,u)=\max_{l=1,\dots,L}g_{l}(x,u),\qquad
g(x,u)=\max_{l=1,\dots,L}(P_{l}(x)u+a_{l}(x)).$

##### Box and Ellipsoidal uncertainty.

The ellipsoidal uncertainty set is defined as

${\mathcal{U}}(\theta)=\left\\{u=Az+b\mid\|z\|_{p}\leq 1,\quad u\in
S\right\\},$

where $A\in{\mbox{\bf R}}^{m\times\hat{m}}$, $b\in{\mbox{\bf R}}^{m}$,
$p\in{\mbox{\bf R}}$ , and $S$ includes additional support information such as
upper and lower bounds. We let $p$ be the order of the norm, which is any
integer value $\geq 1$, or equal to $\infty$. In the latter case, this becomes
the box uncertainty set. The trainable parameter $\theta=(A,b)$. Note that if
$A=\rho I$, $b=0$, and $S={\mbox{\bf R}}^{m}$, this is equivalent to the
traditional ellipsoidal uncertainty set ${\\{u\mid\|u\|_{p}\leq\rho\\}}$. For
the maximum of concave constraint, we have the following reformulation.

$\begin{array}[]{ll}\text{minimize}&f(x)\\\ \mbox{subject
to}&[-g_{l}]^{\star}(x,h_{l})+\sigma_{S}(\omega_{l})-b^{T}(h_{l}+\omega_{l})+\|A^{T}(h_{l}+\omega_{l})\|_{q}\leq
0,\quad l=1,\dots,L,\end{array}$

with auxiliary variables $\omega_{l},h_{l}\in{\mbox{\bf R}}^{m}$ for
$l=1,\dots,L$, and $g^{\star}(x,h)$ the conjugate function. This result
follows from the theorem of infimal convolutions and composites for conjugate
functions (Rockafellar and Wets,, 1998), (Bauschke and Combettes,, 2017). As a
special case, we have the reformulation for the maximum of affine constraint,

$\begin{array}[]{ll}\mbox{minimize}&f(x)\\\ \mbox{subject
to}&a_{l}(x)+\sigma_{S}(\omega_{l})+b^{T}(P_{l}(x)-\omega_{l})+\|A^{T}(\omega_{l}-P_{l}(x))\|_{q}\leq
0,\quad l=1,\dots,L.\end{array}$

##### Budget uncertainty.

The budget uncertainty set requires both box and ellipsoidal constraints, and
is defined as

${\mathcal{U}}_{\text{bud}}(\theta)=\\{u=Az+b\mid\|z\|_{\infty}\leq\rho_{1},\quad\|z\|_{1}\leq\rho_{2},\quad
u\in S\\},$

where $\theta=(A,b,\rho_{1},\rho_{2})$. For the maximum of concave constraint,
we have the following reformulation.

$\begin{array}[]{ll}{\text{minimize}}&f(x)\\\ \mbox{subject
to}&[-g_{l}]^{\star}(x,h_{l})+\sigma_{S}(\omega_{l})-b^{T}(h_{l}+\omega_{l})+\rho_{1}\|r_{l}-A^{T}(h_{l}+\omega_{l})\|_{1}+\rho_{2}\|r_{l}\|_{\infty}\leq
0,\quad l=1,\dots,L,\end{array}$

with $r_{l}\in{\mbox{\bf R}}^{\hat{m}}$, $\omega_{l},h_{l}\in{\mbox{\bf
R}}^{m}$ for $l=1,\dots,L$ as auxiliary variables. The maximum of affine
reformulation is

$\begin{array}[]{ll}\mbox{minimize}&f(x)\\\ \mbox{subject
to}&a_{l}(x)+\sigma_{S}(\omega_{l})+b^{T}(P_{l}(x)-\omega)+\rho_{1}\|r_{l}-A^{T}(P_{l}(x)+\omega_{l})\|_{1}+\rho_{2}\|r_{l}\|_{\infty}\leq
0,\quad l=1,\dots,L.\end{array}$ (19)

##### Polyhedral uncertainty.

The polyhedral uncertainty set is defined as

${\mathcal{U}}(\theta)=\\{u=Az+b\mid Dz\leq d,\quad u\in S\\},$

where $\theta=(D,d)\in{\mbox{\bf R}}^{r\times\hat{m}}\times{\mbox{\bf R}}^{r}$
defines the polyhedron. For the maximum of concave constraint, we have the
following reformulation.

$\begin{array}[]{ll}{\text{minimize}}&f(x)\\\ \mbox{subject
to}&[-g_{l}]^{\star}(x,\omega_{l})-b^{T}(\omega_{l}+\gamma_{l})+\sigma_{S}(\gamma_{l})\leq
0,\quad l=1,\dots,L\\\ &A^{T}(\omega_{l}+\gamma_{l})=-D^{T}h_{l},\quad
l=1,\dots,L.\end{array}$

with auxiliary variables $\gamma_{l},\omega_{l}\in{\mbox{\bf R}}^{m}$,
$h_{l}\in{\mbox{\bf R}}^{r},$ for $l=1,\dots,L$. For the maximum of affine
constraint, we have:

$\begin{array}[]{ll}\mbox{minimize}&f(x)\\\ \mbox{subject
to}&a_{l}(x)+d^{T}h_{l}+b^{T}(P_{l}(x)-\gamma_{l})\leq 0,\quad l=1,\dots,L\\\
&A^{T}(P_{l}(x)-\gamma_{l})=D^{T}h_{l},\quad l=1,\dots,L\\\ &h_{l}\geq 0,\quad
l=1,\dots,L.\end{array}$

##### Mean robust uncertainty.

The Mean Robust Uncertainty (MRO) (Wang et al.,, 2022) set is defined as

$\displaystyle{\mathcal{U}}_{\text{mro}}(\theta)$
$\displaystyle=\left\\{u=A(v_{11},\dots,v_{LK})+b\in S^{K\times
L}\quad\middle|\quad\sum_{k=1}^{K}\sum_{l=1}^{L}\alpha_{lk}\|v_{lk}-\bar{d}_{k}\|^{p}\leq\rho^{p},\quad\alpha\in\Gamma\right\\},$

where $\Gamma=\\{\alpha~{}|~{}\sum_{l=1}^{L}\alpha_{lk}=w_{k},\alpha_{lk}\geq
0~{}\forall k,l\\}$. Here, we assume the constraint function $g$ of the RO to
be a maximum of $L$ pieces, indexed by $l$. We assume A is square, with
dimensions ${\mbox{\bf R}}^{m\times m}$, and we have $\theta=(A,b)$. This set
depends on $U^{N}$, a dataset of $N$ samples. We partition ${U}_{N}$ into $K$
disjoint subsets, with $\bar{d}_{k}$ the centroid of the $k$th subset. The
weight we place upon each subset is $w_{k}$, and is equivalent to the
proportion of points in the subset. We choose $p$ to be an integer exponent,
and $\rho$ to be a positive radius. For the maximum of concave constraint, we
have the following reformulation.

$\displaystyle\begin{array}[]{ll}{\text{minimize}}&f(x)\\\ \mbox{subject
to}&\lambda\rho^{p}+\sum_{k=1}^{K}w_{k}s_{k}\leq 0\\\
&[-g_{l}]^{\star}(x,h_{lk})-b^{T}(h_{lk}+\omega_{lk})+\sigma_{S}(\omega_{lk})-(A^{T}(h_{lk}+\omega_{lk}))^{T}\bar{d}_{k}\\\
&\qquad\qquad+\phi(q)\lambda\left\|A^{T}(h_{lk}+\omega_{lk})/\lambda\right\|^{q}_{*}\leq
s_{k},\quad l=1,\dots,L,\quad k=1,\dots,K\\\ &\lambda\geq 0.\end{array}$

with auxiliary variables $s,h_{lk},\omega_{lk}\in{\mbox{\bf R}}^{m}$ for
$j=1,\dots,J$, $k=1,\dots,K$, and $\lambda\in{\mbox{\bf R}}$. For the special
case of the maximum of affine constraint, we have

$\displaystyle\begin{array}[]{ll}\text{minimize }&f(x)\\\ \mbox{subject
to}&\lambda\rho^{p}+\sum_{k=1}^{K}w_{k}s_{k}\leq 0\\\
&a_{l}(x)+\sigma_{S}(\omega_{lk})-b^{T}(\omega_{lk}-P_{l}(x))-(A^{T}(\omega_{lk}-P_{l}(x)))^{T}\bar{d}_{k}\\\
&\qquad\qquad+\phi(r)\lambda\left\|A^{T}(\omega_{lk}-P_{l}(x))/\lambda\right\|^{r}_{q}\leq
s_{k},\quad l=1,\dots,L,\quad k=1,\dots,K\\\ &\lambda\geq 0.\end{array}$ (20)

###  A.1 Linear transformation to cone programs in standard form

The convex optimization problems above can be transformed to cone programs
(Agrawal et al., 2019a, ) in the form

$\begin{array}[]{ll}\mbox{minimize}&c^{T}x\\\ \mbox{subject to}&Gx+s=h\\\
&s\in\mathcal{K},\end{array}$ (21)

where $x\in{\mbox{\bf R}}^{\bar{n}}$ is the primal variable, $s\in{\mbox{\bf
R}}^{\bar{m}}$ is a slack variable, and the set
$\mathcal{K}\subseteq{\mbox{\bf R}}^{\bar{m}}$ is a nonempty, closed, convex
cone with dual cone $\mathcal{K}^{\star}\subseteq{\mbox{\bf R}}^{\bar{m}}$. So
long as the constituent functions $g_{l}$ of the constraint $g$ from the
robust problem (1) are concave in the uncertain parameters, and that these
functions satisfy the _(disciplined parametrized program) DPP_ rules proposed
by Agrawal et al., 2019a , the cone program’s data $(G,h,c)$ can be expressed
as a linear function (i.e., tensor product) of the parameters $\theta$ and $y$
(Agrawal et al., 2019a, , Lemma 1). The DPP rules depend on well-known
composition theorems for convex functions, and ensure that every function in
the program is convex in the decision variables, or affine if representing
equalities.

##  Appendix B Conservative Jacobians for the augmented Lagrangian

We give the exact forms of the conservative Jacobians from Proposition 3.1. In
this section, for simplicity, we denote the solution of the optimization
problem (1), $x^{j}=x(\theta,y^{j})$, and its conservative Jacobian,
$J_{x^{j}}(\theta)$, for each input pair $(\theta,y^{j})$. For the Lagrangian
function $\hat{L}(\alpha,s,\theta,\lambda,\mu)$ over the batch $B=(U\times
Y)\subseteq(U^{N}\times Y^{J})$, where $|B|$ denotes the cardinality of the
set, we have the following conservative Jacobians.

$\displaystyle J_{\hat{L}}(\theta)$ $\displaystyle=\frac{1}{|Y|}\sum_{y^{j}\in
Y}J_{f}(x^{j})J_{x^{j}}(\theta)$
$\displaystyle+\frac{\lambda}{\eta|B|}\sum_{(u^{i},y^{j})\in
B}\mathbf{1}\left\\{g(x^{j},u^{i},y^{j})\geq\alpha\right\\}J_{g}(x^{j})J_{x^{j}}(\theta)$
$\displaystyle+\frac{\mu}{\eta|B|}\sum_{(u^{i},y^{j})\in
B}((g(x^{j},u^{i},y^{j})-\alpha)/\eta+\alpha+s-\kappa)\mathbf{1}\left\\{g(x^{j},u^{i},y^{j})\geq\alpha\right\\}J_{g}(x^{j})J_{x^{j}}(\theta),$
$\displaystyle J_{\hat{L}}(s)$
$\displaystyle=\lambda+\frac{\mu}{|B|}\sum_{(u^{i},y^{j})\in
B}(g(x^{j},u^{i},y^{j})-\alpha)_{+}/\eta+\alpha+s-\kappa,$ $\displaystyle
J_{\hat{L}}(\alpha)$ $\displaystyle=\frac{\lambda}{|B|}\sum_{(u^{i},y^{j})\in
B}(1-(1/\eta)\mathbf{1}\left\\{g(x^{j},u^{i},y^{j})\geq\alpha\right\\})$
$\displaystyle+\frac{\mu}{|B|}\sum_{(u^{i},y^{j})\in
B}(1-(1/\eta)\mathbf{1}\left\\{g(x^{j},u^{i},y^{j})\geq\alpha\right\\})((g(x^{j},u^{i},y^{j})-\alpha)_{+}/\eta+\alpha+s-\kappa).$

For a condition $A$, we denote its indicator function $\mathbf{1}\\{A\\}$,
taking the value 1 when true and 0 when false. When the Lagrangian function is
defined over the entire dataset, we formulate the conservative Jabocian
$J_{L}(\alpha,s,\theta)$ the same way as above, with $B=(U^{N}\times Y^{J})$.

##  Appendix C Proof of conservative Jacobians for the robust optimization
problem

To prove Theorem 3.1, we characterize the path-differentiability of the
function $x(\theta,y)$ by solving the inner optimization problem and obtaining
conservative Jacobians using the results of Agrawal et al., 2019b  and Bolte
et al., (2021). The optimization problem (1) is first reformulated into a
convex and tractable problem in conic form (21), as shown in Appendix A.1. The
conic program’s data, $(G,h,c)$, is parameterized by $(\theta,y)$, through the
mapping

$(\theta,y)\rightarrow(G,h,c),$

which is linear in $\theta$ and $y$ for all uncertainty sets in Appendix A,
given the disciplined parametrized program (DPP) requirement detailed in
Appendix A.1 (Agrawal et al., 2019a, , Lemma 1).

This formulation has a corresponding dual, with dual variable $r\in{\mbox{\bf
R}}^{m}$, and solving the conic program amounts to solving the homogeneous
self-dual embedding (Ye et al.,, 1994) of the primal-dual problem. In the
following, we define $N=\bar{n}+\bar{m}+1$. The overall procedure is split
into three steps: mapping the problem data ($G,h,c$) onto a skew symmetric
matrix $Q$:

$Q=\left[\begin{array}[]{lll}0&G^{T}&c\\\ -G&0&h\\\
-c^{T}&-h^{T}&0\end{array}\right]\in{\mbox{\bf R}}^{N\times N},$

finding the solution $z=(u,v,w)\in{\mbox{\bf R}}^{N}$ to the homogeneous self-
dual embedding, and mapping $z$ onto the solution of the primal-dual pair
(Agrawal et al., 2019b, ). We can express the chain of mappings as

$(G,h,c)\rightarrow Q\rightarrow z\rightarrow(x,r,s).$

To solve for $z=\nu(Q)$, shown as the mapping $Q\rightarrow z$, we define the
normalized residual map (Busseti et al.,, 2019), given as

$\mathcal{N}(z,Q)=((Q-I)\Pi+I)(z/|w|),$

where $\Pi$ is the projection onto the convex cone
$\mathcal{K}^{\prime}={\mbox{\bf
R}}^{\bar{n}}\times\mathcal{K}^{\star}\times{\mbox{\bf R}}_{+}$. For a given
$Q$, $z=\nu(Q)$ is a solu tion to the homogeneous self-dual embedding if and
only if $\mathcal{N}(z,Q)=0$ and $w>0$. Once this solution is obtained, the
mapping $z\rightarrow(x,r,s)$ is given by the function

$\phi(z)=(u,\Pi(v),\Pi(v)-v)/w.$

We can now invoke proposition 4 of the work by Bolte et al., (2021), together
with the results from Agrawal et al., 2019b  to show that $\nu$ and $\phi$ are
path-differentiable.

######  Proposition C.1.

Assume the functions $\Pi$ and $\mathcal{N}$ are path differentiable, with
conservative Jacobians $J_{\Pi}$ and $J_{\mathcal{N}}$. Assume, for all $Q$
constructed from the problem data $(G,h,c)$, and $z=\nu(Q)$ is the unique
solution to $\mathcal{N}(z,Q)=0$. Furthermore, assume all matrices formed by
the first $N$ columns of $J_{\mathcal{N}}(z,Q)$ are invertible. Then, the
functions $\nu$ and $\phi$ are path-differentiable, with conservative
Jacobians

$J_{\nu}(Q)=\\{-U^{-1}V\mid[U~{}V]\in J_{\mathcal{N}}(z,Q)\\},\qquad
J_{\phi}(z)=\left[\begin{array}[]{ccc}I&0&-x\\\ 0&J_{\Pi}(v)&-r\\\
0&J_{\Pi}(v)-I&-s\end{array}\right].$

The assumptions on $\Pi$ and $\mathcal{N}$ hold, as the cones we consider
(zero cones, free cones, nonnegative cones, and second-order cones) are sub-
differentiable (Busseti et al.,, 2019; Kanzow et al.,, 2009). By Theorem 1 of
the work by Bolte et al., (2021), which states that the Clarke Jacobian is a
minimal conservative Jacobian, we conclude that the conic projections admit
conservative Jacobians, and are thus path differentiable. It follows that
$\mathcal{N}$, as a linear function of $\Pi$ and $Q$, is also path-
differentiable. The invertibility of $U$ holds as $U=((Q-I)J_{\Pi}(z)+I)/w$
(Agrawal et al., 2019b, ), which is nonsingular due to the PSD nature of the
members of the set $J_{\Pi}(z)$. Lastly, the assumption of uniqueness holds
from Assumption 3.3.

With this result, and by the linear nature of the mappings from $\theta$ to
$(G,h,c)$ and to $Q$, which preserves path-differentiablity, we can apply the
chain rule to obtain

$J_{\phi}(\theta)=J_{\phi}(v(Q))J_{v}(Q)J_{Q}(G,h,c)J_{(G,h,c)}(\theta).$

As $(x,r,s)=\phi(z)$, we can partition $J_{\phi}(\theta)$ as
$(J_{x}(\theta),J_{r}(\theta),J_{s}(\theta))$ where the first component is the
conservative Jacobian of interest. Therefore, we conclude that $x(\theta,y)$
is path differentiable, with nonempty conservative Jacobian $J_{x}(\theta)$.

##  Appendix D Proof of convergence of the augmented Lagrangian algorithm

We prove Theorem 3.2, the convergence of the augmented Lagrangian algorithm to
a stationary point. We separate the procedure into three steps. In the first
step, we show for Algorithm 2 with fixed $\lambda,\mu$, the convergence of
$\\{\alpha^{t},s^{t},\theta^{t}\\}$ to a stationary point of the augmented
Lagrangian function (11) as $t\to\infty$. In particular, we note that by
Theorem 1 of Bolte and Pauwels, (2019), conservative Jacobians are gradients
almost everywhere, so by Bolte et al., (2021), the iterates of the algorithm
is a Clarke stochastic subgradient sequence almost surely. We thus use the
Clarke differential (Clarke,, 1975), $\partial^{c}f(x)\coloneqq{\rm
conv}\\{\lim_{x_{i}\rightarrow x}\nabla f(x_{i})\mid x_{i}\in{\rm dom}(\nabla
f)\\}$ to denote the subdifferentials in our analysis, instead of using
$J_{f}(x)$. In the next two steps, we show for Algorithm 1, either the
convergence of $\\{\theta^{k}\\}$ to a feasible point of conditon (9), or the
convergence of $\\{\alpha^{k},s^{k},\theta^{k}\\}$ to a stationary point of
the constraint violation function.

##### Step 1: Convergence of ($\alpha$, $s$, $\theta$).

In this step, we consider the augmented Lagrangian function with fixed
$\lambda$ and $\mu$. By Assumptions 3.4 and 3.5 on the boundedness and step
sizes of the iterates, by the definability of the augmented Lagrangian
function, and by the existence of its conservative Jacobians, we apply Theorem
3 of Bolte et al., (2021), on the convergence of SGD for path-differentiable
functions. The projection onto the nonnegative orthant for $s$ does not affect
the convergence results, as proximal SGD with conservative Jacobians also
converges, subject to the same assumptions (Davis et al.,, 2018, Theorem 6.2).
Let $L(\alpha,s,\theta)=L(\alpha,s,\theta,\lambda,\mu)$, as the suppressed
variables are fixed, and let $(\alpha^{\star},s^{\star},\theta^{\star})$ be a
local minimum point of this function. By the aforementioned theorems, we can
show that $V_{\lambda,\mu}(\alpha,s,\theta)\coloneqq
L(\alpha,s,\theta)-L(\alpha^{\star},s^{\star},\theta^{\star})$ is a Lyapunov
function for the differential inclusion
$(\dot{\alpha},\dot{s},\dot{\theta})\in-\partial^{c}L(\alpha,s,\theta)$.
Specifically, there exist a ball of radius $r$ around
$(\alpha^{\star},s^{\star},\theta^{\star})$ such that for
$(\alpha,s,\theta)\in\mathbf{B}_{\alpha^{\star},s^{\star},\theta^{\star}}(r)$,
the Lyapunov function is locally positive definite. This function is
decreasing along solutions strictly outside of the set of Clarke-critical
points $\\{(\alpha,s,\theta)\in{\mbox{\bf
R}}^{p+2},0\in\partial^{c}L(\alpha,s,\theta)\\}$. Then, for almost all
initializations
$(\alpha^{0},s^{0},\theta^{0})\in\mathcal{D}_{\alpha}\times\mathcal{D}_{s}\times\mathcal{D}_{\theta}$,
the objective $L(\alpha^{t},s^{t},\theta^{t})$ converges almost surely and all
stationary points $(\bar{\alpha},\bar{s},\bar{\theta})$ of
$\\{\alpha^{t},s^{t},\theta^{t}\\}$ are Clarke-critical in the sense that
$0\in\partial^{c}L(\bar{\alpha},\bar{s},\bar{\theta})$ (Bolte et al.,, 2021).
We thus conclude that, with almost any initial point
$(\alpha^{0},s^{0},\theta^{0})\in\mathbf{B}_{\alpha^{\star},s^{\star},\theta^{\star}}(r)\subseteq\mathcal{D}_{\alpha}\times\mathcal{D}_{s}\times\mathcal{D}_{\theta}$,
the sequence $\\{\alpha^{t},s^{t},\theta^{t}\\}$ converges to a stationary
point $(\alpha^{\star},s^{\star},\theta^{\star})$, satisfying $0\in
J_{L}(\alpha^{\star},s^{\star},\theta^{\star})$.

##### Step 2: Convergence of $\lambda$.

From step 1, we have a stationary point $(\alpha^{k},s^{k},\theta^{k})$ of the
augmented Lagrangian function at iteration $k$ of Algorithm 1, with fixed
multiplier and penalty parameters $\lambda^{k-1}$ and $\mu^{k-1}$. Further
properties of this stationary point depend on the convergence, or lack
thereof, of the sequence $\\{\lambda^{k}\\}$. As the step size parameter of
$\lambda$ is the penalty $\mu$, we consider the two cases, (a) where the
sequence of penalty parameters $\\{\mu^{k}\\}$ is bounded, and (b) where
$\\{\mu^{k}\\}$ is unbounded. We first analyze case (a). By lines 6-12 of
Algorithm 1, we note that for $\\{\mu^{k}\\}$ to be bounded, the penalty
parameter must remain constant from some iteration onwards. By the condition
on line 6, with $\tau<1$, this implies
$\lim_{k\rightarrow\infty}\|H(\alpha^{k},\theta^{k})+s^{k}-\kappa\|_{2}=0$,
which in turn implies
$\lim_{k\rightarrow\infty}\mu^{k}(H(\alpha^{k},\theta^{k})+s^{k}-\kappa)=0$.
This sequence of updates for $\lambda$ is geometric and summable, so the
sequence $\\{\lambda^{k}\\}$ is bounded, with a convergent subsequence. Given
the boundedness assumption on the true Lagrangian multiplier $\lambda^{\star}$
in Assumption 3.4, we conclude that the whole sequence converges, ${\it
i.e.},$ $\\{\lambda^{k}\\}\rightarrow\lambda^{\star}$. We now turn to case
(b), where $\\{\mu^{k}\\}$ is unbounded. In this case, $\\{\lambda^{k}\\}$ may
not converge. However, by the $\lambda$ update rule in lines 7 of Algorithm 1,
we note that while $\\{\mu^{k}\\}\rightarrow\infty$, $\\{\lambda^{k}\\}$ is
bounded between $\lambda^{\rm{min}}$ and $\lambda^{\rm{max}}$, the safeguard
values.

##### Step 3: convergence of the outer loop.

We now show that either the sequence $\\{\theta^{k}\\}$ converges to a
feasible point of the $\widehat{\mathop{\bf CVaR}}$ condition, or
$\\{\alpha^{k},s^{k},\theta^{k}\\}$ converges to a stationary point of the
constraint violation function. We begin with the former. Suppose the sequence
of penalty parameters $\\{\mu^{k}\\}$ is bounded, and converges to
$\mu^{\star}$. Then, as $k\to\infty$, $\\{\lambda^{k}\\}$ converges to the
stationary point $\lambda^{\star}$. From step 1, we have that
$\\{\theta^{t}\\}$ in the inner loop converges almost surely to
$\theta^{k}=\theta^{\star}(\lambda^{k-1},\mu^{k-1})$ along each step of
$\lambda$ and $\mu$’s convergence. In order to establish the convergence of
$\\{\theta^{k}\\}$ in the outer loop to a stationary point
$\theta^{\star}(\lambda^{\star},\mu^{\star})$, a feasible point of the
$\widehat{\mathop{\bf CVaR}}$ condition (9), we must ensure that
$\theta^{\star}(\lambda,\mu)$ is well-behaved, by showing its continuity with
respect to $(\lambda,\mu)$.

Recall Assumption 3.3 on the uniqueness of the inner RO solutions
$x(\theta,y)$, Assumption 3.4 on the compactness of the set $\mathcal{X}$, and
the convexity (and therefore continuity) of the objective and constraint
functions $f$ and $g$ in $x$. Suppose these assumptions hold. By Theorem 4.2.1
of Bank et al., (1984), we establish the continuity of
$L(\alpha,s,\theta,\lambda,\mu)$ in $\theta$, for
$\theta\in\mathcal{D}_{\theta}$. As $L$ is bilinear and therefore continuous
in $(\lambda,\mu)$, we obtain our desired stability results, i.e., as
$\\{\lambda^{k}\\}$ and $\\{\mu^{k}\\}$ converge to $\lambda^{\star}$,
$\mu^{\star}$, the optimal solution with respect to $\theta$ is continuous,
and converges almost surely to the local minimizer $\theta^{\star}$. In
addition, by the strong convexity of $L$ in $s$, the optimal solution with
respect to $s$ also converges to the global minimizer $s^{\star}$. Note that
$\\{\alpha^{k}\\}$ does not necessarily converge to a stationary point, as the
minimum over $\alpha$ is not necessarily unique (Uryasev and Rockafellar,,
2001, Theorem 2). However, by the convergence of the inner loop in step 1 and
by the convergence condition of $\mu$ in step 2, as $k\to\infty$ in Algorithm
1, we have a point
$(\alpha^{k},s^{\star},\theta^{\star},\lambda^{\star},\mu^{\star})$ satisfying
$0\in J_{L}(\alpha^{k},s^{\star},\theta^{\star})$ and
$H(\alpha^{k},\theta^{\star})+s^{\star}-\kappa=0$. We can conclude that
$\\{\theta^{k}\\}$ converges almost surely to a feasible point
$\theta^{\star}$ of the $\widehat{\mathop{\bf CVaR}}$ condition, where the
minimizer over $\alpha$ takes some value $\alpha^{\star}\leq\alpha^{k}$.

On the other hand, suppose the sequence of penalty parameters $\\{\mu^{k}\\}$
is unbounded. By step 1, for any $\lambda^{k}$ and $\mu^{k}$, we have obtained
a stationary point of the augmented Lagrangian function, satisfying $0\in
J_{L}(\alpha^{k},s^{k},\theta^{k})$. Note that, with $\\{\mu^{k}\\}$ unbounded
and $\\{\lambda^{k}\\}$ bounded, we observe
$\lim_{k\rightarrow\infty}\lambda^{k}/\mu^{k}=0$. Therefore, dividing both
sides of the conservative Jacobian by $\mu^{k}$ and taking the limit, we
obtain:

$\displaystyle 0$ $\displaystyle\in\lim_{k\rightarrow\infty}\frac{1}{\eta
NJ}\sum_{y^{j}\in
Y^{J}}\mathbf{1}\left\\{g(x^{j},u^{i},y^{j})\geq\alpha^{k}\right\\}((g(x^{j},u^{i},y^{j})-\alpha^{k})_{+}/\eta+\alpha^{k}+s^{k}-\kappa)J_{g}(x^{j})J_{x^{j}}(\theta^{k}),$
$\displaystyle 0$
$\displaystyle\in\lim_{k\rightarrow\infty}\frac{1}{NJ}\sum_{(u^{i},y^{j})\in(U^{N}\times
Y^{J})}(g(x^{j},u^{i},y^{j})-\alpha^{k})_{+}/\eta+\alpha^{k}+s^{k}-\kappa,$
$\displaystyle 0$
$\displaystyle\in\lim_{k\rightarrow\infty}\frac{1}{NJ}\sum_{(u^{i},y^{j})\in(U^{N}\times
Y^{J})}\left(1-(1/\eta)\mathbf{1}\left\\{g(x^{j},u^{i},y^{j})\geq\alpha^{k}\right\\}\right)((g(x^{j},u^{i},y^{j})-\alpha^{k})_{+}/\eta+\alpha^{k}+s^{k}-\kappa).$

These are the stationary conditions for the constraint violation function
$V(\alpha,s,\theta)=(1/2)\|H(\alpha,\theta)+s-\kappa\|^{2}$, i.e., $0\in
J_{V}(\alpha^{k},s^{k},\theta^{k})$. We thus converge to a stationary point of
the constraint violation function, as desired.

##  Appendix E General reformulation proofs

###  E.1 Conjugate of a general function

We first derive the conjugate of general functions, which can then be used in
the upcoming reformulations.

##### Linear function.

For a linear function $g(x,u)=(Pu+a)^{T}x$, we compute the conjugate as

$[-g]^{\star}(x,h)=\sup_{u}h^{T}u-(a+Pu)^{T}x=\begin{dcases}a^{T}x&\text{
if}\;h=-P^{T}x\\\ \infty&\text{ otherwise}.\end{dcases}$

##### Sum of concave functions.

For a general function $g(x,u)=\sum_{j=1}^{J}g_{j}(u,x)$, where each $g_{j}$
is concave in $u$, by the theorem of infimal convolutions (Rockafellar and
Wets,, 1998, Theorem 11.23(a), p. 493), we compute the conjugate as

$\displaystyle[-g]^{\star}(x,h)=\begin{cases}\underset{h_{j},\dots,h_{J}}{\inf}&\sum_{j=1}^{J}[-g_{j}]^{\star}(x,h_{j})\\\
\mbox{subject to}&\sum_{j=1}^{J}h_{j}=h.\end{cases}$

##### Maximum of concave functions.

We also consider the maximum of concave functions,

$g(x,u)=\max_{l=1,\dots,L}g_{l}(x,u).$

In practice, these functions $g(x,u)$ are often affine functions. For all the
reformulations below, we assume $g(x,u)$ to be of this maximum of concave
form.

###  E.2 Box and Ellipsoidal reformulations

Rewriting the robust constraint as optimization problem, we have

$\displaystyle\begin{array}[t]{ll}\underset{Az+b\in S}{\sup}&g(x,Az+b)\\\
\mbox{subject to}&\|z\|_{p}\leq 1\end{array}$
$\displaystyle=\begin{array}[t]{ll}\underset{Az+b\in
S}{\sup}~{}\underset{\lambda\geq
0}{\inf}~{}g(x,Az+b)-\lambda\|z\|_{p}+\lambda\\\ \end{array}$
$\displaystyle=\begin{array}[t]{ll}\underset{Az+b\in
S}{\sup}~{}\underset{\lambda\geq
0}{\inf}~{}\underset{l}{\max}~{}g_{l}(x,Az+b)-\lambda\|z\|_{p}+\lambda\\\
\end{array}$
$\displaystyle=\begin{array}[t]{ll}\underset{l}{\max}~{}\underset{Az+b\in
S}{\sup}~{}\underset{\lambda\geq
0}{\inf}~{}g_{l}(x,Az+b)-\lambda\|z\|_{p}+\lambda\\\ \end{array}$
$\displaystyle=\begin{array}[t]{ll}\underset{l}{\max}~{}\underset{\lambda\geq
0}{\inf}~{}\underset{Az+b\in
S}{\sup}~{}g_{l}(x,Az+b)-\lambda\|z\|_{p}+\lambda\\\ \end{array}$
$\displaystyle=\begin{array}[t]{ll}\underset{\lambda_{l}\geq
0,\tau}{\inf}&\tau\\\ \mbox{subject to}&\lambda_{l}+\underset{Az+b\in
S}{\sup}g_{l}(x,Az+b)-\lambda_{l}\|z\|_{p}\leq\tau,\quad
l=1,\dots,L,\end{array}$

where $p=\infty$ for box uncertainty and $p\geq 1$ otherwise. We used the
Lagrangian in the first equality, and decomposed $g(x,Az+b)$ into its
constituents for the second equality. As $l$ and $\lambda$ are separable, in
the third step we interchanged the $\inf$ and the $\max$. Next, for each $l$,
as the expression is convex in $\lambda$ and concave in $u$, we used the Von
Neumann-Fan minimax theorem to interchange the $\inf$ and the $\sup$. We then
express the $\max$ in epigraph form by introducing an auxiliary variable
$\tau$.

Now, if we define the function $c(z)=\|z\|_{p}$, then using the definition of
conjugate functions, we have, for each $l$, the left-hand-side of the
constraint equal to

$\displaystyle\lambda_{l}+[-g_{l}(x,A\cdot+b)+\chi_{S}(A\cdot+b)+\lambda_{l}c]^{\star}(0).$

We then borrow results from (Rockafellar and Wets,, 1998, Theorem 11.23(a), p.
493), (Zhen et al.,, 2021, Lemma B.8), and (Bauschke and Combettes,, 2017,
Proposition 13.24) with regards to the inf-convolution of conjugate functions,
the conjugate functions of norms, and the conjugate of linear compositions, to
note

$\displaystyle[-g_{l}(x,A\cdot+b)+\chi_{S}(A\cdot+b)+\lambda_{l}c]^{\star}(0)=\underset{\omega_{l},h_{l},y_{l}}{\inf}([-g_{l}]^{\star}(x,h_{l})$
$\displaystyle-b^{T}(h_{l}+\omega_{l})+\sigma_{S}(\omega_{l})+\mathcal{I}_{\|y_{l}\|_{q}\leq\lambda_{l}}$
$\displaystyle+\mathcal{I}_{A^{T}(h_{l}+\omega_{l})=y_{l}}),$

where $h_{l}\in{\mbox{\bf R}}^{m}$, $y_{l}\in{\mbox{\bf R}}^{\hat{m}}$,
$\omega_{l}\in{\mbox{\bf R}}^{m}$, and $q$ is the conjugate number of $p$.
Substituting these back, we arrive at

$\displaystyle\begin{array}[t]{ll}\underset{x,\tau,h,\omega,y}{\inf}&\tau\\\
\mbox{subject
to}&[-g_{l}]^{\star}(x,h_{l})+\sigma_{S}(\omega_{l})-b^{T}(h_{l}+\omega_{l})+\|y_{l}\|_{q}\leq\tau,\quad
l=1,\dots,L\\\ &A^{T}(h_{l}+\omega_{l})=y_{l},\quad l=1,\dots,L.\end{array}$

Recall that this optimization problem upper bounds our original robust
constraint, which we constrain to be nonpositive. Therefore, we only need
$\tau\leq 0$. For the final step, we can now write the original minimization
problem (1) in minimization form for both the objective and constraints, to
arrive at

$\displaystyle\begin{array}[t]{ll}\underset{x,h,\omega}{\text{minimize}}&f(x)\\\
\mbox{subject
to}&[-g_{l}]^{\star}(x,h_{l})+\sigma_{S}(\omega_{l})-b^{T}(h_{l}+\omega_{l})+\|A^{T}(h_{l}+\omega_{l})\|_{q}\leq
0,\quad l=1,\dots,L.\end{array}$

###  E.3 Budget reformulation

Rewriting the robust constraint as optimization problem, we have

$\displaystyle\begin{array}[t]{ll}\underset{Az+b\in S}{\sup}&g(x,Az+b)\\\
\mbox{subject to}&\|z\|_{\infty}\leq\rho_{1}\\\ &\|z\|_{1}\leq\rho_{2}\\\
\end{array}$ $\displaystyle=\begin{array}[t]{ll}\underset{Az+b\in
S}{\sup}~{}\underset{\lambda\geq
0}{\inf}~{}\underset{l}{\max}~{}g_{l}(x,Az+b)-\lambda_{1}\|z\|_{\infty}+\lambda_{1}\rho_{1}-\lambda_{2}\|z\|_{1}+\lambda_{2}\rho_{2}\\\
\end{array}$
$\displaystyle=\begin{array}[t]{ll}\underset{l}{\max}~{}\underset{Az+b\in
S}{\sup}~{}\underset{\lambda\geq
0}{\inf}~{}g_{l}(x,Az+b)-\lambda_{1}\|z\|_{\infty}+\lambda_{1}\rho_{1}-\lambda_{2}\|z\|_{1}+\lambda_{2}\rho_{2}\\\
\end{array}$
$\displaystyle=\begin{array}[t]{ll}\underset{l}{\max}~{}\underset{\lambda\geq
0}{\inf}~{}\underset{Az+b\in
S}{\sup}~{}g_{l}(x,z)-\lambda_{1}\|z\|_{\infty}+\lambda_{1}\rho_{1}-\lambda_{2}\|z\|_{1}+\lambda_{2}\rho_{2},\\\
\end{array}$

where we again used the Lagrangian, separated $g$ into its constituents, and
applied the minimax theorem. We now express the max in epigraph form,

$\displaystyle\begin{array}[t]{ll}\underset{\lambda_{l}\geq
0,\tau}{\inf}&\tau\\\ \mbox{subject
to}&\lambda_{l1}\rho_{1}+\lambda_{l2}\rho_{2}+\underset{Az+b\in
S}{\sup}g_{l}(x,Az+b)-\lambda_{l1}\|z\|_{\infty}-\lambda_{l2}\|z\|_{1}\leq\tau,\quad
l=1,\dots,L.\end{array}$

If we define the functions $c_{1}(z)=\|z\|_{\infty}$, $c_{2}(z)=\|z\|_{1}$,
then using the definition of conjugate functions, we have for each $l$ the
left-hand-side of the constraint equal to

$\displaystyle\lambda_{l1}\rho_{1}+\lambda_{l2}\rho_{2}+[-g_{l}(x,A\cdot+b)+\chi_{S}(A\cdot+b)+\lambda_{l1}c_{1}+\lambda_{l2}c_{2}]^{\star}(0)$
$\displaystyle=\underset{h_{l},\omega_{l},y_{l}}{\inf}\lambda_{l1}\rho_{1}+\lambda_{l2}\rho_{2}+[-g_{l}]^{\star}(x,h_{l})-b^{T}(h_{l}+\omega_{l})+\sigma_{S}(\omega_{l})+[\lambda_{l1}c_{1}+\lambda_{l2}c_{2}]^{\star}(-y_{l})$
$\displaystyle\hskip 341.43306pt+\mathcal{I}_{A^{T}(h_{l}+\omega_{l})=y_{l}}$
$\displaystyle=\underset{h_{l},\omega_{l},y_{l},r_{l}}{\inf}\lambda_{l1}\rho_{1}+\lambda_{l2}\rho_{2}+[-g_{l}]^{\star}(x,h_{l})-b^{T}(h_{l}+\omega_{l})+\sigma_{S}(\omega_{l})+\mathcal{I}_{\|r_{l}-y_{l}\|_{1}\leq\lambda_{l1}}+\mathcal{I}_{\|r_{l}\|_{\infty}\leq\lambda_{l2}}$
$\displaystyle\hskip 341.43306pt+\mathcal{I}_{A^{T}(h_{l}+\omega_{l})=y_{l}}$

We can express this as

$\displaystyle\begin{array}[t]{ll}\underset{x,h_{l},\omega_{l},y_{l},r_{l},\tau}{\inf}&\tau\\\
\mbox{subject
to}&[-g_{l}]^{\star}(x,h_{l})+\sigma_{S}(\omega_{l})-b^{T}(h_{l}+\omega_{l})+\rho_{1}\|r_{l}-y_{l}\|_{1}+\rho_{2}\|r_{l}\|_{\infty}\leq\tau,\\\
&\hfill k=1,\dots,K\\\ &A^{T}(h_{l}+\omega_{l})=y_{l},\quad
l=1,\dots,L.\end{array}$

We can then let $\tau=0$ and substitute the new constraints in for the robust
constraint in (1),

$\begin{array}[t]{ll}{\text{minimize}}&f(x)\\\ \mbox{subject
to}&[-g_{l}]^{\star}(x,h_{l})+\sigma_{S}(\omega_{l})-b^{T}(h_{l}+\omega_{l})+\rho_{1}\|r_{l}-A^{T}(h_{l}+\omega_{l})\|_{1}+\rho_{2}\|r_{l}\|_{\infty}\leq
0,\\\ &\hfill l=1,\dots,L.\end{array}$

###  E.4 Polyhedral reformulation

Rewriting the robust constraint as an optimization problem, we have

$\displaystyle\begin{array}[t]{ll}\underset{Az+b\in S}{\sup}&g(x,Az+b)\\\
\mbox{subject to}&Dz\leq d\end{array}$
$\displaystyle=\begin{array}[t]{ll}\underset{Az+b\in
S}{\sup}~{}\underset{h\geq 0}{\inf}~{}g(x,Az+b)-(Dz)^{T}h+d^{T}h\\\
\end{array}$ $\displaystyle=\begin{array}[t]{ll}\underset{Az+b\in
S}{\sup}~{}\underset{h\geq
0}{\inf}~{}\underset{l}{\max}~{}g_{l}(x,Az+b)-(Dz)^{T}h+d^{T}h\\\ \end{array}$
$\displaystyle=\begin{array}[t]{ll}\underset{l}{\max}~{}\underset{Az+b\in
S}{\sup}~{}\underset{h\geq 0}{\inf}~{}g_{l}(x,Az+b)-(Dz)^{T}h+d^{T}h\\\
\end{array}$
$\displaystyle=\begin{array}[t]{ll}\underset{l}{\max}~{}\underset{h\geq
0}{\inf}~{}\underset{Az+b\in S}{\sup}~{}g_{l}(x,Az+b)-(Dz)^{T}h+d^{T}h\\\
\end{array}$ $\displaystyle=\begin{array}[t]{ll}\underset{h\geq
0,\tau}{\inf}&\tau\\\ \mbox{subject to}&d^{T}h_{l}+\underset{Az+b\in
S}{\sup}g_{l}(x,Az+b)-(Dz)^{T}h_{l}\leq\tau,\quad l=1,\dots,L\end{array}$
$\displaystyle=\begin{array}[t]{ll}\underset{h\geq
0,\omega,\tau,y,\gamma}{\inf}&\tau\\\ \mbox{subject
to}&d^{T}h_{l}+[-g_{l}]^{\star}(x,\omega_{l})-b^{T}(\omega_{l}+\gamma_{l})+\sigma_{S}(\gamma_{l})+\mathcal{I}_{D^{T}h_{l}=-y_{l}}\\\
&\hfill+\mathcal{I}_{A^{T}(\omega_{l}+\gamma_{l})=y_{l}}\leq\tau,\quad
l=1,\dots,L\end{array}$ $\displaystyle=\begin{array}[t]{ll}\underset{z\geq
0,\omega,\tau,\gamma}{\inf}&\tau\\\ \mbox{subject
to}&[-g_{l}]^{\star}(x,\omega_{l})-b^{T}(\omega_{l}+\gamma_{l})+\sigma_{S}(\gamma_{l})\leq\tau,\quad
l=1,\dots,L\\\ &A^{T}(\omega_{l}+\gamma_{l})=-D^{T}h_{l},\quad
l=1,\dots,L,\end{array}$

where we again used the Lagrangian, the Von Neumann-Fan minimax theorem, the
infimal convolution of conjugates, and the conjugate of affine compositions.
The final reformulation (1) is then

$\begin{array}[t]{ll}{\text{minimize}}&f(x)\\\ \mbox{subject
to}&[-g_{l}]^{\star}(x,\omega_{l})-b^{T}(\omega_{l}+\gamma_{l})+\sigma_{S}(\gamma_{l})\leq
0,\quad l=1,\dots,L\\\ &A^{T}(\omega_{l}+\gamma_{l})=-D^{T}h_{l},\quad
l=1,\dots,L.\end{array}$

###  E.5 Mean Robust Optimization reformulation

We adapt the proof from (Wang et al.,, 2022, Appendix A.3). To simplify
notation, we define
$c_{k}(v_{lk})\coloneqq\|v_{lk}-\bar{d}_{k}\|^{p}-\epsilon^{p}$. We also note
that $u_{lk}=Av_{lk}+b$. Then, we have

$\displaystyle\begin{array}[t]{ll}\underset{u_{11},\dots,u_{LK}\in
S,\alpha\in\Gamma}{\sup}&\sum_{k=1}^{K}\sum_{l=1}^{L}\alpha_{lk}g_{l}(x,u_{lk})\\\
\mbox{~{}~{}~{}~{}subject
to}&\sum_{k=1}^{K}\sum_{l=1}^{L}\alpha_{lk}c_{k}(v_{lk})\leq 0\end{array}$
$\displaystyle=\begin{array}[t]{ll}\underset{u_{11},\dots,u_{LK}\in
S,\alpha\in\Gamma}{\sup}&\underset{\lambda\geq
0}{\inf}~{}\sum_{k=1}^{K}\sum_{l=1}^{L}\alpha_{lk}g_{l}(x,u_{lk})-\lambda\sum_{k=1}^{K}\sum_{l=1}^{L}\alpha_{lk}c_{k}(v_{lk})\\\
\end{array}$
$\displaystyle=\begin{array}[t]{ll}\underset{\alpha\in\Gamma}{\sup}~{}\underset{\lambda\geq
0}{\inf}~{}\underset{u_{11},\dots,u_{LK}\in
S}{\sup}&\sum_{k=1}^{K}\sum_{l=1}^{L}\alpha_{lk}g_{l}(x,u_{lk})-\lambda\sum_{k=1}^{K}\sum_{l=1}^{L}\alpha_{lk}c_{k}(v_{lk}).\\\
\end{array}$

We applied the Lagrangian in the first equality. Then, as the summation is
over upper-semicontinuous functions $g_{l}(x,u_{lk})$ concave in $u_{lk}$, we
applied the Von Neumann-Fan minimax theorem to interchange the inf and the
sup. Next, we rewrite the formulation using an epigraph trick, and make a
change of variables.

$\displaystyle=\begin{array}[t]{ll}\underset{\alpha\in\Gamma}{\sup}~{}\underset{\lambda\geq
0}{\inf}&\sum_{k=1}^{K}s_{k}\\\ \mbox{subject
to}~{}\underset{u_{11},\dots,u_{LK}\in
S}{\sup}&\sum_{l=1}^{L}\alpha_{lk}(g_{l}(x,u_{lk})-\lambda c_{k}(v_{lk}))\leq
s_{k},\quad k=1,\dots,K\end{array}$
$\displaystyle=\begin{array}[t]{ll}\underset{\alpha\in\Gamma}{\sup}~{}\underset{\lambda\geq
0}{\inf}&\sum_{k=1}^{K}s_{k}\\\ \mbox{subject
to}~{}\underset{\alpha_{11}u_{11},\dots,\alpha_{LK}u_{LK}\in
S^{\prime}}{\sup}&\sum_{l=1}^{L}\alpha_{lk}(g_{l}(x,(\alpha_{lk}u_{lk})/\alpha_{lk})\\\
&-\lambda c_{k}((\alpha_{lk}v_{lk})/\alpha_{lk}))\leq s_{k},\quad
k=1,\dots,K.\end{array}$

In the last step, we rewrote $u_{lk}=(\alpha_{lk}u_{lk})/\alpha_{lk}$,
$v_{lk}=(\alpha_{lk}v_{lk})/\alpha_{lk}$, and maximized over
$\alpha_{lk}u_{lk}\in S^{\prime}$, where $S^{\prime}$ is the support
transformed by $\alpha_{lk}$. See Wang et al., (2022) for details. Next, we
make substitutions $h_{ij}=\alpha_{lk}u_{lk}$, $z_{ij}=\alpha_{lk}v_{lk}$, and
define functions
$g_{l}^{\prime}(x,h_{lk})=\alpha_{lk}g_{l}(x,h_{lk}/\alpha_{lk})$,
$c_{k}^{\prime}(z_{lk})=\alpha_{lk}c_{k}^{\prime}(z_{lk}/\alpha_{lk})$.

$\displaystyle=\begin{array}[t]{ll}\underset{\alpha\in\Gamma}{\sup}~{}\underset{\lambda\geq
0}{\inf}&\sum_{k=1}^{K}s_{k}\\\ \mbox{subject
to}~{}\underset{h_{11},\dots,h_{lk}\in
S^{\prime}}{\sup}&\sum_{l=1}^{L}g_{l}^{\prime}(x,h_{lk})-\lambda
c_{k}^{\prime}(z_{lk})\leq s_{k},\quad k=1,\dots,K\end{array}$
$\displaystyle=\begin{array}[t]{ll}\underset{\alpha\in\Gamma}{\sup}~{}\underset{\lambda\geq
0}{\inf}&\sum_{k=1}^{K}s_{k}\\\ \mbox{subject
to}&\sum_{l=1}^{L}[-g_{l}^{\prime}(x,A\cdot+b)+\chi_{S^{\prime}}(A\cdot+b)+\lambda
c_{k}^{\prime}]^{\star}(0)\leq s_{k},\quad k=1,\dots,K.\end{array}$

For the new functions defined, we applied the definition of conjugate
functions. Now, using the conjugate form $f^{\star}(y)=\alpha g^{\star}(y)$ of
a right-scalar-multiplied function $f(x)=\alpha g(x/\alpha)$, we rewrite the
above as

$\displaystyle=\begin{array}[t]{ll}\underset{\alpha\in\Gamma}{\sup}~{}\underset{\lambda\geq
0}{\inf}&\sum_{k=1}^{K}s_{k}\\\ \mbox{subject
to}&\sum_{l=1}^{L}\alpha_{lk}[-g_{l}(x,A\cdot+b)+\chi_{S}(A\cdot+b)+\lambda
c_{k}]^{\star}(0)\leq s_{k},\quad k=1,\dots,K.\end{array}$

We note that while the support of $h_{lk}$ is $S^{\prime}$, the support of
$h_{lk}/\alpha_{lk}$, which is the input of $g_{l}$, is still $S$. Now, again
using properties of conjugate functions, we note that:

$\displaystyle\alpha_{lk}[(-g_{l}(x,A\cdot+b)+\chi_{S}(A\cdot+b)+\lambda
c_{k})]^{\star}(0)$
$\displaystyle=\alpha_{lk}\underset{\omega_{lk},h_{lk},y_{lk}}{\inf}([-g_{l}]^{\star}(x,h_{lk})-b^{T}(h_{lk}+\omega_{lk})$
$\displaystyle+\mathcal{I}_{A^{T}(h_{lk}+\omega_{lk})=y_{lk}}+\sigma_{S}(\omega_{lk})+[\lambda
c_{k}]^{\star}(-y_{lk}))$

and

$\displaystyle{}^{\star}(-y_{lk})=-y_{lk}^{T}\bar{d}_{k}+\phi(q)\lambda{\|}y_{lk}/\lambda{\|}^{q}_{*}+\lambda\epsilon^{p}.$

Substituting these back, we arrive at

$\displaystyle\begin{array}[t]{ll}\underset{\alpha\in\Gamma}{\sup}~{}\underset{\lambda\geq
0,h_{lk},\omega_{lk},s_{k}}{\inf}&\sum_{k}^{K}s_{k}\\\ \mbox{subject
to}&\sum_{l=1}^{L}\alpha_{lk}([-g_{l}]^{\star}(x,h_{lk})-b^{T}(h_{lk}+\omega_{lk})\\\
&+\sigma_{S}(\omega_{lk})-(A^{T}(h_{lk}+\omega_{lk}))^{T}\bar{d}_{k}+\phi(q)\lambda{\|}A^{T}(h_{lk}+\omega_{lk})/\lambda{\|}^{q}_{*}+\lambda\epsilon^{p})\leq
s_{k},\\\ &\hskip 56.9055pt\quad k=1,\dots,K,\quad j=1,\dots,J.\end{array}$

Taking the supremum over $\alpha$, noting that
$\sum_{l=1}^{L}\alpha_{lk}=w_{k}$ for all $k$, we arrive at

$\displaystyle\begin{array}[t]{ll}\underset{\lambda\geq
0,z_{lk},\omega_{lk},s_{k}}{\inf}&\lambda\epsilon^{p}+\sum_{k}^{K}w_{k}s_{k}\\\
\mbox{s.t.}&[-g_{l}]^{\star}(x,h_{lk})-b^{T}(h_{lk}+\omega_{lk})+\sigma_{S}(\omega_{lk})-(A^{T}(h_{lk}+\omega_{lk}))^{T}\bar{d}_{k}\\\
&\hfill+\phi(q)\lambda\left\|A^{T}(h_{lk}+\omega_{lk})/\lambda\right\|^{q}_{*}\leq
s_{k},\\\ &\hfill l=1,\dots,L,\quad k=1,\dots,K,\end{array}$

for which the final reformulation of (1) is

$\displaystyle\begin{array}[t]{ll}{\text{minimize}}&f(x)\\\ \mbox{subject
to}&\lambda\epsilon^{p}+\sum_{k}^{K}w_{k}s_{k}\leq 0\\\
&[-g_{l}]^{\star}(x,h_{lk})-b^{T}(h_{lk}+\omega_{lk})+\sigma_{S}(\omega_{lk})-(A^{T}(h_{lk}+\omega_{lk}))^{T}\bar{d}_{k}\\\
&\hfill+\phi(q)\lambda\left\|A^{T}(h_{lk}+\omega_{lk})/\lambda\right\|^{q}_{*}\leq
s_{k},\\\ &\hfill l=1,\dots,L,\quad k=1,\dots,K.\end{array}$