# Adjustable robust optimization with objective uncertainty

Boris Detienne, Henri Lefebvre, Enrico Malaguti, Michele Monaci

► To cite this version:

Boris Detienne, Henri Lefebvre, Enrico Malaguti, Michele Monaci. Adjustable robust optimization with objective uncertainty. European Journal of Operational Research, 2024, 312 (1), pp.373-384. 10.1016/j.ejor.2023.06.042. hal-04235571

HAL Id: hal-04235571

https://hal.science/hal-04235571v1

Submitted on 2 Sep 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

ARTICLE INFO

Article history:

Received 15 July 2022

Accepted 30 June 2023

Available online 4 July 2023

Keywords:

Uncertainty modelling

Two-stage robust optimization

Reformulation

Fenchel duality

Branch-and-price

Computational experiments

A B S T R A C T

In this work, we study optimization problems where some cost parameters are not known at decision time and the decision flow is modeled as a two-stage process within a robust optimization setting. We address general problems in which all constraints (including those linking the first and the second stages) are defined by convex functions and involve mixed-integer variables, thus extending the existing literature to a much wider class of problems. We show how these problems can be reformulated using Fenchel duality, allowing to derive an enumerative exact algorithm, for which we prove asymptotic convergence in the general case, and finite convergence for cases where the first-stage variables are all integer.

An implementation of the resulting algorithm, embedding a column generation scheme, is then computationally evaluated on a variant of the Capacitated Facility Location Problem with uncertain transportation costs, using instances that are derived from the existing literature. To the best of our knowledge, this is the first approach providing results on the practical solution of this class of problems.

© 2023 The Authors. Published by Elsevier B.V.

This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/)

1. Introduction

Robust Optimization (RO) has emerged as a solution approach to deal with uncertainty in optimization problems. Contrary to stochastic optimization, another popular approach, robust optimization does not rely on probability distributions. Indeed, RO considers an uncertainty set for the unknown parameters, against which the taken decision should be immune. In that sense, constraints have to be respected in every possible realization of the parameters and the objective function evaluated in the least advantageous case. The concept was first introduced in Soyster (1973) and received considerable attention in the scientific literature. Recent advances in RO can be found in Bertsimas, Brown, & Caramanis (2010), Hassene, Bazgan, & Vanderpooten (2009), Ben-Tal, Ghaoui, & Nemirovski (2009), Leyffer, Menickelly, Munson, Vanaret, & Wild (2020) and Yanıkoğlu, Gorissen, & den Hertog (2019), among others.

More formally, a basic (one-stage) robust optimization problem can be cast as follows

$$ \inf_{z} \sup_{\xi \in \Xi} f(\xi, z) $$

$$ \text{subject to} \quad g(\xi, z) \le 0 \quad \forall \xi \in \Xi \\ z \in Z. \qquad (1SR-P) $$

Here, the unknown data is represented by variables $\xi$ that belong to the so-called uncertainty set $\Xi$. As mentioned above, decision $z$ has to be feasible in every possible occurrence of the uncertainty, hence robust solutions tend to be overly conservative. To tackle this drawback, Ben-Tal, Goryashko, & Guslitzer (2004) introduced the so-called adjustable robust optimization, also known as two-stage robust optimization. As its name suggests, in a two-stage context, part of the decisions are made in a first stage (i.e., herein-and-now, before uncertainty reveals), while recourse decisions can be taken in a second stage (i.e., once the actual values of the uncertain data are known) as an attempt to react to the outcome of the uncertain process. Typically, the feasible region of (1SR-P) can, indeed, be recast to embed a two-stage decision process by splitting variables $z$ in $(x, y)$. Here, $x \in \mathcal{X}$ are decisions to be made here and now, while $y \in \mathcal{Y}$ may be taken at a later instant. Accordingly, set $Z$ is defined as $\mathcal{X} \times \mathcal{Y}$. With the convention that the minimum objective function value for an infeasible problem is $+\infty$, a two-stage robust problem can be formulated as follows

$$ \inf_{x \in \mathcal{X}} \inf_{\xi \in \Xi} f(\xi, x, y), \qquad (2SR-P) $$

where $\mathcal{Y}(x, \xi) = \{y : y \in \mathcal{Y}, g(\xi, x, y) \le 0\}$, and $g(\xi, x, y) \le 0$ are the so-called linking constraints. Set $\mathcal{X}$ is now referred to as the

$$ \inf_{x \in \mathcal{X}} \inf_{\xi \in \Xi} f(\xi, x, y), \qquad (2SR-P) $$

<footnote>ªIMB UMR CNRS 5251, Inria Bordeaux Sud-Ouest, Université de Bordeaux, 200 Avenue de la Vieille Tour, Talence, 33405, France</footnote>

<footnote>ᵇDipartimento di Ingegneria dell'Energia Elettrica e dell'Informazione "Guglielmo Marconi", Universit di Bologna, Viale del Risorgimento, 2, Bologna, 40136, Italy</footnote>

<footnote>\* Corresponding author.</footnote>

<footnote>E-mail addresses: boris detienne@u-bordeaux.fr (B. Detienne), henri lefebvre@unibo.it (H. Lefebvre), enrico malaguti@unibo.it (E. Malaguti), michele monaci@unibo.it (M. Monaci).</footnote>
first-stage feasible region. Given $\bar{\mathbf{x}} \in \mathcal{X}$ and $\bar{\xi} \in \Xi$, the resulting second-stage feasible region is $\mathcal{Y}(\bar{\mathbf{x}}, \bar{\xi})$, and the second-stage problem is $\inf\{f(\bar{\xi}, \bar{\mathbf{x}}, \mathbf{y}) : \mathbf{y} \in \mathcal{Y}(\bar{\mathbf{x}}, \bar{\xi})\}$. It is known (see, Ben-Tal et al., 2004) that most of the problems which can be cast as two-stage robust problems are at least NP-hard. This result even holds for cases where first and second-stage variables are continuous and all the involved functions are linear. Several approaches have been developed to tackle this class of problems. Assuming that the second stage is continuous and exhibits strong duality, it can be replaced by its dual. This way, the inner maximization problem can be reformulated using its epigraph, leading to a constraint-generation algorithm in the spirit of Benders' decomposition (see, e.g., Bertsimas, Litvinov, & Sun, 2013; Jiang, Zhang, Li, & Guan, 2014; Terry, Epelman, & Thiele, 2009 and Gabrel, Lacroix, Murat, & Remli, 2011). A column-and-constraint-generation scheme has been proposed in Zeng & Zhao (2013), which consists in adding one set of second-stage decision variables and constraints associated with each realization of the uncertainty. These realizations are dynamically generated by solving a so-called adversarial problem which identifies the worst-case scenario for a current estimate of the first-stage decisions. The algorithm stops when no such scenario can be found. Later, the same approach was used in Ayoub & Poss (2016), where the authors model the adversarial problem as a mixed integer program, derived using Farkas' lemma and standard linearization techniques. Note that the column-and-constraint generation approach can handle mixed-integer second-stage decisions, which is not the case for classical Benders-type approaches. Unfortunately, this method seems to be of practical relevance only when a small number of variables has to be added for reaching optimality.

The inherent difficulty of this class of problems motivated the development of approximate solution methods. In the *affine decision rule* approach (Ben-Tal et al., 2004), the second-stage decisions are expressed as affine functions of the uncertainty. Another relevant approach, introduced in Bertsimas & Caramanis (2010), is the *finite adaptability* (also known as K-adaptability) in which the number of second-stage decisions is restricted to some finite number. An MILP formulation for the case of binary second-stage decisions and objective uncertainty was proposed in Hanasusanto, Kuhn, & Wiesemann (2015) and a branch-and-bound algorithm was later proposed in Subramanyam, Gounaris, & Wiesemann (2019) to address cases with uncertain linear constraints.

An important special case of (2SR-P) arises when uncertainty affects the objective function only, i.e., $\mathcal{Y}(\mathbf{x}, \xi) = \mathcal{Y}(\mathbf{x}), \forall \xi \in \Xi$. For this specific case, Kämmerling & Kurtz (2020) proposed an oracle-based algorithm relying on a hull relaxation combining the first- and second-stage feasible regions embedded within a branch-and-bound framework. However, this approach applies to first-stage binary variables and linear constraints only. On the other hand, Arslan & Detienne (2022) proposed an exact MILP reformulation of the problem in case of linear linking constraints that involve binary variables only. Besides solving the problem by means of a branch-and-price algorithm, a further contribution of Arslan & Detienne (2022) is proving the NP-completeness of the problem in this setting.

Our analysis shows that, in the setting where uncertainty affects the objective function only, no contribution has been presented in the literature for tackling problems where linking constraints are defined by nonlinear functions or involve both integer and continuous variables. Similarly, to the best of our knowledge, the case in which the objective function is nonlinear has not been considered yet. This paper contributes in filling this gap, as we consider two-stage robust problems with objective uncertainty, convex constraints and mixed-integer first and second stage. By extending in a non-trivial way some recent results from the two-stage stochastic optimization literature (see Sherali & Frat-

icelli, 2002; Sherali & Zhu, 2006 and Li & Grossmann, 2019), we obtain a relaxation of the problem, and analyze its tightness for different special cases. This relaxation can be embedded within an enumerative scheme thus producing an exact solution approach, for which we prove asymptotic convergence in the general case, and finite convergence in the integer case. Besides the theoretical analysis, we also show that, from a computational viewpoint, the proposed algorithm is able to solve instances of practical relevance arising from the logistic field. We also point out that the class of problems which can be addressed by our solution approach is quite large since we only require mild assumptions on the nature of the involved optimization problem.

The article is organized as follows. In Section 2 we formally introduce the considered class of problems, whereas in Section 3 we present a relaxation of the problem. We then derive sufficient conditions for the relaxation to coincide with the original problem in a mixed-integer context. So as to close the optimality gap, we introduce an enumerative algorithm which embeds a spatial branching mechanism on continuous first-stage variables. We prove asymptotic convergence of the overall algorithm in presence of continuous first-stage decisions and finite $\varepsilon$-convergence in case of integer first-stage decisions. In Section 3.4, we propose a column-generation algorithm to solve the relaxation problem. Finally, Section 4 applies the proposed method to a capacitated facility location problem with congestion.

**Notations** Throughout this paper, matrices and vectors are written in bold case, e.g., $\mathbf{x} \in \mathbb{R}^n$ or $\mathbf{A} \in \mathbb{R}^{n \times m}$, while components are written in normal font, e.g., $x_i$ or $a_{ij}$. Columns of $\mathbf{A}$ are written in bold case with exponent indexing, e.g., $\mathbf{a}^i$. Let $f: \mathbb{R}^n \to \mathbb{R}$ be a given function with $\text{dom}(f) = \{\mathbf{x} \in \mathbb{R}^n : f(\mathbf{x}) < +\infty\}$; its convex conjugate is denoted by $f^*: \mathbb{R}^n \to \mathbb{R}$ and is given by

$$f^*(\pi) = \sup_{\mathbf{x} \in \text{dom}(f)} \{\pi^T \mathbf{x} - f(\mathbf{x})\}.$$

Similarly, we denote by $f_*$ the concave conjugate of $f$. Let $X \subseteq \mathbb{R}^n \times \mathbb{Z}^{n-p}$ be a given set, described in terms of constraints and integer restrictions. We denote by $\bar{X}$ its continuous relaxation and by $\text{conv}(X)$ its convex hull, i.e., the smallest convex set $C$ satisfying $X \subseteq C$.

The indicator function of $X$ is noted $\delta(\cdot|X)$ and equals zero if its argument belongs to $X$ and $+\infty$ otherwise. Its convex conjugate is therefore given by $\delta^*(\pi^T X : \mathbf{x} \in X)$. Basic results on conjugate calculus are summarized in Appendix A. If $X$ is a convex polytope, we note $\text{vert}(X)$ the set of its extreme points. Finally, for a logical proposition $\mathcal{E}$, function $\mathbf{1}(\mathcal{E})$ equals one if $\mathcal{E}$ is true and zero otherwise.

## 2. Problem description

### 2.1. General setting

As anticipated, our goal is to solve problem (2SR-P) with objective uncertainty, convex constraints and mixed-integer first and second stages.

For the sake of clarity, let us first introduce several sets. Set $I = \{1, \dots, n_1\}$ denotes the set of indices for the first-stage variables, partitioned into two sets $I_l$ and $I_c$. Variables whose index belongs to $I_l$ are required to take integer values, while those whose index belongs to $I_c$ are continuous variables, i.e., $wlog, \mathcal{X} \subset \mathbb{R}^{|I_l|} \times \mathbb{Z}^{|I_c|}$. Similarly, we introduce set $J = \{1, \dots, n_2\}$ as the indices for the second-stage variables and partition this set into $J_l$ and $J_c$, i.e., $wlog, \mathcal{Y} \subset \mathbb{R}^{|J_l|} \times \mathbb{Z}^{|J_c|}$. Finally, we introduce set $U = \{1, \dots, n_3\}$ as the index set for the uncertain variables, i.e., $\Xi \subset \mathbb{R}^{n_3}$.

We now explicit some assumptions on the problem.

**Assumption 1 (Objective uncertainty).** For all $\xi \in \Xi$ and $\mathbf{x} \in \bar{\mathcal{X}}$, $\mathcal{Y}(\xi, \mathbf{x}) = \mathcal{Y}(\mathbf{x})$.
**Assumption 2 (Convexity).**

1. $\bar{x}$ is compact and convex;

2. The uncertainty set $\Xi$ is a finite-dimensional, bounded convex set;

3. For all $\mathbf{x} \in \bar{\mathcal{X}}$, $\bar{\mathcal{Y}}(\mathbf{x})$ is a finite-dimensional, bounded convex set;

4. The objective function $f$ is a concave function of the uncertainty and a convex function of the first- and second-stage decisions, i.e., $f_{\mathbf{x},\mathbf{y}}: \xi \mapsto f(\xi, \mathbf{x}, \mathbf{y})$ is a concave function for all fixed $\mathbf{x} \in \bar{\mathcal{X}}$ and $\mathbf{y} \in \bar{\mathcal{Y}}(\mathbf{x})$, and $f_{\xi}: (\mathbf{x}, \mathbf{y}) \mapsto f(\xi, \mathbf{x}, \mathbf{y})$ is a convex function for all fixed $\xi \in \Xi$.

**Assumption 3 (Complete recourse).** For every (relaxed) first-stage decision, there exists at least one feasible second-stage decision, i.e., for every $\mathbf{x} \in \bar{\mathcal{X}}$, $\mathcal{Y}(\mathbf{x})$ is a non-empty set.

**Assumption 4 (Boundedness).**

1. The objective function $f$ is bounded over the first- and second-stage feasible region, i.e., for all fixed $\xi \in \Xi$, $\{(\mathbf{x}, \mathbf{y}) : \mathbf{x} \in \bar{\mathcal{X}}, \mathbf{y} \in \bar{\mathcal{Y}}(\mathbf{x})\} \subseteq \text{dom}(f_{\xi})$;

2. For all $(\mathbf{x}, \mathbf{y}) : \mathbf{x} \in \bar{\mathcal{X}}$ and $\mathbf{y} \in \bar{\mathcal{Y}}(\mathbf{x})$, $\text{relint}(\Xi) \cap \text{dom}(f_{\mathbf{x},\mathbf{y}}) \neq \emptyset$.

**Assumption 5 (Separability).** Let $Q = \{1, \dots, q\}$.

1. The objective function $f$ can be expressed as a sum of $q$ functions, i.e., there exist $q$ functions $(\psi_i : \mathbb{R}^{|U|+|V|} \to \mathbb{R})_{i \in Q}$ such that $f(\xi, \mathbf{x}, \mathbf{y}) = \sum_{i \in Q} \psi_i(\xi, \mathbf{x}, \mathbf{y})$ for all $\mathbf{x} \in \mathcal{X}, \mathbf{y} \in \mathcal{Y}(\mathbf{x})$ and all $\xi \in \Xi$;

2. For all $i \in Q$, $\psi_i$ is separable in $\xi$ and $(\mathbf{x}, \mathbf{y})$ meaning that there exists functions $(w_i : \mathbb{R}^{|U|} \to \mathbb{R})_{i \in Q}$ and $(\varphi_i : \mathbb{R}^{|V|} \to \mathbb{R})_{i \in Q}$ such that $\psi_i(\xi, \mathbf{x}, \mathbf{y}) = w_i(\xi)\varphi_i(\mathbf{x}, \mathbf{y})$. In addition, we assume that $w_i(\cdot)$ is a concave function and $\varphi_i(\cdot)$ is a convex function.

A few remarks regarding these assumptions are necessary. First, note that **Assumptions 1 and 2** are here to define what we refer to as convex mixed-integer robust problems with objective uncertainty. We highlight that the word “convex” is here to suggest that all involved functions are convex with respect to the first- and second-stage variables. Yet, in general, even under these assumptions, problem (2SR-P) may fail to have a straightforward convex MINLP formulation. Indeed, function $h : \mathbf{x} \mapsto \max_{\xi \in \Xi} \min_{y \in \mathcal{Y}(\mathbf{x})} f(\xi, \mathbf{x}, \mathbf{y})$ is not necessarily a convex function over the continuous relaxation of $\mathcal{X}$. We give here a small example.

**Example 1 (nonconvex MINLP).** Consider the following first- and second-stage feasible regions.

$$ \mathcal{X} = [0, 1] \text{ and } \mathcal{Y}(\mathbf{x}) = \left\{ \mathbf{y} \in \{0, 1\}^2 \middle| \begin{array}{l} y_1 + y_2 \le 1 \\ y_1 \le 1 - x \end{array} \right\}. $$

By inspection, we have that $(y_1, y_2) = (0, 0)$ and $(y_1, y_2) = (0, 1)$ are always feasible second-stage solutions, while $(y_1, y_2) = (1, 0)$ is feasible only when $x = 0$. Fixing the uncertainty set $\Xi = [0, 1]$, we take interest in the following convex mixed-integer two-stage robust problem

$$ \min_{x \in [0,1]} h(\mathbf{x}) \text{ with } h : \mathbf{x} \mapsto \max_{\xi \in [0,1]} \min_{(y_1, y_2) \in \mathcal{Y}(\mathbf{x})} \xi(-2y_1 + y_2 + 1). $$

Though every involved functions are convex (in fact, affine) with respect to $\xi, \mathbf{x}$ and $\mathbf{y}$, we have that

$$ h(\mathbf{x}) = \begin{cases} \max_{\xi \in [0,1]} \min\{\xi; 2\xi; -\xi\} = 0 & \text{if } x = 0 \\ \max_{\xi \in [0,1]} \min\{\xi; 2\xi\} = 1 & \text{if } x > 0 \end{cases} = \mathbf{1}(x > 0). $$

Clearly, $h$ fails to be convex over $[0,1]$ which ends our example. More generally, function $h$ can be seen as the optimal value function of a minimization problem whose objective function is convex

but in which the feasible region is defined by a non-convex set-valued map (this will be later proven in **Lemma 2**).

**Assumptions 3** is a standard assumption in the two-stage optimization literature, and is known to be easy to enforce as soon as the considered problem is bounded, which is implied by **Assumption 4.1. Assumption 4.2** is not restrictive in practice, and will be used in the proof of **Lemma 2**.

Finally, **Assumption 5** is structural to our work, and implies the following remarks.

**Remark 1.** The assumption that $\varphi_i(\cdot)$ is a convex function (at most affine) for all $i \in Q$ is without loss of generality.

**Proof.** Let $i \in Q$ such that $\varphi_i(\cdot)$ is concave, then, to fulfill **Assumption 2.4**, $w_i(\xi)$ must be negative for all $\xi \in \Xi$. Thus, one may equivalently replace $w_i(\cdot)$ by $-w_i(\cdot)$ and $\varphi_i(\cdot)$ by $-\varphi_i(\cdot)$. $\square$

**Remark 2.** For all $i \in Q$ such that $\varphi_i(\cdot)$ (resp. $w_i(\cdot)$) is not single-signed, then $w_i(\cdot)$ (resp. $\varphi_i(\cdot)$) is affine.

**Remark 3.** For all $i \in Q$ such that $\varphi_i(\cdot)$ (resp. $w_i(\cdot)$) is not affine, then $w_i(\cdot)$ (resp. $\varphi_i(\cdot)$) is a non-negative function.

Note that **Assumption 5** could be relaxed to address situations in which, for $i \in Q$ such that $\varphi_i(\cdot)$ (resp. $w_i(\cdot)$) is affine, there is no restriction on the concavity (resp. convexity) of the associated $w_i(\cdot)$ (resp. $\varphi_i(\cdot)$).

**Example 2 (Fulfilling Assumption 5).** We give here some examples of functions which satisfy **Assumption 5**. For simplicity, we denote $z = (\mathbf{x}, \mathbf{y})$.

* Uncertain linear functions of the form $(\xi, z) \mapsto \xi Az$ where $A$ is a given real matrix;

* Diagonal uncertain convex quadratic form $(\xi, z) \mapsto z^T \text{diag}(\xi) z$ where $\xi \ge 0$;

* Uncertain positively weighted sum of convex functions of the form $(\xi, z) \mapsto \sum_{i \in Q} \xi_i \varphi_i(\mathbf{z})$ with $\Xi \subset \mathbb{R}_+^{|U|}$, e.g., $(\xi, \mathbf{x}, \mathbf{y}) \mapsto \sum_{i \in Q} \xi_i x_i^2 / y_i$ with $\mathbf{y} \ge \mathbf{0}$.

**Example 3 (Violating Assumption 5).** We give here some examples of functions which do not satisfy **Assumption 5**.

* Non-concave functions of the uncertainty, e.g., $(\xi, z) \mapsto \|z - \xi\|$ for any given norm;

* General uncertain quadratic form $(\Sigma, z) \mapsto z^T \Sigma z$ even with $\Sigma \ge 0$ (unless $\Xi \cap \mathbb{R}_+^{|U|} \neq \emptyset$).

In the following lemma, we finally state the class of problems we consider.

**Lemma 1.** Under **Assumptions 1–5**, there exists $[\mathbf{l}, \mathbf{u}] \subset \mathbb{R}^{|U|+|V|}$ such that (2SR-P) is equivalent to the following two-stage optimization problem with convex objective function and objective uncertainty

$$ \inf_{\mathbf{x} \in \mathcal{X}} \sup_{[\mathbf{l}, \mathbf{u}] \in \Xi} \inf_{(\mathbf{t}, \mathbf{y}) \in \mathcal{Y}(\mathbf{x})} \sum_{i \in Q} w_i(\xi) t_i \quad (2SRO-P) $$

with $\mathcal{Y}(\mathbf{x})$ such that $\mathcal{Y}(\mathbf{x}) = \text{proj}_{\mathbf{y}}(\mathcal{Y}'(\mathbf{x}))$ and $\bar{\mathcal{Y}}'(\mathbf{x})$ is a convex and finite-dimensional set.

**Proof.** The existence of the hyper-rectangle $[\mathbf{l}, \mathbf{u}]$ is trivial as $\mathcal{X}$ is assumed to be bounded (**Assumption 2.1**). Moreover, the following equality holds.

$$ \begin{aligned} & \inf_{\mathbf{y}} \left\{ \sum_{i \in Q} w_i(\xi) \varphi_i(\mathbf{x}, \mathbf{y}) : \mathbf{y} \in \mathcal{Y}(\mathbf{x}) \right\} \\ & = \inf_{\mathbf{y}, t} \left\{ \sum_{i \in Q} w_i(\xi) t_i : \mathbf{y} \in \mathcal{Y}(\mathbf{x}), t_i = \varphi_i(\mathbf{x}, \mathbf{y}), \forall i \in Q \right\} \end{aligned} $$
However, the optimization problem on the right side of the equality may fail to be convex if there exists $i \in \mathcal{Q}$ such that $\varphi_i$ is not affine. Let $\mathcal{Q}^A \subseteq \mathcal{Q}$ be the set of indices for which $\varphi_i$ is affine. By Assumption 5, for all $i \in \mathcal{Q} \setminus \mathcal{Q}^A$, we have $w_i(\cdot) \ge 0$ and thus constraint $"t_i = \varphi_i(\mathbf{x}, \mathbf{y})"$ may be equivalently replaced by $"t_i \ge \varphi_i(\mathbf{x}, \mathbf{y})"$, which is convex. We therefore can choose

$$ \mathcal{Y}'(\mathbf{x}) = \left\{ \begin{array}{ll} \mathbf{y} \in \mathcal{Y}(\mathbf{x}) & \\ (\mathbf{t}, \mathbf{y}) : t_i = \varphi_i(\mathbf{x}, \mathbf{y}) & \forall i \in \mathcal{Q}^A \\ t_i \ge \varphi_i(\mathbf{x}, \mathbf{y}) & \forall i \in \mathcal{Q} \setminus \mathcal{Q}^A \end{array} \right\}. $$

For every $\mathbf{x} \in \bar{\mathcal{X}}$, the continuous relaxation of $\mathcal{Y}'(\mathbf{x})$ is convex and non-empty (Assumption 3); by construction, it is also finite dimensional. □

In what remains, we will assume to know a hyper-rectangle $[\mathbf{l}, \mathbf{u}]$ as described in Lemma 1.

## 2.2. Special case: Linear and binary setting

We complete the introduction by discussing the special case of (2SRO-P) under the following additional assumptions:

1. $\mathcal{X}, \Xi$ and $\mathbf{x} \mapsto \mathcal{Y}(\mathbf{x})$ are defined by linear constraints;

2. there exists a matrix $\mathbf{A} \in \mathbb{R}^{[U] \times [Q]}$ such that $\forall i \in \mathcal{Q}, w_i(\xi) = \xi^T \mathbf{d}^i$; and

3. linking constraints are defined by functions $g(\mathbf{x}, \mathbf{y})$ that do not depend on first-stage variables in $I_c$.

In a recent paper Arslan & Detienne (2022), the authors observed that, for this variant of the problem, the inner minimization $\min_{\mathbf{y} \in \text{conv}(\mathcal{Y}(\mathbf{x}))} \xi^T \mathbf{A} \mathbf{y}$ can be equivalently replaced by $\min_{\mathbf{y} \in \text{conv}(\mathcal{Y}(\mathbf{x}))} \xi^T \mathbf{A} \mathbf{y}$, i.e., the second-stage feasible region can be substituted by its convex hull. This allows to transform the min-max-min problem into a min-max problem by the well known minimax theorem. Assuming that $\Xi$ is expressed as $\{\xi \in \mathbb{R}_+^{|U|} : F\xi \le d\}$, the inner maximization problem is dualized so as to obtain the following equivalent problem

$$ \min_{\mathbf{x}, \mathbf{y}, \mathbf{\lambda}} \mathbf{d}^T \mathbf{\lambda} \quad (1) $$

$$ \text{subject to} \quad \mathbf{x} \in \mathcal{X} \quad (2) $$

$$ \mathbf{y} \in \text{conv}(\mathcal{Y}'(\mathbf{x})) \quad (3) $$

$$ \mathbf{F}^T \mathbf{\lambda} \ge \mathbf{A} \mathbf{y} \quad (4) $$

$$ \mathbf{\lambda} \ge \mathbf{0}, \quad (5) $$

where $\mathbf{\lambda}$ are the dual variables associated to the inner maximization problem. Note that, besides the integrality requirements on the variables, the only nonconvex constraint is (3). By exploiting a reformulation already used in Sherali & Fratilcelli (2002), Sherali & Zhu (2006) and Li & Grossmann (2019) for two-stage stochastic optimization problems with mixed-integer first and second stage, Arslan & Detienne (2022) showed that, for each fixed $\bar{\mathbf{x}} \in \mathcal{X}$, $\{\bar{\mathbf{x}}\} \times \text{conv}(\mathcal{Y}(\bar{\mathbf{x}})) = \text{conv}(\mathcal{S}) \cap \{\mathbf{x}, \mathbf{t}, \mathbf{y} : \mathbf{x} = \bar{\mathbf{x}}\}$ where $\mathcal{S} = \{\mathbf{x}, \mathbf{y} : \mathbf{x} \in \{0, 1\}, \mathbf{y} \in \mathcal{Y}(\mathbf{x})\}$. Hence, constraint (3) may be equivalently enforced as

$$ (\mathbf{x}, \mathbf{y}) \in \text{conv}(\mathcal{S}). \quad (6) $$

The obtained reformulation is then solved by means of a branch-and-price algorithm where branching is performed on the first-stage variables only.

## 3. A hull-relaxation-based branch-and-price algorithm

In this section we present our main contribution and its theoretical foundations. We first turn problem (2SRO-P) from a min-max-min problem to a min-max problem in our mixed-integer and convex context. Then, since linear duality does not apply in our

setting, we resort to Fenchel duality to obtain a reformulation of the problem. Similarly to the linear and binary case, we then replace the counterpart of (3) by constraints which play the same role as (6). This only provides a relaxation of the problem in the general setting. This relaxation is thus embedded into an enumerative scheme to obtain an optimal solution of (2SRO-P).

### 3.1. Problem reformulation

The following lemma extends the result given in Arslan & Detienne (2022) to the mixed-integer and convex context.

**Lemma 2 (Single-stage reformulation).** *Problem (2SRO-P) is equivalent to the following problem*

$$ \inf_{(\mathbf{x}, \mathbf{t}, \mathbf{y}) \in \mathcal{F}} \sup_{\xi \in \Xi} \sum_{i \in \mathcal{Q}} w_i(\xi) t_i \quad (7) $$

$$ \text{with } F = \{ (\mathbf{x}, \mathbf{t}, \mathbf{y}) : \mathbf{x} \in \mathcal{X} \cap [\mathbf{l}, \mathbf{u}], (\mathbf{t}, \mathbf{y}) \in \text{conv}(\mathcal{Y}'(\mathbf{x})) \}. $$

**Proof.** This lemma relies on the same arguments as those employed in Arslan & Detienne (2022): first, the feasible region of the inner minimization problem is replaced by its convex hull. This is valid by linearity of the objective function. By Assumption 2.2 and Lemma 2, both $\Xi$ and $\text{conv}(\mathcal{Y}'(\mathbf{x}))$ (for all $\mathbf{x} \in \mathcal{X}$) are convex and finite dimensional sets. Thus, the minmax theorem in Perchet & Vigeral (2015) can be used to turn the inner $\sup$-inf into an $\inf$-sup problem. This achieves the proof. □

The inner maximization problem may be turned into a minimization problem by use of Fenchel duality, as done in Ben-Tal et al. (2009). In the following proposition, we therefore derive a general convex reformulation of problem (2SRO-P).

**Proposition 1 (Deterministic reformulation).** *Problem (2SRO-P) is equivalent to the following problem*

$$ \inf_{\mathbf{x}, \mathbf{y}, \mathbf{t}, (\mathbf{v})_{i \in \mathcal{Q}}} \delta^*(\xi | \Xi) - \sum_{i \in \mathcal{Q}} (t_i w_i)_* (\mathbf{v}^i) \quad (8) $$

$$ \text{subject to} \quad \mathbf{x} \in \mathcal{X} \cap [\mathbf{l}, \mathbf{u}] \quad (9) $$

$$ (\mathbf{t}, \mathbf{y}) \in \text{conv}(\mathcal{Y}'(\mathbf{x})) \quad (10) $$

$$ \sum_{i \in \mathcal{Q}} \mathbf{v}^i = \xi \quad (11) $$

$$ \xi \in \mathbb{R}^{|U|} \quad (12) $$

$$ \mathbf{v}^i \in \mathbb{R}^{|U|} \quad \forall i \in \mathcal{Q}. \quad (13) $$

**Proof.** By a direct application of Fenchel duality and some conjugate calculus results, the following holds.

$$ \begin{aligned} \sup_{\xi \in \Xi} \sum_{i \in \mathcal{Q}} t_i w_i(\xi) &= \sup_{\xi \in \mathbb{R}^{|U|}} \left\{ \sum_{i \in \mathcal{Q}} t_i w_i(\xi) - \delta(\xi | \Xi) \right\} \\ &= \inf_{\xi \in \mathbb{R}^{|U|}} \left\{ \delta^*(\xi | \Xi) - \left( \sum_{i \in \mathcal{Q}} t_i w_i(\xi) \right)_* \right\} \\ &= \inf_{\xi \in \mathbb{R}^{|U|}} \left\{ \delta^*(\xi | \Xi) - \sup_{\mathbf{v}^i \in \mathbb{R}^{|U|}, i \in \mathcal{Q}} \left\{ \sum_{i \in \mathcal{Q}} (t_i w_i)_* (\mathbf{v}^i) : \sum_{i \in \mathcal{Q}} \mathbf{v}^i = \xi \right\} \right\} \\ &= \inf \left\{ \delta^*(\xi | \Xi) - \sum_{i \in \mathcal{Q}} (t_i w_i)_* (\mathbf{v}^i) : \sum_{i \in \mathcal{Q}} \mathbf{v}^i = \xi, \mathbf{v}^i \in \mathbb{R}^{|U|}, i \in \mathcal{Q}, \xi \in \mathbb{R}^{|U|} \right\}. \end{aligned} $$

See also Appendix A for more details on conjugate calculus. □

The following results show that, although the reformulation for the general case adds $|Q| \times |U|$ continuous variables, these additional variables can be omitted for some relevant cases. In particular this is true in case all the $w_i(\cdot)$ functions are either separable or affine.
<figure><img src="image_4.png" /><figcaption>Fig. 1. Graphical representation of different sets from Example 1.</figcaption></figure>

**Remark 4.** Assume wlog that $|Q| = |U|$. If, for all $i \in Q$, $w_i(\xi) = w_i(\xi_i)$, then problem (2SRO-P) is equivalent to

$$ \inf_{(\mathbf{x}, \mathbf{t}, \mathbf{y}) \in \mathcal{F}} \left\{ \delta^*(\xi | \Xi) - \sum_{i \in Q} (t_i w_i)_*(\xi) \right\}. \quad (14) $$

**Proof.** By assumption, we have $(t_i w_i)_*(\mathbf{v}^i) = (t_i w_i)_*(v_{ii})$. In addition, (11) is as follows

$$ \sum_{i \in Q} v_{ij} = \xi_j \quad j = 1, \dots, |U|, \quad (15) $$

and $v_{ij}$, $i \neq j$, does not appear in the objective function. Thus, there always exists an optimal solution such that $\xi_j = v_{jj}$, $j = 1, \dots, |U|$. □

**Remark 5.** Let $i \in Q$ such that $w_i(\cdot)$ is affine, i.e., $w_i(\xi) = (\mathbf{r}^i)^T \xi + r_{i0}$. Problem (2SRO-P) is equivalent to

$$ \inf_{(\mathbf{x}, \mathbf{t}, \mathbf{y}) \in \mathcal{F}} \{ \delta^*(\mathbf{Rt} | \Xi) + r_{0i}^T \mathbf{t} \}. \quad (16) $$

**Proof.** Indeed, we have

$$ (t_i w_i)_*(\mathbf{v}) = \inf_{\xi \in \mathbb{R}^{|U|}} \{ \mathbf{v}^T \xi - t_i ((\mathbf{r}^i)^T \xi + r_{i0}) \} = \begin{cases} -t_i r_{i0} & \text{if } \mathbf{v} = t_i \mathbf{r}^i \\ -\infty & \text{otherwise.} \end{cases} \quad \square $$

## 3.2. Relaxation

Note that the deterministic reformulation presented above still is not, in general, a convex MINLP. Indeed, $\mathrm{conv}(\mathcal{Y}'(\mathbf{x}))$ in constraints (10) depends on variables $\mathbf{x}$. Since no tractable compact form is known in the general case, we replace constraint $(\mathbf{t}, \mathbf{y}) \in \mathrm{conv}(\mathcal{Y}'(\mathbf{x}))$ by the following relaxed requirement.

$$ (\mathbf{x}, \mathbf{t}, \mathbf{y}) \in \mathrm{conv}(S) \text{ with } S = \left\{ \begin{array}{ll} l_j \le x_j \le u_j & \forall j \in I \\ (\mathbf{x}, \mathbf{t}, \mathbf{y}) : x_j \in \mathbb{Z} & \forall j \in I_l \\ (\mathbf{t}, \mathbf{y}) \in \mathcal{Y}'(\mathbf{x}) & \end{array} \right\}. \quad (17) $$

The substitution yields the following problem, which is a relaxation of (8)-(13)

$$ \begin{array}{ll} \min_{\mathbf{x}, \mathbf{y}, (\mathbf{v})_{i \in Q}, \xi} & \delta^*(\xi | \Xi) - \sum_{i \in Q} (t_i w_i)_*(\mathbf{v}^i) \\ \text{subject to} & \mathbf{x} \in \mathcal{X} \cap [\mathbf{l}, \mathbf{u}] \\ & (\mathbf{x}, \mathbf{t}, \mathbf{y}) \in \mathrm{conv}(S) \\ & \sum_{i \in Q} \mathbf{v}^i = \xi \\ & \mathbf{v}^i \in \mathbb{R}^{|U|} \quad \forall i \in Q \\ & \xi \in \mathbb{R}^{|U|}. \end{array} \quad (P) $$

It is clear that, for any fixed $\tilde{\mathbf{x}} \in \mathcal{X}$, we have $\{\tilde{\mathbf{x}}\} \times \mathcal{Y}'(\tilde{\mathbf{x}}) = S \cap \{(\mathbf{x}, \mathbf{t}, \mathbf{y}) : \mathbf{x} = \tilde{\mathbf{x}}\}$, and that the same holds even for $\tilde{\mathbf{x}} \in \tilde{\mathcal{X}}$. However,

as shown, e.g., in Sherali & Zhu (2006), the convexified counterpart does not hold, in the sense that the inclusion $\{\tilde{\mathbf{x}}\} \times \mathrm{conv}(\mathcal{Y}'(\tilde{\mathbf{x}})) \subseteq \mathrm{conv}(S) \cap \{(\mathbf{x}, \mathbf{t}, \mathbf{y}) : \mathbf{x} = \tilde{\mathbf{x}}\}$" may be strict. Example 4 below illustrates this case.

**Example 4** (Hull relaxation). We consider the first- and second-stage feasible sets introduced in Example 1. In Fig. 1a, we represent the convex hull of S. For a fixed first-stage decision $\tilde{\mathbf{x}}$ (here, $\tilde{\mathbf{x}} = 0.4$), Fig. 1b reports the feasible points for constraint (17), whereas Fig. 1c describes the exact shape of $\mathrm{conv}(\mathcal{Y}'(\tilde{\mathbf{x}}))$. The figure shows an example in which inclusion is strict. In addition, note that, whenever $\tilde{\mathbf{x}}$ attains its bounds (i.e., $\tilde{\mathbf{x}} \in \{0, 1\}$), $\{\tilde{\mathbf{x}}\} \times \mathrm{conv}(\mathcal{Y}'(\tilde{\mathbf{x}})) = \mathrm{conv}(S) \cap \{(\mathbf{x}, \mathbf{y}) : \mathbf{x} = \tilde{\mathbf{x}}\}$ holds.

The following Lemma follows from the considerations above.

**Lemma 3** (Lower-bounding property). Denoting by $v(\bullet)$ the optimal objective value of problem $\bullet$, we have

$$ v(P) \le v(2SRO - P). $$

In other words, (P) is a relaxation of (2SRO-P). In the next proposition, we introduce a condition under which a feasible solution for problem (P) is feasible for problem (2SRO-P) as well.

**Proposition 2.** If $\tilde{\mathbf{x}} \in \mathrm{vert}([\mathbf{l}, \mathbf{u}])$, then

$$ \{\tilde{\mathbf{x}}\} \times \mathrm{conv}(\mathcal{Y}'(\tilde{\mathbf{x}})) = \mathrm{conv}(S) \cap \{(\mathbf{x}, \mathbf{t}, \mathbf{y}) : \mathbf{x} = \tilde{\mathbf{x}}\}. $$

**Proof.** Let $\tilde{\mathbf{x}} \in \mathrm{vert}([\mathbf{l}, \mathbf{u}])$ and let $(\hat{\mathbf{x}}, \hat{\mathbf{t}}, \hat{\mathbf{y}}) \in \mathrm{conv}(S) \cap \{(\mathbf{x}, \mathbf{t}, \mathbf{y}) : \mathbf{x} = \tilde{\mathbf{x}}\}$. Then, $(\hat{\mathbf{x}}, \hat{\mathbf{t}}, \hat{\mathbf{y}})$ can be expressed as a (finite) convex combination of points of $\mathrm{conv}(S)$ (Carathéodory's theorem), i.e.,

$$ (\hat{\mathbf{x}}, \hat{\mathbf{t}}, \hat{\mathbf{y}}) = \sum_{e \in E} (\tilde{\mathbf{x}}^e, \tilde{\mathbf{t}}^e, \tilde{\mathbf{y}}^e) \alpha_e, $$

where E is a given index list of such elements of $\mathrm{conv}(S)$. Assume that there exists $j \in I$ and $i \in E$ such that $\tilde{\mathbf{x}}_j^i \neq \tilde{\mathbf{x}}_j$. If $\tilde{\mathbf{x}}_j^i > \tilde{\mathbf{x}}_j$, condition $\tilde{\mathbf{x}}_j^i \in \mathrm{conv}(S)$ implies that $\tilde{\mathbf{x}}_j = l_j$. Hence, $\alpha_i = 0$ since $\tilde{\mathbf{x}}_j^k \ge l_j \forall k \in E$. The same argument shows that $\tilde{\mathbf{x}}_j^i < \tilde{\mathbf{x}}_j$ implies $\alpha_i = 0$. Thus, for each $e \in E$ such that $\alpha_e > 0$, we must have $\tilde{\mathbf{x}}_e^i = \tilde{\mathbf{x}}$. This implies that $(\tilde{\mathbf{t}}^e, \tilde{\mathbf{y}}^e) \in \mathcal{Y}'(\tilde{\mathbf{x}})$ and thus $\sum_{e \in E} (\tilde{\mathbf{t}}^e, \tilde{\mathbf{y}}^e) \alpha_e \in \mathrm{conv}(\mathcal{Y}'(\tilde{\mathbf{x}}))$. □

**Corollary 1** (Tightness condition). Let $X^*$ be the set of optimal first-stage decisions of problem (P). Then

$$ X^* \cap \mathrm{vert}([\mathbf{l}, \mathbf{u}]) \neq \emptyset \Rightarrow v(P) = v(2SRO - P). $$

**Proof.** Let $(\mathbf{x}^*, \mathbf{t}^*, \mathbf{y}^*)$ be an optimal solution of (P) with $\mathbf{x}^* \in \mathrm{vert}([\mathbf{l}, \mathbf{u}])$. From Proposition 2, it is also feasible for problem (2SRO-P). Thus, Lemma 3 implies optimality for problem (2SRO-P). □

This result directly implies Corollary 2 which states that, in the special case where the first-stage variables are all binary, problem (P) is always an exact reformulation of (2SRO-P).
**Corollary 2** (Tightness condition/binary case). If the first-stage decisions are all binary, i.e., $l_C = \emptyset$ and $\forall j \in I_l, l_j = 0, u_j = 1$, then

$$v(P) = v(2SRO - P).$$

**Proof.** In this case, $[l, u] = [0, 1]$, hence any optimal first-stage solution $\mathbf{x}^*$ satisfies $\mathbf{x}^* \in \{0, 1\}^{|I_l|} = \text{vert}([l, u])$ which, by Corollary 1, proves the result. □

### 3.3. Enumerative algorithm

We now present an exact method for solving problem (2SRO-P) which works as follows:

* we exploit the deterministic reformulation (8)-(13) of the problem;

* we relax requirement (10), i.e., $(\mathbf{t}, \mathbf{y}) \in \text{conv}(\mathcal{Y}^{\prime}(\mathbf{x}))$, and replace it by constraint (17) imposing that $(\mathbf{x}, \mathbf{t}, \mathbf{y}) \in \text{conv}(S)$, thus obtaining the lower-bounding problem (P);

* we relax the integrality of the $\mathbf{x}$, and impose that $\mathbf{x} \in \bar{\mathcal{X}}$;

* we solve this relaxation, which requires to provide a description of the convex hull of set $S$, by means of column generation techniques;

* as the resulting solution may violate the relaxed requirements, the column generation scheme is embedded within an enumerative algorithm which branches on first-stage variables only;

* our algorithm first branches on integer $\mathbf{x}$ variables having a fractional value, and then possibly resorts to spatial branching on continuous $\mathbf{x}$ variables, until each continuous variable attains either its lower or upper bound;

* when this is the case, according to Corollary 1, the current solution is optimal for the actual subproblem.

The resulting branch-and-price algorithm stores the best feasible solution found (the incumbent solution) which is returned when the method stops.

#### 3.3.1. Node solution

Let $p$ denote a generic node of the branching tree, associated with bounds $I^p$ and $\mathbf{u}^p$ on first-stage variables.

A lower bound on the optimal solution value of node $p$ can be computed solving the following problem

$$\min_{\mathbf{x}, \mathbf{t}, \mathbf{y}} \delta^*(\xi | \Xi) - \sum_{i \in Q} (t_i w_i)_* (\mathbf{v}^i)$$

$$\text{subject to} \quad \begin{array}{l} \mathbf{x} \in \bar{\mathcal{X}} \cap [I^p, \mathbf{u}^p] \\ (\mathbf{x}, \mathbf{t}, \mathbf{y}) \in \text{conv}(S^p) \\ \sum_{i \in Q} \mathbf{v}^i = \xi \\ \mathbf{v}^i \in \mathbb{R}^{|U|}, \quad \forall i \in Q \\ \xi \in \mathbb{R}^{|U|}, \quad \text{LB}^p \end{array}$$

where $S^p = \{(\mathbf{x}, \mathbf{t}, \mathbf{y}) : I^p \le \mathbf{x} \le \mathbf{u}^p, x_j \in \mathbb{Z}, \forall j \in I_l, (\mathbf{t}, \mathbf{y}) \in \mathcal{Y}^{\prime}(\mathbf{x})\}$. This problem is exactly the continuous relaxation of problem (P) where the bounds $l$ and $u$ have been replaced by $I^p$ and $\mathbf{u}^p$. Note that at the root node we have $l^0 = l$ and $u^0 = u$.

Let $(\mathbf{x}^{p*}, \mathbf{t}^{p*}, \mathbf{y}^{p*}, (\mathbf{v}^{p*})_{i \in Q}, \xi^{p*})$ be an optimal solution of problem $LB^p$. If $v(LB^p)$ is greater than or equal to the cost of the incumbent, the node is fathomed by bounding. Otherwise, we distinguish three cases:

* if $\mathbf{x}^{p*} \in \text{vert}([I^p, \mathbf{u}^p])$, by Proposition 2, this solution is optimal for the current node. Hence, the node is fathomed by optimality and the incumbent is updated;

* if $\mathbf{x}^{p*} \in \mathcal{X} \setminus \text{vert}([I^p, \mathbf{u}^p])$, we compute a feasible solution for (2SRO-P) by solving the following model

$$\begin{array}{ll} \min_{\mathbf{t}, \mathbf{y}, (\mathbf{v}^i)_{i \in Q}, \xi} & \delta^*(\xi | \Xi) - \sum_{i \in Q} (t_i w_i)_* (\mathbf{v}^i) \\ \text{subject to} & (\mathbf{t}, \mathbf{y}) \in \text{conv}(\mathcal{Y}^{\prime}(\mathbf{x}^{p*})) \\ & \sum_{i \in Q} \mathbf{v}^i = \xi \\ & \mathbf{v}^i \in \mathbb{R}^{|U|}, \quad \forall i \in Q \\ & \xi \in \mathbb{R}^{|U|}, \quad \text{LB}^p \end{array}$$

in which the first-stage variables are fixed to $\mathbf{x}^{p*}$. Note that, in this case, $\mathbf{x}^{p*}$ is a feasible first-stage solution; hence, by Assumption 3, problem $\text{LB}^p$ is always feasible, and possibly the incumbent is updated. If $v(\text{LB}^p) = v(\text{LB}^p)$ then node $p$ is solved; otherwise, we perform a branching;

* if $\mathbf{x}^{p*} \in \bar{\mathcal{X}} \setminus \mathcal{X}$, we branch.

In the last case, before branching, one can try to round $\mathbf{x}^{p*}$; if the resulting point is in $\mathcal{X}$, a feasible solution for (2SRO-P) can be computed.

#### 3.3.2. Branching

We now describe how to select the branching variable at node $p$. For each first-stage variable, say with index $j \in I_l$, we compute the minimum distance of $\mathbf{x}_j^{p*}$ from one of its bounds at the node, i.e., we evaluate

$$\theta_j^p = \begin{cases} \min\{\mathbf{x}_j^{p*} - \lfloor \mathbf{x}_j^{p*} \rfloor; \lfloor \mathbf{x}_j^{p*} \rfloor - \mathbf{x}_j^{p*} \rfloor \} & \text{if } j \in I_l \\ \min\{\mathbf{x}_j^{p*} - I_l^p; \mathbf{u}_j^p - \mathbf{x}_j^{p*} \} & \text{otherwise.} \end{cases}$$

For branching, we give priority to integer variables that do not attain their bound. Otherwise, we resort to spatial branching on continuous variables. In both cases, we select the variable with maximum $\theta_j^p$ value, i.e., we select variable $\mathbf{x}_j$ such that

$$\bar{j} \in \begin{cases} \operatorname{argmax}\{\theta_j^p : j \in I_l\} & \text{if } \exists j \in I_l : \theta_j^p > 0 \\ \operatorname{argmax}\{\theta_j^p : j \in I_c\} & \text{otherwise.} \end{cases}$$

If $\bar{j} \in I_l$, then a standard integer branching is executed. Otherwise, spatial branching generates two descendant nodes by imposing $\mathbf{x}_j^* \le \mathbf{x}_j^{p*}$ for the left node and $\mathbf{x}_j^* \ge \mathbf{x}_j^{p*}$ for the right one. We associate to each node the lower bound value of the current node $v(LB^p)$ and insert them in a list of open nodes. At each iteration, we extract from the list one node with minimum lower bound value, halting the algorithm when the list is empty.

**Example 5.** Fig. 2 illustrates the feasible region of the left and right child obtained by spatial branching on $x \le \beta$ and $x \ge \beta$, respectively, from Example 1 (here, $\beta = 0.4$). Clearly, the right child allows the same second-stage decisions as in $\mathcal{Y}^{\prime}(x)$ for all $x \ge \beta$. The left child, however, still allows second-stage decisions that could end up being infeasible in the original problem. In particular, $(\mathbf{x}, \mathbf{y}) = (\varepsilon, 1 - \varepsilon, 0)$ with $\varepsilon \in (0, \beta]$ is feasible for $(LB^p)$ but not for (2SRO-P).

#### 3.3.3. Convergence

In the following, we analyze the convergence of the branch-and-price algorithm. While finite convergence is ensured if all first-stage variables are integer, this may not be the case when the first-stage includes continuous variables. We now consider the case where our algorithm has an infinite number of nodes. Note that, in this case, there exists at least one infinite branch to the branching tree since the number of variables which can be selected for branching is finite. We consider one such branch and denote it by $P$. For each node $p \in P$, we denote by $(I^p, \mathbf{u}^p)$ the associated bounds
<figure><img src="image_5.png" /><figcaption>Fig. 2. Branching on continuous variable x from Example 1.</figcaption></figure>

<figure><img src="image_6.png" /></figure>

for the **x** variables and by $(\mathbf{x}^{p*}, \mathbf{t}^{p*}, \mathbf{y}^{p*}, \mathbf{V}^{p*}, \xi^{p*})$ an optimal solution to the lower-bounding problem. Additionally, we introduce the following function

$$f_{LB}(\mathbf{t}, \mathbf{V}, \xi) := \delta^*(\xi | \Xi) - \sum_{i \in Q} (t_i w_i)_* (\mathbf{v}^i), \quad (18)$$

which gives the value of the lower bounding problem (P) as a function of its arguments. Similarly, from Lemma 2, function

$$f_{ST}(\mathbf{t}) := \sup_{\xi \in \Xi} \sum_{i \in Q} w_i(\xi) t_i. \quad (19)$$

gives the value of the single stage reformulation as a function of argument **t**.

**Remark 6.** For each node $p \in P$ it holds $f_{LB}(\mathbf{t}^{p*}, \mathbf{V}^{p*}, \xi^{p*}) = f_{ST}(\mathbf{t}^{p*})$.

**Proof.** This directly follows from the definition of $(\mathbf{x}^{p*}, \mathbf{t}^{p*}, \mathbf{y}^{p*}, \mathbf{V}^{p*}, \xi^{p*})$ and Proposition (1). □

**Lemma 4.** Let P be a sequence of nodes of any infinite branch of the branching tree. Then,

(i) The sequence $\{(\mathbf{l}^p, \mathbf{u}^p)\}_{p \in P}$ has a unique accumulation point, which we denote by $(\mathbf{l}^*, \mathbf{u}^*)$;

(ii) The sequence $\{(\mathbf{x}^{p*}, \mathbf{t}^{p*}, \mathbf{y}^{p*})\}_{p \in P}$ has at least one accumulation point;

(iii) Let $\mathbf{x}^*$ be any accumulation point of $\{\mathbf{x}^{p*}\}_{p \in P}$, then, for each $j \in I_C$ which is infinitely selected for branching, there exists a sub-sequence $P^j \subseteq P$ such that either $\{u_j^p\}_{p \in P^j} \to x_j^*$ or $\{l_j^p\}_{p \in P^j} \to x_j^*$;

(iv) Every accumulation point $\mathbf{x}^*$ of $\{\mathbf{x}^{p*}\}_{p \in P}$ satisfies $\mathbf{x}^* \in \text{vert}([\mathbf{l}^*, \mathbf{u}^*])$.

**Proof.**

(i) This follows from the fact that $\mathbf{l}^p$ (resp. $\mathbf{u}^p$) is a bounded, non-decreasing (resp. non-increasing) sequence.

(ii) This follows from the Bolzano-Weierstrass theorem since the sequence $\{\mathbf{x}^{p*}\}_{p \in P}$ is generically bounded by $[\mathbf{l}, \mathbf{u}]$, $\bar{\mathcal{X}}$ is compact and $\text{conv}(S)$ is closed and bounded, thus compact (indeed, for all $i \in Q$, $t_i^p$ is trivially bounded by $\sup\{\phi_i(\mathbf{x}, \mathbf{y}) : \mathbf{x} \in \bar{\mathcal{X}}, \mathbf{y} \in \bar{\mathcal{Y}}(\mathbf{x})\}$ which is finite by Assumption 4.1).

(iii) Consider any accumulation point $\mathbf{x}^*$ of $\{\mathbf{x}^{p*}\}_{p \in P}$ with its associated convergent sub-sequence $P' \subseteq P$, i.e., $\{\mathbf{x}^{p*}\}_{p \in P'} \to \mathbf{x}^*$. Let $j \in I_C$ be as described in the lemma, and consider the sub-sequence $P^{u_j} \subseteq P'$ such that, for all $p \in P^{u_j}$, $u_j^{p+1} = x_j^{p*}$. Assume $P^{u_j}$ is not finite. Then, we have that $\{x_j^{p*}\}_{p \in P^{u_j}} \to x_j^*$ since $P^{u_j} \subseteq P'$. And thus, by definition of $P^{u_j}$, we have that $\{u_j^{p+1}\}_{p \in P^{u_j}} \to x_j^*$. We therefore chose $P^j = \{p+1 : p \in P^{u_j}\}$ and have $\{u_j^p\}_{p \in P^j} \to x_j^*$. If instead $P^{u_j}$ is finite, the sub-sequence $P^j \subseteq P'$ defined by nodes $p$ for which $l_j^{p+1} = x_j^{p*}$ is infinite; therefore, the similar argument can be applied.

(iv) We have just shown that, for any accumulation point $\mathbf{x}^*$ of $\{\mathbf{x}^{p*}\}_{p \in P}$, with its associated convergent sub-sequence

$P' \subseteq P$, and any infinitely branched index $j \in I_C$, there exists $P^j \subseteq P'$ such that either $\{u_j^p\}_{p \in P^j} \to x_j^*$ or $\{l_j^p\}_{p \in P^j} \to x_j^*$. Assume $\{l_j^p\}_{p \in P'} \to l_j^*$. Thus, $x_j^* = l_j^*$ holds, since any sub-sequence of a converging sequence converges to the same point. The same argument can be applied when $\{u_j^p\}_{p \in P^j} \to x_j^*$. □

**Theorem 1.** Let P be a sequence of nodes of any infinite branch of the branching tree. Then, every accumulation point of $\{\mathbf{x}^p, \mathbf{t}^p, \mathbf{y}^p\}_{p \in P}$, say $(\mathbf{x}^*, \mathbf{t}^*, \mathbf{y}^*)$, is an optimal solution of problem (7), and, thus, $\mathbf{x}^*$ is an optimal solution of (2SRO-P).

**Proof.** By Lemma 4 (ii), there exists a sub-sequence $P' \subseteq P$ such that $\{\mathbf{x}^{p*}, \mathbf{t}^{p*}, \mathbf{y}^{p*}\}_{p \in P'} \to (\mathbf{x}^*, \mathbf{t}^*, \mathbf{y}^*)$. Note that $\bar{\mathcal{X}}$ and $\text{conv}(S)$ are compact sets, hence we have that $(\mathbf{x}^*, \mathbf{t}^*, \mathbf{y}^*) \in \bar{\mathcal{X}} \times \text{conv}(S)$. Moreover, by Lemma 4 (iv), we know that $\mathbf{x}^* \in \text{vert}([\mathbf{l}^*, \mathbf{u}^*])$ which, by Proposition 2, ensures that $(\mathbf{t}^*, \mathbf{y}^*) \in \text{conv}(\mathcal{Y}'(\mathbf{x}^*))$. Hence, $(\mathbf{x}^*, \mathbf{t}^*, \mathbf{y}^*)$ is feasible for (7). Note that $f_{ST}$ is a continuous function since it is the point-wise supremum of continuous (affine) functions. Thus, by Remark 6, we have $\{f_{LB}(\mathbf{t}^{p*}, \mathbf{V}^{p*}, \xi^{p*})\}_{p \in P'} \to f_{ST}(\mathbf{t}^*)$. In other words, the objective value of the feasible solution $(\mathbf{x}^*, \mathbf{t}^*, \mathbf{y}^*)$ to (7) is $f_{ST}(\mathbf{t}^*)$. Yet, by Lemma 2, we know that (7) and (2SRO-P) have the same objective value. This makes $\mathbf{x}^*$ a feasible solution to (2SRO-P) of value $f_{ST}(\mathbf{t}^*)$. Since our node selection strategy always picks a node with minimum lower bound, for each node $p \in P$, we have $f_{LB}(\mathbf{t}^{p*}, \mathbf{V}^{p*}, \xi^{p*}) = \nu(LB^p) \le \nu(2SRO-P) \le f_{ST}(\mathbf{t}^*)$. As $\nu(LB^p)$ converges to $f_{ST}(\mathbf{t}^*)$, we also have $f_{LB}(\mathbf{t}^{p*}, \mathbf{V}^{p*}, \xi^{p*}) = f_{ST}(\mathbf{t}^*) = \nu(2SRO-P)$. □

We conclude this section by observing that, at each node of the branch-and-price algorithm, the lower bounding problem can be solved with $\varepsilon$-tolerance in a finite number of operations. Indeed, as shown in Ceria & Soares (1999) and Grossmann & Ruiz (2012), one can reformulate a convex disjunctive program as a convex MINLP by introducing an exponential number of auxiliary variables that model the disjunctions $\bigcup_{k \in \mathbb{Z} \cap [l_j, u_j]} \{x_j = k\}$ for each $j \in I_l$. The resulting model can thus be solved in finite number of steps by using any algorithm designed for convex optimization. However, relying on such a large reformulation is only acceptable for small problems with a limited number of integer variables. To avoid this, we introduce a nonlinear column-generation algorithm able to solve the lower bounding problem at each node by sequentially improving an inner approximation for the convex hull of $S^p$. This is presented in the next section.

### 3.4. A convexification scheme based on column-generation

In this section, we propose a nonlinear column-generation algorithm to be used, at each node $p$, to solve problem $(LB^p)$ to $\varepsilon$-optimality in a finite number of iterations. According to this scheme, we approximate $\text{conv}(S^p)$ by the convex hull of a finite set of points belonging to $S^p$.

**Restricted Master Problem:** To determine this set, we use an iterative approach. At each iteration $k$, let $K = \{1, \dots, k\}$ and de-
note by $H^{pk} = \{(\bar{\mathbf{x}}^{pj}, \bar{\mathbf{t}}^{pj}, \bar{\mathbf{y}}^{pj}) : j \in K\}$ the current set of points in the restricted master. As $H^{pk} \subseteq S^p$, we have $\mathrm{conv}(H^{pk}) \subseteq \mathrm{conv}(S^p)$, thus the optimal solution of the problem obtained by substituting $\mathrm{conv}(S^p)$ with $\mathrm{conv}(H^{pk})$ in $(LB^p)$ gives an upper bound of $(LB^p)$. The resulting problem, denoted as $(LB^{pk})$, is called the Restricted Master, and is formulated as follows.

$$ \min_{\mathbf{x}, \mathbf{t}, \mathbf{y}, \mathbf{v}, \mathbf{\xi}, \alpha} \delta^*(\mathbf{\xi} | \Xi) - \sum_{i \in Q} (t_i w_i)_* (\mathbf{v}^i) \quad (21) $$

$$ \text{subject to} \quad \mathbf{x} \in \bar{\mathcal{X}} \cap [\mathcal{I}^p, \mathcal{U}^p] \quad (22) $$

$$ \mathbf{x} = \sum_{j \in K} \alpha_j \bar{\mathbf{x}}^{pj} \quad (23) $$

$$ \mathbf{t} = \sum_{j \in K} \alpha_j \bar{\mathbf{t}}^{pj} \quad (24) $$

$$ \mathbf{y} = \sum_{j \in K} \alpha_j \bar{\mathbf{y}}^{pj} \quad (\widehat{LB}^{pk}) \quad (25) $$

$$ \sum_{j \in K} \alpha_j = 1 \quad (26) $$

$$ \sum_{i \in Q} \mathbf{v}^i = \mathbf{\xi} \quad (27) $$

$$ \mathbf{v}^i \in \mathbb{R}^{|\mathcal{I}|} \quad \forall i \in \mathcal{Q} \quad (28) $$

$$ \mathbf{\xi} \in \mathbb{R}^{|\mathcal{I}|} \quad (29) $$

$$ \alpha_j \ge 0 \quad \forall j \in K. \quad (30) $$

Following the classical column-generation framework, the current approximation can be improved by means of a so-called Pricing Problem, aimed at identifying new points to be added to the Restricted Master, and defined as follows.

**Pricing Problem:** Let $\lambda^{pk*}$, $\mu^{pk*}$, $\pi^{pk*}$ and $\eta^{pk*}$ be the values of the dual variables associated with constraints (23), (24), (25), and (26) in an optimal solution of problem $(LB^{pk})$.

Pricing asks to solve the following problem

$$ (\bar{\mathbf{x}}^{p,k+1}, \bar{\mathbf{t}}^{p,k+1}, \bar{\mathbf{y}}^{p,k+1}) \in \underset{(\mathbf{x}, \mathbf{t}, \mathbf{y}) \in S^p}{\operatorname{argmin}} - \lambda^{pk*} \mathbf{x} - \mu^{pk*} \mathbf{t} - \pi^{pk*} \mathbf{y} - \eta^{pk*} \quad (pp^{pk}) $$

and generates a new point $(\bar{\mathbf{x}}^{p,k+1}, \bar{\mathbf{t}}^{p,k+1}, \bar{\mathbf{y}}^{p,k+1})$ belonging to $S^p$. If $v(pp^{pk}) \ge -\epsilon$, we have an $\epsilon$-optimal solution to $(LB^p)$, and hence the algorithm terminates. Otherwise, we set $H^{k+1} = H^k \cup \{(\bar{\mathbf{x}}^{p,k+1}, \bar{\mathbf{t}}^{p,k+1}, \bar{\mathbf{y}}^{p,k+1})\}$, $k = k+1$ and iterate. Note that, at each iteration $k$, a lower bound on the optimal solution value of $(LB^p)$ is given by $v(\widehat{LB}^{pk}) - v(pp^{pk})$. This lower bound, combined with an upper bound, can allow us to early terminate the solution of problem $(LB^p)$.

# 4. Computational experiments

In this section, we report computational results of our solution algorithm when applied to an uncertain Capacitated Facility Location Problem with congestion.

## 4.1. Problem definition

We consider a variant of the Facility Location Problem, in which we are given a set $V_1$ of candidate sites for opening facilities, as well as a set $V_2$ of clients to be served with some product. Each client $j \in V_2$ has a demand $d_j$ representing the quantity of product that she/he wants to receive. Each site $i \in V_1$ can be activated at a given fixed cost $f_i > 0$. In this case, one has to decide the capacity to be installed, at cost $u_i$ per unit of capacity. Each site $i$ has an upper bound $\bar{q}_i$ on the maximum capacity that can be installed. Each connection $(i, j) \in V_1 \times V_2$ is associated with a fixed cost $c_{ij}$,

and a variable cost $t_{ij}$ per unit of product which is transported. In our setting, we explicitly model congestion at each site $i$ by means of an additional cost which depends on the total amount of product, say $o_i$, leaving the facility. As in the congested Facility Location Problem considered in Desrochers, Marcotte, & Stan (1995) and in Fischetti, Ljubić, & Sinnl (2016), the congestion cost for site $i$ is given by

$$ F_i(o_i) = (\alpha_i + \beta_i o_i^{\gamma_i}) o_i, \quad (31) $$

where $\alpha_i \ge 0$, $\beta_i > 0$ and $\gamma_i \ge 1$ are input parameters. Note that each function $F_i$ is convex for non-negative arguments $o_i$. The problem asks to determine the facilities to be opened, the capacity to be installed at each facility and the flow of product from facilities to clients, so as to serve all clients at minimum cost. This problem can be reduced to the one addressed in Desrochers et al. (1995) in case there are no capacity constraints at the sites (i.e., $\bar{q}_i = \infty$ and $u_i = 0$ for each $i \in V_1$) and transportation costs only include a variable component (i.e., $c_{ij} = 0$ for each $(i, j) \in V_1 \times V_2$).

In our context, connection costs are not known when deciding the capacities to be installed. Formally, for each connection $(i, j) \in V_1 \times V_2$, we denote by $\tilde{c}_{ij}$ and $\tilde{t}_{ij}$ the nominal fixed and variable costs from $i$ to $j$, and by $\tilde{c}_{ij}$ and $\tilde{t}_{ij}$ their maximal deviations. Without loss of generality, we assume that, for each connection $(i, j) \in V_1 \times V_2$, the actual realizations for the costs are determined by the same variable $\tilde{\xi}_{ij}$. In other words, we have $c_{ij} = \tilde{c}_{ij} + \tilde{\xi}_{ij}\tilde{t}_{ij}$ and $t_{ij} = \tilde{t}_{ij} + \tilde{\xi}_{ij}\tilde{t}_{ij}$, with $\tilde{\xi} \in \Xi$ and $\Xi \subseteq [0, 1]^{|V_1| \times |V_2|}$ is a given uncertainty set (see Section 4.2).

We consider the adjustable robust version of this uncertain problem, where capacity installation is determined at the first stage whereas product flows are determined after uncertainty reveals. We denote the resulting problem as ARCCFLP (for Adjustable Robust Congested Capacitated Facility Location Problem).

## 4.2. Mathematical formulation

To model ARCCFLP, we introduce, for each site $i \in V_1$, first-stage variables $x_i$ and $q_i$; the former takes the value 1 if site $i$ is activated, whereas the latter denotes the actual capacity installed. The feasible set $X$ for the first-stage variables is defined as

$$ X = \{ (\mathbf{x}, \mathbf{q}) : x_i \in \{0, 1\} \text{ and } 0 \le q_i \le \bar{q}_i x_i \quad \forall i \in V_1 \}. \quad (32) $$

Once the actual realization of uncertainty $\mathbf{\xi} \in \Xi$ is known, thus defining the transportation costs, the remaining decisions concerning the flow of product from opened sites to clients must be taken. To this aim we introduce, for each connection $(i, j) \in V_1 \times V_2$, variables $b_{ij}$ and $y_{ij}$ denoting if the connection is activated and the fraction of request of client $j$ that is served by site $i$, respectively. For each site $i$, we also denote by $o_i$ the total amount of product leaving the site. Accordingly, the feasible set $Y(\mathbf{x}, \mathbf{q})$ associated with a given pair $(\mathbf{x}, \mathbf{q})$ is defined by the following constraints.

$$ o_i = \sum_{j \in V_2} d_j y_{ij} \le q_i \quad \forall i \in V_1 \quad (33) $$

$$ \sum_{i \in V_1} y_{ij} = 1 \quad \forall j \in V_2 \quad (34) $$

$$ y_{ij} \le b_{ij} \quad \forall i \in V_1, \forall j \in V_2 \quad (35) $$

$$ y_{ij} \ge 0 \quad \forall i \in V_1, \forall j \in V_2 \quad (36) $$

$$ b_{ij} \in \{0, 1\} \quad \forall i \in V_1, \forall j \in V_2. \quad (37) $$

Constraints (33) enforce that the total demand $o_i$ leaving each site $i$ does not exceed the installed capacity, while constraints (34) impose that, for each client, all the demand is served. Constraints
(35) activate connections with a positive flow. Finally, (36) and (37) give the domain of the variables.

In order to ensure complete recourse, we introduce a dummy facility k with $f_k = u_k = 0$, $\bar{q}_k = \sum_{j \in V_2} d_j$, and with very large values for $\tilde{c}_{kj}$ and $\tilde{t}_{kj}$ for each $j \in V_2$. By adding constraints $x_k = 1$ and $q_k = \bar{q}_k$ in the definition of set X, we force facility k to be opened at maximum capacity. This choice has zero cost in the first stage, and allows a feasible solution with very large cost in the second stage regardless the values of the remaining first-stage variables.

Then, ARCCFLP is formulated as

$$ \min_{(\mathbf{x}, \mathbf{q}) \in X} \left\{ \sum_{i \in V_1} (f_i x_i + u_i q_i) + \max_{\xi \in \Xi} \min_{(\mathbf{b}, \mathbf{y}, \mathbf{o}) \in Y(\mathbf{x}, \mathbf{q})} \sum_{i \in V_1} \left( F_i(o_i) + \sum_{j \in V_2} \left( (\tilde{c}_{ij} + \xi_{ij} \tilde{c}_{ij}) b_{ij} + (\tilde{t}_{ij} + \xi_{ij} \tilde{t}_{ij}) d_j y_{ij} \right) \right) \right\}. \quad (38) $$

By applying the methodology introduced in this paper, the resulting lower-bounding problem is given as follows

$$ \min_{\mathbf{x}, \mathbf{q}, \mathbf{r}, \mathbf{o}, \mathbf{y}, \mathbf{b}} \sum_{i \in V_1} \left( f_i x_i + u_i q_i + r_i + \sum_{j \in V_2} (\tilde{c}_{ij} b_{ij} + \tilde{t}_{ij} d_j y_{ij}) \right) + \max_{\xi \in \Xi} \sum_{i \in V_1} \sum_{j \in V_2} \xi_{ij} (\tilde{c}_{ij} b_{ij} + \tilde{t}_{ij} d_j y_{ij}) \quad (39) $$

$$ \text{subject to } (\mathbf{x}, \mathbf{q}) \in X \quad (40) $$

$$ (\mathbf{q}, \mathbf{r}, \mathbf{o}, \mathbf{y}, \mathbf{b}) \in \text{conv} \left( \left\{ (\mathbf{q}, \mathbf{r}, \mathbf{o}, \mathbf{y}, \mathbf{b}) : \begin{array}{ll} 0 \le q_i \le \bar{q}_i & \forall i \in V_1 \\ r_i \ge \tilde{F}_i(o_i) & \forall i \in V_1 \\ (33) - (37) & \end{array} \right\} \right). \quad (41) $$

The inner maximization problem can then be expressed by using Fenchel duality, and the resulting formulation depends on the uncertain set. In our experiments, we consider two widely used uncertainty sets, namely, the $\Gamma$-uncertainty set and the ellipsoidal uncertainty set.

$\Gamma$-uncertainty is a paradigm introduced in Bertsimas & Sim (2004) to model uncertain situations in which the robustness of the solution can be controlled by an input parameter $\Gamma > 0$. By using this uncertainty set in our context, we obtain

$$ \Xi_{\Gamma}^{\text{d}} = \left\{ \xi \in [0, 1]^{|V_1| \times |V_2|} : \sum_{i \in V_1} \sum_{j \in V_2} \xi_{ij} \le \Gamma \right\}. \quad (42) $$

In this case, Fenchel duality reduces to LP duality as follows.

$$ \max_{\xi \in \Xi_{\Gamma}^{\text{d}}} \sum_{i \in V_1} \sum_{j \in V_2} \xi_{ij} (\tilde{c}_{ij} b_{ij} + \tilde{t}_{ij} d_j y_{ij}) = \min_{\Gamma} \Gamma \lambda + \sum_{i \in V_1} \sum_{j \in V_2} \pi_{ij} \quad (43) $$

$$ \text{subject to } \lambda + \pi_{ij} \ge \tilde{c}_{ij} b_{ij} + \tilde{t}_{ij} d_j y_{ij} \quad \forall i \in V_1, \forall j \in V_2 \quad (44) $$

$$ \lambda \ge 0 \quad (45) $$

$$ \pi_{ij} \ge 0 \quad \forall i \in V_1, \forall j \in V_2. \quad (46) $$

The Ellipsoidal uncertainty set is defined as

$$ \Xi_{\Gamma}^{\text{e}} = \left\{ \xi \in [0, 1]^{|V_1| \times |V_2|} : \sqrt{\sum_{i \in V_1} \sum_{j \in V_2} \xi_{ij}^2} \le \Gamma \right\}. \quad (47) $$

where again $\Gamma$ is a control parameter. In this case, one obtains the following formulation.

$$ \max_{\xi \in \Xi_{\Gamma}^{\text{e}}} \sum_{i \in V_1} \sum_{j \in V_2} \xi_{ij} (\tilde{c}_{ij} b_{ij} + \tilde{t}_{ij} d_j y_{ij}) = \min_{\Gamma} \Gamma \lambda + \sum_{i \in V_1} \sum_{j \in V_2} \pi_{ij} \quad (48) $$

$$ \text{subject to } v_{ij} + \pi_{ij} \ge \tilde{c}_{ij} b_{ij} + \tilde{t}_{ij} d_j y_{ij} \quad \forall i \in V_1, \forall j \in V_2 \quad (49) $$

$$ \sqrt{\sum_{i \in V_1} \sum_{j \in V_2} v_{ij}^2} \le \lambda \quad (50) $$

$$ \lambda \ge 0 \quad (51) $$

$$ \pi_{ij} \ge 0 \quad \forall i \in V_1, \forall j \in V_2 \quad (52) $$

$$ v_{ij} \ge 0 \quad \forall i \in V_1, \forall j \in V_2. \quad (53) $$

An interested reader may refer to Li, Ding, & Floudas (2011) for associated theoretical properties of both uncertainty sets, including their robust counterparts and probabilistic guarantees for linear constraints.

### 4.3. Test bed

**Instance generation** We tested our solution method on random instances, that were generated by following the guidelines of the extensive computational study by Cornuejols, Sridharan, & Thizy (1991). Accordingly, for each facility $i \in V_1$, the maximum capacity $\bar{q}_i$ and the fixed opening cost $f_i$ follow uniform distributions in [10, 160] and [0, 180], respectively, whereas the variable coefficient $u_i$ was generated in $[200/\sqrt{\bar{q}_i}, 220/\sqrt{\bar{q}_i}]$. Moreover, locations for each facility $i \in V_1$ and each client $j \in V_2$ were generated in the unit square, and nominal transportation costs were set to the Euclidean distance multiplied by 10 and rounded up. The demands were uniformly generated between 0 and 1 and scaled so that $\sum_{i \in V_1} \bar{q}_i / \sum_{j \in V_2} d_j = \nu$ where $\nu$ is a parameter taking value in $\{1.1, 1.2, 1.3\}$. In addition, following Desrochers et al. (1995), for each $i \in V_1$, we used $\gamma_i = 1$ and $\alpha_i = \beta_i = 0.75$, i.e., each function $F_i$ is quadratic with respect to the amount of product leaving site i. Concerning the parameters affected by uncertainty, the maximum deviation $\tilde{t}_{ij}$ was set to $0.50 \times \tilde{t}_{ij}$, thus allowing a maximum of 50% deviation. As to the opening cost of each arc, we randomly generated the nominal value between 50 and 100, allowing a maximum of 50% deviation with respect to this value.

Finally, the number of sites and clients take values (4, 8), (5, 10) and (6, 12). For each combination of $|V_1|$, $|V_2|$, and $\nu$, we generated 5 test-cases. Each test-case was solved for $\Gamma = 1, 2, 3, 4$, both in the $\Gamma$-uncertainty and in the Ellipsoidal uncertainty settings, thus producing 360 instances.

### 4.4. Implementation details

We implemented our branch-and-price algorithm with spatial branching in C++ and run all the experiments on an AMD Ryzen 5 PRO 4650GE at 3.3 GHz, with a time limit equal to 10,800 CPU seconds per run (3 hours).

At the root node, an initial upper bound is computed by solving the single-stage version of ARCCFLP where all decisions are taken here and now. This bound is obtained by solving (39)–(41) without convexifying the second-stage feasible region in constraint (41). At each node of the algorithm, we solve the restricted master problem by using Mosek 10.0.36, and the pricing problem by means of Gurobi 10.0. This combination of solvers turned out to be the most numerically stable on our instances. Our code and instances are publicly available on GitHub.<fnref n="1" />

<footnote n="1">¹ https://github.com/hlefebvr/AB_AdjustableRobustOptimizationWithObjectiveUncertainty.</footnote>
<table><caption>**Table 1**
Computational experiments on ARCCFLP instances. Each row refers to 15 instances.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">|V<sub>1</sub>|</th><th rowspan="2">|V<sub>2</sub>|</th><th rowspan="2">Γ</th><th rowspan="2">solved</th><th colspan="4">Time (s)</th><th colspan="2">Gap (%)</th></tr><tr><th>total</th><th>RMP</th><th>PP</th><th>nodes</th><th>columns</th><th>root</th><th>end</th></tr></thead><tbody><tr><td rowspan="12">Γ-unc.</td><td rowspan="4">4</td><td rowspan="4">8</td><td>1</td><td>15</td><td>2.45</td><td>0.12</td><td>2.04</td><td>4.73</td><td>42.60</td><td>0.21</td><td></td></tr><tr><td>2</td><td>15</td><td>3.19</td><td>0.14</td><td>2.68</td><td>5.00</td><td>56.53</td><td>0.21</td><td></td></tr><tr><td>3</td><td>15</td><td>4.75</td><td>0.26</td><td>4.03</td><td>5.27</td><td>84.60</td><td>0.22</td><td></td></tr><tr><td>4</td><td>15</td><td>4.74</td><td>0.22</td><td>4.10</td><td>5.27</td><td>78.60</td><td>0.23</td><td></td></tr><tr><td rowspan="4">5</td><td rowspan="4">10</td><td>1</td><td>15</td><td>44.83</td><td>0.78</td><td>42.35</td><td>5.80</td><td>110.00</td><td>0.29</td><td></td></tr><tr><td>2</td><td>15</td><td>21.83</td><td>0.47</td><td>20.27</td><td>6.60</td><td>82.20</td><td>0.31</td><td></td></tr><tr><td>3</td><td>15</td><td>36.59</td><td>0.64</td><td>34.70</td><td>7.53</td><td>125.87</td><td>0.33</td><td></td></tr><tr><td>4</td><td>15</td><td>38.33</td><td>0.80</td><td>36.16</td><td>9.40</td><td>166.80</td><td>0.34</td><td></td></tr><tr><td rowspan="4">6</td><td rowspan="4">12</td><td>1</td><td>11</td><td>319.17</td><td>1.92</td><td>313.31</td><td>4.82</td><td>139.55</td><td>0.17</td><td>0.24</td></tr><tr><td>2</td><td>13</td><td>572.80</td><td>1.64</td><td>562.48</td><td>5.46</td><td>102.69</td><td>0.19</td><td>0.30</td></tr><tr><td>3</td><td>13</td><td>1117.35</td><td>2.35</td><td>1110.73</td><td>6.23</td><td>162.92</td><td>0.22</td><td>0.23</td></tr><tr><td>4</td><td>11</td><td>1261.83</td><td>4.63</td><td>1249.46</td><td>6.45</td><td>239.00</td><td>0.22</td><td>0.26</td></tr><tr><td rowspan="12">Ellips.</td><td rowspan="4">4</td><td rowspan="4">8</td><td>1</td><td>13 (2)</td><td>6.64</td><td>0.35</td><td>5.70</td><td>5.00</td><td>67.92</td><td>0.22</td><td>0.11</td></tr><tr><td>2</td><td>14 (1)</td><td>7.46</td><td>0.52</td><td>6.25</td><td>5.00</td><td>88.79</td><td>0.22</td><td>0.30</td></tr><tr><td>3</td><td>14 (1)</td><td>12.47</td><td>0.67</td><td>10.90</td><td>5.29</td><td>121.07</td><td>0.20</td><td>0.46</td></tr><tr><td>4</td><td>13 (2)</td><td>11.16</td><td>0.40</td><td>10.08</td><td>4.54</td><td>86.77</td><td>0.21</td><td>0.09</td></tr><tr><td rowspan="4">5</td><td rowspan="4">10</td><td>1</td><td>12 (3)</td><td>26.96</td><td>1.27</td><td>24.59</td><td>6.00</td><td>117.58</td><td>0.26</td><td>0.44</td></tr><tr><td>2</td><td>12 (3)</td><td>83.59</td><td>2.34</td><td>79.28</td><td>7.17</td><td>206.75</td><td>0.32</td><td>0.38</td></tr><tr><td>3</td><td>15</td><td>82.47</td><td>3.19</td><td>76.72</td><td>8.73</td><td>278.80</td><td>0.34</td><td></td></tr><tr><td>4</td><td>15</td><td>151.89</td><td>2.17</td><td>147.62</td><td>6.33</td><td>195.27</td><td>0.31</td><td></td></tr><tr><td rowspan="4">6</td><td rowspan="4">12</td><td>1</td><td>14 (1)</td><td>324.66</td><td>3.45</td><td>315.11</td><td>5.29</td><td>169.57</td><td>0.18</td><td>0.21</td></tr><tr><td>2</td><td>14 (1)</td><td>420.62</td><td>4.01</td><td>412.94</td><td>5.86</td><td>202.21</td><td>0.22</td><td>0.21</td></tr><tr><td>3</td><td>13</td><td>859.62</td><td>5.66</td><td>847.91</td><td>7.00</td><td>280.08</td><td>0.24</td><td>0.10</td></tr><tr><td>4</td><td>13 (2)</td><td>1574.26</td><td>10.83</td><td>1551.17</td><td>10.69</td><td>526.38</td><td>0.22</td><td>0.15</td></tr></tbody></table>

The column-generation procedure includes stabilization by dual price smoothing, as described in Pessoa, Sadykov, Uchoa, & Vanderbeck (2013); and at most one column is added to the restricted master problem at each iteration.

For the branching strategy of continuous variables, we used a tolerance of $\epsilon = 10^{-3}$ for comparing real numbers in finite precision. At each node, local bounds derived from branching decisions are applied to the column generation sub-problem. The restricted master inherits all columns from the father node. For each column, we compare the value of the first-stage variables in the column with the actual bound at the node, and possibly remove it from the master when locally infeasible.

At a given node, to check the optimality of a first-stage decision and possibly fathom the node, we use the sufficient condition from Proposition 2. If the latter does not hold, we check if all active columns at optimality are built on the same values for the continuous variables; in this case, Proposition 2 can be exploited to ensure local optimality of the solution. If the first-stage solution is not feasible, an upper bound is computed as follows. We detect the variable with largest value in the RMP current solution, recover the values of variables $(x, q)$ that were used for generating this column and fix variables $(x, q)$ to those values, possibly rounding up integer variables.

Finally, observe that branching may induce infeasibility in the second stage. To early detect this situation, we consider a problem-specific improvement: at a given node of the branching tree, we check if $\sum_{i \in V_1} x_i^l q_i^l < \sum_{j \in V_2} d_j$ holds, where $x_i^l$ and $q_i^l$ denote the local upper bound of variable $x_i$ and $q_i$, respectively. In this case the node is declared infeasible.

## 4.5. General results

Table 1 reports our computational results on ARCCFLP. The upper part relates to experiments done with the $\Gamma$-uncertainty set ($\Gamma$-unc.), while the lower part addresses those with the ellipsoidal uncertainty set (Ellips.). Columns $|V_1|$, $|V_2|$ and $\Gamma$ give the number of sites, the number of clients and the value for the uncertainty parameter $\Gamma$, respectively. Column “solved” reports the number of

instances (out of 15) which could be solved to proven optimality within the given time limit. Into brackets we report the number of instances for which the computation was stopped due to numerical issues of the used solvers. For the sake of consistency, all remaining columns but the last one are relative to instances which could be solved within the time limit. In particular, columns “time” report, from left to right, the average time needed to prove optimality (“total”), the average time spent solving the RMP (“RMP”) and the average time spent solving the pricing problem (“PP”) during the execution of our branch-and-price algorithm. All times are expressed in seconds. Column “nodes” reports the average number of explored nodes, while “columns” gives the average number of generated columns throughout the entire execution of our algorithm. Finally, we report the average gap of the root node lower bound with respect to the optimal solution value (“root”) and the average optimality gap at time limit (“end”) (or when the algorithm was stopped for numerical troubles). Optimality gap is computed only over those instances which could not be solved to optimality.

The table shows that our method is able to solve a large fraction of the instances, namely 168 in the $\Gamma$-uncertainty setting and 162 in the Ellipsoidal uncertainty setting. In most cases, the solution time required by the algorithm is quite small and solving the RMP is very fast in practice (below 7% of the total time, on average). Indeed, the most challenging subproblem solved is the pricing problem; when increasing the size of the instances the number of columns that are needed increases and each execution of the pricing problem is more time consuming. However, the solved relaxation allows to compute a very tight approximation, as the gap between lower and upper bounds at the root that is always below 0.35%, and producing small enumeration trees, in which the number of generated nodes is always below 11. Moreover, the performance of the algorithm is satisfactory also for the instances that were not solved to proven optimality, as the average residual gap at the end of the enumeration is always quite small (below 0.5%). Finally, observe that numerical issues arise only when uncertainty belongs to the Ellipsoidal uncertainty set, a nonlinear setting in which Mosek may encounter numerical instability on some instances.
<table><caption>**Table 2**
Comparison of exact and linearized approaches.</caption><thead><tr><th rowspan="2"></th><th rowspan="2">|V<sub>1</sub>|</th><th rowspan="2">|V<sub>2</sub>|</th><th rowspan="2">exact<br>solved</th><th colspan="3">linearized</th></tr><tr><th>solved</th><th>avg. err.</th><th>max. err.</th></tr></thead><tbody><tr><td rowspan="3">Γ-unc.</td><td>4</td><td>8</td><td>60</td><td>60</td><td>0.37</td><td>0.61</td></tr><tr><td>5</td><td>10</td><td>60</td><td>60</td><td>0.38</td><td>0.63</td></tr><tr><td>6</td><td>12</td><td>48</td><td>42</td><td>0.39</td><td>0.56</td></tr><tr><td rowspan="3">Ellips.</td><td>4</td><td>8</td><td>54</td><td>58</td><td>0.36</td><td>0.59</td></tr><tr><td>5</td><td>10</td><td>54</td><td>59</td><td>0.39</td><td>0.58</td></tr><tr><td>6</td><td>12</td><td>54</td><td>59</td><td>0.40</td><td>0.56</td></tr></tbody></table>

### 4.5.1. Linearized costs

A classical approach for dealing with non-linearities in optimization problems is to approximate them via piecewise linear functions. Thus, as an alternative approach for both uncertainty sets described in Section 4.2, we considered solving a linearized approximation of ARCCFLP obtained by replacing each function $F_i$ ($i \in V_1$) by a piecewise linear approximation. By introducing L approximation points {$\bar{o}_{il}$}$_{l=1,...,L}$, function $F_i$ is underestimated by the following one:

$$ \bar{F}_i(o_i) = \max_{l=1,...,L} \{F_i(\bar{o}_{il}) + F'_i(\bar{o}_{il})(o_i - \bar{o}_{il})\}. \quad (54) $$

In our experiments, we chose L = 10 and, for all $i \in V_1$, defined the approximation points to be equally distributed in the interval $[0, \bar{q}_i]$, i.e., $\bar{o}_{il} = \bar{q}_i(l-1)/(L-1)$ for $l = 1, ..., L$.

Table 2 reports the results obtained by using the linearized approach. For each combination of $|V_1|$ and $|V_2|$, we give for both the exact and the linearized approaches, the number of instances solved to optimality. Moreover, for the latter we report the average and maximum percentage error introduced by the linearization, computed as $\frac{z^*-z^L}{z^*}$, where $z^*$ and $z^L$ denote the optimal values of an ARCCFLP instance and of its linearized counterpart, respectively. These figures are computed with respect to instances solved by both approaches only.

The table shows that the linearized approach turns out to be computationally harder in the Γ-uncertainty setting, while it gives some improvement when the Ellipsoidal setting is considered; this is mainly due to the reduced number of instances for which we encountered numerical troubles. However, in both settings, linearization introduces a nonnegligible error when underestimating the true cost of a solution. The average percentage error, over all instances, is 0.38% and can be as large as 0.63%.

## 5. Conclusion

In this work, we studied optimization problems where part of the cost parameters are not known at decision time, and the decision flow is modeled as a two-stage process. In particular, we addressed general problems in which all constraints (including those linking the first and the second stages) are defined by convex functions and involve mixed-integer variables. To the best of our knowledge, this is the first attempt to extend the existing literature to tackle this wide class of problems.

To this aim, we derive a relaxation of the problem which can be formulated as a convex optimization problem, and embed it within an enumerative algorithm where branching occurs on integer and continuous variables. By combining enumeration and on-the-fly generation of the variables, we obtain a branch-and-price scheme, for which we prove asymptotic convergence in the general mixed-integer case and give sufficient conditions for finite convergence.

In addition to the theoretical analysis, we applied our method to an optimization problem affected by objective uncertainty arising in the logistic field, namely a variant of the congested Capacitated Facility Location problem with uncertain transportation costs.

Our computational experiments showed that the proposed method is able to solve non-trivial instances for this problem. In addition, we provide a comparison with a natural approach based on linearization of the congestion function, showing that this alternative solution method would give marginal improvements in terms of performances though introducing a nonnegligible error in terms of cost of the provided solution.

## Funding

This research was supported by "Mixed-Integer Non Linear Optimization: Algorithms and Application" consortium, which has received funding from the European Union's EU Framework Programme for Research and Innovation Horizon 2020 under the Marie Skłodowska-Curie Actions Grant Agreement No 764759.

## Acknowledgments

The authors are grateful to three anonymous reviewers for their constructive comments and remarks.

## Appendix A. Recalls of convex and concave conjugate

In this appendix we review some basic results on conjugate functions and Fenchel duality. For a detailed treatment we refer to Rockafellar (1970).

Let $f: \mathbb{R}^n \to \mathbb{R}$ be a given function, its convex conjugate is denoted by $f^*: \mathbb{R}^n \to \mathbb{R}$ and is given by

$$ f^*(\pi) = \sup_{\mathbf{x} \in \text{dom}(f)} \{\pi^T \mathbf{x} - f(\mathbf{x})\}. $$

Similarly, we denote by $g_*$ the concave conjugate of a given function $g: \mathbb{R}^n \to \mathbb{R}$, given by

$$ g_*(\pi) = \inf_{\mathbf{x} \in \text{dom}(g)} \{\pi^T \mathbf{x} - g(\mathbf{x})\}. $$

Note that, if $f$ is a proper convex function and $g$ a proper concave function, we have that $f^{**} = f$ and $g^{**} = g$. We now state the following Fenchel duality theorem.

**Theorem 2.** Let $f: \mathbb{R}^n \to \mathbb{R}$ be a proper convex function and $g: \mathbb{R}^n \to \mathbb{R}$ be a proper concave function, then

$$ \inf_{\mathbf{x} \in \text{dom}(f) \cap \text{dom}(g)} \{f(\mathbf{x}) - g(\mathbf{x})\} = \sup_{\pi \in \text{dom}(f^*) \cap \text{dom}(g_*)} \{g_*(\pi) - f^*(\pi)\} $$

or equivalently,

$$ \sup_{\mathbf{x} \in \text{dom}(f) \cap \text{dom}(g)} \{g(\mathbf{x}) - f(\mathbf{x})\} = \inf_{\pi \in \text{dom}(g_*) \cap \text{dom}(f^*)} \{f^*(\pi) - g_*(\pi)\}. $$

**Corollary 3** (Maximizing a concave function over a convex set). Let $\mathcal{X} \subseteq \mathbb{R}^n$ be a non-empty convex set, $g: \mathbb{R}^n \to \mathbb{R}$ be a proper concave function, then

$$ \sup_{\mathbf{x} \in \mathcal{X}} g(\mathbf{x}) = \inf_{\pi} \{\delta^*(\pi | \mathcal{X}) - g_*(\pi)\}, $$

$$ \text{where } \delta(\mathbf{x}|\mathcal{X}) = \begin{cases} 0 & \text{if } \mathbf{x} \in \mathcal{X} \\ +\infty & \text{otherwise.} \end{cases} $$

**Proof.** The result holds from the fact that $\sup\{g(\mathbf{x}) : \mathbf{x} \in \mathcal{X}\} = \sup\{g(\mathbf{x}) - \delta(\mathbf{x}|\mathcal{X})\}$ and by application of Fenchel duality. More precisely, $\delta(\mathbf{x}|\mathcal{X})$ is convex and, by non-emptiness of $\mathcal{X}$, is proper. $\square$

Notice that Fenchel duality allows the reformulation of an optimization problem which consists in maximizing a concave function over a convex set as an unconstrained convex problem since $\delta^*(\cdot|\mathcal{X})$ and $(-g_*(\cdot)|\mathcal{X})$ are convex functions and positively weighted sums preserve convexity.
<table><caption>**Table A.1**
Some convex conjugate calculus rules.</caption><thead><tr><th>h(x)</th><th>h<sup>*</sup>(π)</th></tr></thead><tbody><tr><td>Separable sums<br>h(x<sup>1</sup>, x<sup>2</sup>) = f<sub>1</sub>(x<sup>1</sup>) + f<sub>2</sub>(x<sup>2</sup>)</td><td>h<sup>*</sup>(π<sup>1</sup>, π<sup>2</sup>) = f<sub>1</sub><sup>*</sup>(π<sup>1</sup>) + f<sub>2</sub><sup>*</sup>(π<sup>2</sup>)</td></tr><tr><td>Scalar multiplications (α > 0)<br>h(x) = αf(x)</td><td>h<sup>*</sup>(π) = αf<sup>*</sup>(π/α)</td></tr><tr><td>Affine mapping composition (det A ≠ 0)<br>h(x) = f(Ax + b)</td><td>h<sup>*</sup>(π) = f<sup>*</sup>(A<sup>-T</sup>π) - b<sup>T</sup>A<sup>-T</sup>π</td></tr><tr><td>Sum with affine functions<br>h(x) = f(x) + a<sup>T</sup>x + b</td><td>h<sup>*</sup>(π) = f<sup>*</sup>(π - a) - b</td></tr><tr><td>Sum of functions<br>h(x) = ∑<sub>i=1</sub><sup>m</sup> f<sub>i</sub>(x)</td><td>h<sup>*</sup>(π) = inf<sub>x</sub> {∑<sub>i=1,...,m</sub><sup>m</sup> f<sub>i</sub><sup>*</sup>(x') | ∑<sub>i=1</sub><sup>m</sup> x' = π}</td></tr></tbody></table>

**Proposition 3.** Let $f$ be a convex function, we have $(-f)_*(\pi) = -f^*(-\pi)$.

**Proof.**

$$(-f)_*(\pi) = \inf_{x} \{ \pi^T x - (-f)(x) \} = -\sup_{x} \{-\pi^T x - f(x)\} \\ = -f^*(-\pi). \quad \square$$

Table A.1 reports some calculus rules regarding convex conjugates. The extension to concave conjugates is straightforward.

**References**

Arslan, A. N., & Detienne, B. (2022). Decomposition-based approaches for a class of two-stage robust binary optimization problems. *INFORMS Journal on Computing*, 34(2), 857–871.

Ayoub, J., & Poss, M. (2016). Decomposition for adjustable robust linear optimization subject to uncertainty polytope. *Computational Management Science*, 13, 219–239.

Ben-Tal, A., Ghaoui, L. E., & Nemirovski, A. (2009). *Robust optimization*. Princeton University Press.

Ben-Tal, A., Goryashko, A., Guslitzer, E., et al., (2004). Adjustable robust solutions of uncertain linear programs. *Mathematical Programming*, 99(2), 351–376.

Bertsimas, D., Brown, D., & Caramanis, C. (2010). Theory and applications of robust optimization. *SIAM Review*, 53, 464–501.

Bertsimas, D., & Caramanis, C. (2010). Finite adaptability in multistage linear optimization. *IEEE Transactions on Automatic Control*, 55(12), 2751–2766.

Bertsimas, D., Litvinov, E., Sun, X., et al., (2013). Adaptive robust optimization for the security constrained unit commitment problem. *IEEE Transactions on Power Systems*, 28, 52–63.

Bertsimas, D., & Sim, M. (2004). The price of robustness. *Operations Research*, 52(1), 35–53.

Ceria, S., & Soares, J. (1999). Convex programming for disjunctive convex optimization. *Mathematical Programming*, 86, 595–614.

Cornuejols, G., Sridharan, R., & Thizy, J. (1991). A comparison of heuristics and relaxations for the capacitated plant location problem. *European Journal of Operational Research*, 50(3), 280–297.

Desrochers, M., Marcotte, P., & Stan, M. (1995). The congested facility location problem. *Location Science*, 3(1), 9–23.

Fischetti, M., Ljubić, I., & Sinnl, M. (2016). Benders decomposition without separability: A computational study for capacitated facility location problems. *European Journal of Operational Research*, 253(3), 557–569.

Gabrel, V., Lacroix, M., Murat, C., & Remli, N. (2011). Robust location transportation problems under uncertain demands. *Discrete Applied Mathematics*, 164, 100–111.

Grossmann, I. E., & Ruiz, J. P. (2012). Generalized disjunctive programming: A framework for formulation and alternative algorithms for MINLP optimization. In J. Lee, & S. Leyffer (Eds.), *Mixed integer nonlinear programming* (pp. 93–115). New York, NY: Springer.

Hanasusanto, G. A., Kuhn, D., & Wiesemann, W. (2015). K-adaptability in two-stage robust binary programming. *Operations Research*, 63(4), 877–891.

Hassene, A., Bazgan, C., & Vanderpooten, D. (2009). Min-max and min-max regret versions of combinatorial optimization problems: A survey. *European Journal of Operational Research*, 197, 427–438.

Jiang, R., Zhang, M., Li, G., & Guan, Y. (2014). Two-stage network constrained robust unit commitment problem. *European Journal of Operational Research*, 234, 751–762.

Kämmerling, N., & Kurtz, J. (2020). Oracle-based algorithms for binary two-stage robust optimization. *Computational Optimization and Applications*, 77(2), 539–569.

Leyffer, S., Menickelly, M., Munson, T., Vanaret, C., & Wild, S. M. (2020). A survey of nonlinear robust optimization. *INFOR: Information Systems and Operational Research*, 58(2), 342–373.

Li, C., & Grossmann, I. E. (2019). A finite ε-convergence algorithm for two-stage stochastic convex nonlinear programs with mixed-binary first and second-stage variables. *Journal of Global Optimization*, 75(4), 921–947.

Li, Z., Ding, R., & Floudas, C. A. (2011). A comparative theoretical and computational study on robust counterpart optimization: I. robust linear optimization and robust mixed integer linear optimization. *Industrial & Engineering Chemistry Research*, 50(18), 10567–10603.

Perchet, V., & Vigeral, G. (2015). A minmax theorem for concave-convex mappings with no regularity assumptions. *Journal of Convex Analysis*, 22, 537–540.

Pessoa, A., Sadykov, R., Uchoa, E., & Vanderbeck, F. (2013). In-out separation and column generation stabilization by dual price smoothing. *SEA 2013*. In *Lecture notes in computes science: Vol. 7933* (pp. 354–365). Springer Berlin Heidelberg.

Rockafellar, R. (1970). *Convex analysis*. Princeton University Press.

Sherali, H. D., & Fraticelli, B. (2002). A modification of benders' decomposition algorithm for discrete subproblems: An approach for stochastic programs with integer recourse. *Journal of Global Optimization*, 22, 319–342.

Sherali, H. D., & Zhu, X. (2006). On solving discrete two-stage stochastic programs having mixed-integer first- and second-stage variables. *Mathematical Programming*, 108(2), 597–616.

Soyster, A. L. (1973). Technical note-convex programming with set-inclusive constraints and applications to inexact linear programming. *Operations Research*, 21(5), 1154–1157.

Subramanyam, A., Gounaris, C. E., & Wiesemann, W. (2019). K-adaptability in two-stage mixed-integer robust optimization. *Mathematical Programming Computation*, 12(2), 193–224.

Terry, T., Epelman, M., & Thiele, A. (2009). Robust linear optimization with recourse. *Technical Report*. Lehigh University.

Yanikoğlu, I., Gorissen, B. L., & den Hertog, D. (2019). A survey of adjustable robust optimization. *European Journal of Operational Research*, 277(3), 799–813.

Zeng, B., & Zhao, L. (2013). Solving two-stage robust optimization problems by a column-and-constraint generation method. *Operations Research Letters*, 41, 457–461.