# Two-stage and Lagrangian Dual Decision Rules for Multistage Adaptive Robust Optimization

Maryam Daryalal, Ayşe N Arslan, Merve Bodur

► To cite this version:

Maryam Daryalal, Ayşe N Arslan, Merve Bodur. Two-stage and Lagrangian Dual Decision Rules for Multistage Adaptive Robust Optimization. 2023. hal-04090602

HAL Id: hal-04090602

https://hal.science/hal-04090602

Preprint submitted on 5 May 2023

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_1.png" /></figure>

HAL
open science

# Two-stage and Lagrangian Dual Decision Rules for Multistage Adaptive Robust Optimization

Maryam Daryalal, Ayşe N Arslan, Merve Bodur

► To cite this version:

Maryam Daryalal, Ayşe N Arslan, Merve Bodur. Two-stage and Lagrangian Dual Decision Rules for Multistage Adaptive Robust Optimization. 2024. hal-04090602v2

HAL Id: hal-04090602

https://hal.science/hal-04090602v2

Preprint submitted on 12 Sep 2024

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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Two-stage and Lagrangian Dual Decision Rules for Multistage Adaptive Robust Optimization

Maryam Daryalal

Department of Decision Sciences, HEC Montréal, Montréal, Québec H3T 2A7, Canada,
maryam.daryalal@hec.ca

Ayşe N. Arslan

Univ. Bordeaux, CNRS, INRIA, Bordeaux INP, IMB, UMR 5251, F-33400 Talence, France,
Centre Inria de l'Université de Bordeaux, F-33405 Talence, France,
ayse-nur.arslan@inria.fr

Merve Bodur

Department of Mechanical and Industrial Engineering, University of Toronto, Toronto, Ontario M5S 3G8, Canada,
bodur@mie.utoronto.ca

In this work, we design primal and dual bounding methods for multistage adjustable robust optimization (MSARO) problems by adapting two decision rules rooted in the stochastic programming literature. This approach approximates the primal and dual formulations of an MSARO problem with two-stage models. From the primal perspective, this is achieved by applying two-stage decision rules that restrict the functional forms of a certain subset of decision variables. We present sufficient conditions under which the well-known constraint-and-column generation algorithm can be used to solve the primal approximation with finite convergence guarantees. From the dual side, we introduce a distributionally robust dual problem for MSARO models using their nonanticipative Lagrangian dual and then apply linear decision rules on the Lagrangian multipliers. For this dual approximation, we present a monolithic bilinear program valid for continuous recourse problems, and a cutting-plane method for mixed-integer recourse problems. Our framework is general-purpose and does not require strong assumptions such as a stage-wise independent uncertainty set, and can consider integer recourse variables. Computational experiments on newsvendor, location-transportation, and capital budgeting problems show that our bounds yield considerably smaller optimality gaps compared to the existing methods.

**Key words:** Optimization under uncertainty, Robust optimization, Decision rules

## 1. Introduction

Many practical planning, design and operational problems involve making decisions under uncertainty at consecutive stages, where the decisions in one stage affect the decisions of the future stages. In such sequential decision-making problems, first-stage (here-and-now) decisions are the
Two-stage and Lagrangian Dual Decision Rules for
Multistage Adaptive Robust Optimization

Maryam Daryalal

Department of Decision Sciences, HEC Montréal, Montréal, Québec H3T 2A7, Canada,
maryam.daryalal@hec.ca

Ayşe N. Arslan

Univ. Bordeaux, CNRS, INRIA, Bordeaux INP, IMB, UMR 5251, F-33400 Talence, France,
Centre Inria de l'Université de Bordeaux, F-33405 Talence, France,
ayse-nur.arslan@inria.fr

Merve Bodur

School of Mathematics, University of Edinburgh, Edinburgh EH9 3FD, UK,
merve.bodur@ed.ac.uk

In this work, we design primal and dual bounding methods for multistage adaptive robust optimization (MSARO) problems motivated by two decision rules rooted in the stochastic programming literature. From the primal perspective, this is achieved by applying decision rules that restrict the functional forms of only a certain subset of decision variables resulting in an approximation of MSARO as a two-stage adjustable robust optimization problem. We leverage the two-stage robust optimization literature in the solution of this approximation. From the dual perspective, decision rules are applied to the Lagrangian multipliers of a Lagrangian dual of MSARO, resulting in a two-stage stochastic optimization problem. As the quality of the resulting dual bound depends on the distribution chosen when developing the dual formulation, we define a distribution optimization problem with the aim of optimizing the obtained bound and develop solution methods tailored to the nature of the recourse variables. Our framework is general-purpose and does not require strong assumptions such as a stage-wise independent uncertainty set, and can consider integer recourse variables. Computational experiments on newsvendor, location-transportation, and capital budgeting problems show that our bounds yield considerably smaller optimality gaps compared to the existing methods.

Key words: Optimization under uncertainty, Robust optimization, Decision rules

# 1. Introduction

Many practical planning, design and operational problems involve making decisions under uncer-
tainty at consecutive stages, where the decisions in one stage affect the decisions of the future
stages. In such sequential decision-making problems, first-stage (here-and-now) decisions are the
ones that are immediately implementable. Subsequent recourse (wait-and-see) decisions depend on the state of the system, which is a result of previous decisions and observations of the uncertain parameters. A solution is then an adaptable policy or decision rule that takes the previous decisions and history of uncertainty realizations as an input, and returns a new implementable decision. The dynamics of a sequential decision-making problem is depicted in Figure 1.

<figure><img src="image_2.png" /><figcaption>**Figure 1** Sequential decision-making under uncertainty</figcaption></figure>

There are several modeling frameworks for sequential decision-making problems under uncertainty. With a known probability distribution, these problems may be addressed by the *multistage stochastic programming* (MSP) paradigm, with the goal of optimizing some statistical performance measure over the planning horizon. There is an extensive body of research on MSP problems of various structures, with a rich literature on problems with continuous decision variables. Without a precise distribution, Bertsimas and Thiele (2006) point out that an *assumed* probability distribution might heavily impact the results if it is different than the *true* one, even if both distributions share the same first and second moments. To mitigate this effect, *distributionally robust optimization* (DRO) models are proposed for making decisions that are based on a family of probability distributions, often defined by using historical data (Goh and Sim 2010). These models aim to hedge against tuning decisions to a perceived distribution. However, while some recent studies have proposed tractable solution methods for linear DRO problems under certain conditions (Philpott et al. 2018, Bertsimas et al. 2019), they remain largely challenging to solve. To the best of our knowledge, the only such methodology for a multistage mixed-integer problem is developed under the assumption of endogenous uncertainty (Yu and Shen 2020), meaning that the decisions impact the underlying stochastic structure of the problem, an assumption that might not always hold.

*Multistage adaptive robust optimization* (MSARO) is another framework for modeling sequential decision-making problems under uncertainty, when we do not have any knowledge about the underlying probability distribution or when the underlying uncertainty is not stochastic in nature. In designing such robust solutions, we protect ourselves against worst-case outcomes. While *static* (i.e., single-stage) robust optimization is tractable for a considerable number of problem structures (e.g. a linear problem with ellipsoidal uncertainty set, yielding a conic quadratic program (Ben-Tal
ones that are immediately implementable. Subsequent recourse (wait-and-see) decisions depend on the state of the system, which is a result of previous decisions and observations of the uncertain parameters. A solution is then an adaptable policy or decision rule that takes the previous decisions and history of uncertainty realizations as an input, and returns a new implementable decision. The dynamics of a sequential decision-making problem is depicted in Figure 1.

<figure><img src="image_3.png" /><figcaption>**Figure 1** Sequential decision-making under uncertainty</figcaption></figure>

There are several modeling frameworks for sequential decision-making problems under uncertainty. When the probability distribution governing the uncertain parameters is known, these problems may be addressed by the *multistage stochastic programming* (MSP) paradigm, with the goal of optimizing some statistical performance measure over the planning horizon. There is an extensive body of research on MSP problems of various structures, with a rich literature on problems with continuous decision variables. However, Bertsimas and Thiele (2006) point out that implemented solutions may perform poorly if the probability distribution used in the MSP model is different than the *true* distribution, even if both distributions share the same first and second moments. To mitigate this effect, *distributionally robust optimization* (DRO) models are proposed for making decisions that are based on a family of probability distributions, often defined by using historical data (Goh and Sim 2010). These models aim to hedge against tuning decisions to a perceived distribution. While the DRO framework has received significant attention (Mohajerin Esfahani and Kuhn (2018), Cheramin et al. (2022)) from the research community and some recent studies have proposed tractable solution methods for linear DRO problems under certain conditions (Philpott et al. 2018, Bertsimas et al. 2019), they remain largely challenging to solve especially in the multi-stage setting.

*Multistage adaptive/adjustable robust optimization* (MSARO) is another framework for modeling sequential decision-making problems under uncertainty that does not require any knowledge about the probability distribution governing the uncertain parameters. This framework is also adapted to contexts where the underlying uncertainty is not stochastic in nature, for instance, in the case of adversarial participants. In MSARO, the uncertainty is represented as belonging to a pre-structured (often compact) set, called the *uncertainty set*, and the decisions are optimized with respect to the worst-case outcome in this set.
and Nemirovski 1999); also see Ben-Tal et al. (2009), Bertsimas et al. (2011a, 2015)), MSARO in general remains computationally difficult to solve, even for two stages (Guslitser 2002, Ben-Tal et al. 2004, Feige et al. 2007). Consequently, most solution methods rely on an approximation of MSARO. In this work we present a new solution paradigm for MSARO problems by employing primal and dual decision rules that result in two-stage approximations of MSARO and its dual problem, respectively. This leads to a solution framework that returns an adaptable policy, and an optimality gap to assess the quality of the proposed policy. In what follows, we provide an overview of the existing literature on MSARO, followed by a summary of our contributions.

## 1.1. Literature Review

Static robust optimization is computationally less demanding than MSARO. As such, a number of studies attempt to find problem structures where the static optimal solution is also guaranteed to be an optimal solution of the MSARO (Ben-Tal et al. 2004, Marandi and Den Hertog 2018, Bertsimas et al. 2015). For example, this is the case for linear MSARO with constraint-wise uncertainty and a compact uncertainty set, provided that it is also linear with respect to uncertain parameters. In general, however, the solutions to the static and adaptive variants differ and the MSARO solutions are preferable by virtue of being more flexible. Figure 2 presents a summary of existing solution methods for obtaining exact, approximate, and dual bounds on the MSARO, with methods developed specifically for two-stage adaptive robust optimization (2ARO) separately categorized. In the following, we briefly discuss each method and the specific problem structure it can address.

Exact solution methods are scarce in the MSARO literature and the existing studies mostly focus on two-stage adaptive robust optimization (2ARO). Bertsimas et al. (2012) designed a Benders decomposition and outer approximation approach for a two-stage unit commitment problem. Zeng and Zhao (2013) developed constraint-and-column generation for 2ARO problems with fixed recourse and finite or polyhedral uncertainty set. For small linear problems, Zhen et al. (2018) used the Fourier-Motzkin elimination iteratively to remove the second-stage decisions in a 2ARO problem with continuous fixed recourse, eventually forming a static robust optimization problem that recovers optimal first-stage solutions. They extended their method to the multistage setting by modifying the algorithm to make sure that at every stage, no decision depends on the future information (these are the so-called nonanticipativity constraints). In the case of mixed-binary recourse and only objective uncertainty, Arslan and Detienne (2022) proposed an exact method based on a Dantzig-Wolfe reformulation of the recourse problem. In a rare instance of studies addressing multistage decision-making in the literature of robust optimization, Georghiou et al. (2019) developed robust dual dynamic programming (RDDP) for continuous MSARO problems with a stage-wise rectangular uncertainty set, and proved finite/asymptomatic convergence for
In this paper, we focus on the MSARO framework. Throughout, we use $[a] := \{1, 2, \dots, a\}$ and $[a, b] := \{a, a+1, \dots, b\}$ for positive integers $a$ and $b$ (with $a \le b$), and $(\cdot)^\top$ for the transpose operator. For a problem with $T$ decision-making stages, we denote with $\Xi^T \subseteq \mathbb{R}^{\ell^T}$ the uncertainty set which governs the set of uncertain parameters $(\xi_1, \xi_2, \dots, \xi_T)$, where $\xi_t$ denotes the vector of parameters associated with stage $t$, and $\xi_1 = 1$, by convention. For ease of presentation, we also define the sequence of uncertain parameter vectors up to stage $t$ along with their (projected) support as $\xi^t := (\xi_1, \dots, \xi_t) \in \Xi^t := \text{proj}_{\xi^t}(\Xi^T) \subseteq \mathbb{R}^{\ell^t}$. We assume $\Xi^t$ is compact for all $t \in [T]$. Then, we study the following general MSARO problem:

$$ \min_{x_1 \in X_1(\xi_1)} c_1(\xi_1)^\top x_1 + \sup_{\xi_2:(\xi^1, \xi_2) \in \Xi^2} \min_{\substack{x_2 \in X_2(\xi^2): \\ A_2(\xi^2)x_2 + B_2(\xi^2)x_1 \le b_2(\xi^2)}} c_2(\xi^2)^\top x_2 + \dots \quad (1) $$

$$ \dots + \sup_{\xi_t:(\xi^{t-1}, \xi_t) \in \Xi^t} \min_{\substack{x_t \in X_t(\xi^t): \\ A_t(\xi^t)x_t + B_t(\xi^t)x_{t-1} \le b_t(\xi^t)}} c_t(\xi^t)^\top x_t + \dots + \sup_{\xi_T:(\xi^{T-1}, \xi_T) \in \Xi^T} \min_{\substack{x_T \in X_T(\xi^T): \\ A_T(\xi^T)x_T + B_T(\xi^T)x_{T-1} \le b_T(\xi^T)}} c_T(\xi^T)^\top x_T $$

where $X_t(\xi^t) := \{x_t \in \mathbb{R}^{n_t-n_t^i} \times \mathbb{Z}^{n_t^i}: D_t(\xi^t)x_t \le d_t(\xi^t)\}$ and $c_t: \mathbb{R}^{\ell^t} \to \mathbb{R}^{n_t}$, $b_t: \mathbb{R}^{\ell^t} \to \mathbb{R}^{m_t^s}$, $d_t: \mathbb{R}^{\ell^t} \to \mathbb{R}^{m_t^R}$, $A_t: \mathbb{R}^{\ell^t} \to \mathbb{R}^{m_t^s \times n_t}$, $B_t: \mathbb{R}^{\ell^t} \to \mathbb{R}^{m_t^s \times n_{t-1}}$, $D_t: \mathbb{R}^{\ell^t} \to \mathbb{R}^{m_t^R \times n_t}$ for $t \in [T]$. The main output of model (1) is the first-stage deterministic (here-and-now) decisions $x_1$ which minimize the worst-case objective value over $T$ stages taking into account the sequential uncertainty realizations and optimal wait-and-see decisions. In the sequential framework, at each stage $t \in T$, the worst-case realization $\xi_t$ that is consistent with the history of the realizations up to stage $t-1$, $\xi^{t-1}$, is revealed. The vector $\xi_t$ combined with the history $\xi^{t-1}$ yields the history of realizations up to stage $t$, that is, $\xi^t = (\xi^{t-1}, \xi_t)$ from the support $\Xi^t$. This determines the parameters of the stage-t minimization problem, from which the optimal wait-and-see decision vector $x_t$ is obtained. As such, the wait-and-see decisions, $x_2, \dots, x_T$, also known as recourse decisions, are adapted to the history of the uncertain parameter realizations up to their decision-making stage, $\xi^2, \dots, \xi^T$, respectively. We remark that, by definitions of $X_t$ and $\Xi^t$, we allow for the possibility of mixed-integer wait-and-see decisions and dependence between the uncertain parameters of different stages. In the following, we assume, for the data, that all uncertain vectors and matrices are affine functions of the associated uncertain parameters $\xi^t$, as the majority of the literature mentioned makes this assumption; further assumptions will be specified explicitly when necessary. When the set $X_t(\xi^t), t \in [2, T]$ does not have integrality restrictions ($n_t^i = 0$), MSARO problem has continuous recourse, otherwise, it has (mixed-)integer recourse.

MSARO problems are highly challenging to solve. Indeed, as has been recently proven by Goerigk et al. (2024), they are, in general, harder than NP-hard problems, lying at the higher levels of the polynomial hierarchy, and their complexity increases with the number of decision stages. Specifically, $T$-stage MSARO problems with certain uncertainty set structures are $\Sigma_{2T-1}^P$-hard. However,
<figure><img src="image_3.png" /><figcaption>**Figure 2** Solution methods for the MSARO</figcaption></figure>

various problem sub-classes. RDDP is an adaptation of the stochastic dual dynamic programming algorithm from the stochastic programming literature (Pereira and Pinto 1991), to the robust setting. Georghiou et al. (2020) proposed a convergent method based on enumeration of the extreme points of the uncertainty set combined with affine decision rules to provide gradually improving primal and dual bounds for 2ARO problems with continuous recourse.

A commonly used approximation approach restricts adaptable/adjustable decisions to follow a certain functional form. These are known in the literature as *decision rules*. Ben-Tal et al. (2004) proposed the first decision rule for the MSARO, the affine/linear decision rule (LDR), where continuous recourse decisions are expressed as affine functions of uncertain parameters where the parameters of this function are to be optimized. This approach was further generalized by introducing polynomial decision rules to construct more flexible policies. For the MSARO with continuous fixed recourse, Bertsimas et al. (2011b) showed that polynomial policies with a predetermined degree can be computed using a semidefinite program (comparable to the work of Bampou and Kuhn (2011) for MSP). Xu and Hanasusanto (2021) proposed tighter approximations than that
some sub-classes and special cases of MSARO problems are theoretically and/or computationally
more tractable, depending on the number of stages, the structure of the uncertainty set, and the
nature of recourse decisions. Among these, the most well-known is *static* robust optimization, which
considers that all decisions are here-and-now. For a considerable number of problem structures,
e.g., when the uncertainty set is a polyhedron or an ellipsoid, it is possible to derive a mono-
lithic reformulation of the static robust optimization problem through a compact reformulation
of the adversarial problem, usually relying on duality techniques (Ben-Tal and Nemirovski 1999,
Ben-Tal et al. 2009, Bertsimas et al. 2011a, 2015). However, this paradigm cannot capture the
flexibility offered by the possibility of adapting some decisions to the realization of uncertainty,
thus often producing overly conservative decisions. As such, there has been significant research
effort on developing exact and approximate solution methods for two-stage adjustable robust opti-
mization (2ARO) problems, i.e., MSARO with $T = 2$ (Ben-Tal et al. 2004, Zeng and Zhao 2013,
Postek and Hertog 2016, Subramanyam et al. 2020). In developing these methods, the presence
of discrete recourse variables, known as (mixed-)integer recourse, poses additional challenges in
ensuring exact or high-quality solutions within a reasonable computational effort compared to the
continuous recourse case.

On the other hand, scientific progress on general MSARO has been much more limited. Given
the aforementioned theoretical complexity of these problems, the focus of most existing studies is
the approximate solution of these problems. Approximations proposed for MSARO mostly rely on
reducing the multi-stage problem to a static problem, with a view to leverage the tractability of
these problems. While these approximations can produce feasible solutions for MSARO problems,
they can lead to a significant degradation in solution quality. Furthermore, some of these methods
are quite restrictive being only applicable to special classes such as MSARO with continuous
recourse.

To address these limitations, this paper aims to develop approximations for general MSARO
problems of form (1). More specifically, inspired by the recent developments in the MSP literature,
namely the two-stage linear decision rules (Bodur and Luedtke 2018), we propose applying decision
rule approximations to only a certain subset of decision variables, resulting in an approximation of
MSARO problems as 2ARO problems. In so doing we have three motivations: (i) the strength of the
added flexibility in adapting recourse decisions to uncertainty in 2ARO compared to static robust
optimization, (ii) the significantly reduced theoretical complexity of 2ARO compared to MSARO,
and (iii) the progress made in computationally viable exact and approximate solution methodologies
for 2ARO. As a result of point (iii), our proposed framework is capable of considering a large
variety of MSARO classes, most notably the mixed-integer recourse case, and will directly benefit
of Ben-Tal et al. (2004) and Bampou and Kuhn (2011) (when using the same degree), by presenting exact copositive programming reformulations of the LDR-restricted MSARO with random recourse, and the quadratic decision rule-restricted problem with fixed recourse. Nonlinear decision rules were also proposed for the MSARO with continuous fixed recourse, such as deflected and segregated affine (Chen et al. 2008), extended affine (Chen and Zhang 2009), piecewise affine (Goh and Sim 2010), truncated linear (See and Sim 2010), and piecewise affine with exponentially many pieces (Ben-Tal et al. 2020). For a comprehensive list of nonlinear decision rules, interested reader may refer to a recent survey by Yanıkoğlu et al. (2019).

With mixed-integer recourse, LDRs and most of its extensions are not applicable. For the mixed-binary case, Bertsimas and Georghiou (2015) endogenously designed piecewise linear/constant decision rules. Later, Bertsimas and Georghiou (2018) introduced binary decision rules for the binary MSARO problems with fixed or random recourse. Using Dantzig-Wolfe reformulation, for a subclass of 2ARO problems with mixed-integer recourse and a finite uncertainty set, Hashemi Doulabi et al. (2021) proposed a static formulation which is amenable to Benders decomposition. A popular line of research for general mixed-integer problems uses the notion of *finite adaptability* first introduced by Bertsimas and Caramanis (2010) for linear MSARO. In finite or *K*-adaptability, the decision-maker a priori commits to *K* recourse decisions (while making the first-stage decisions), and then chooses among them after observing the uncertain parameters. Hanasusanto et al. (2015) tailored the method for problems with binary recourse, and Subramanyam et al. (2020) later generalized it for mixed-integer 2ARO. A related approach is based on the concept of uncertainty set partitioning which returns a sequence of improving approximations (Bertsimas and Dunning 2016, Postek and Hertog 2016, Romeijnders and Postek 2020). Although these latter methods can be used to approximate mixed-integer MSARO problems, they do not provide a fully adaptable policy.

The computational advantage obtained from the described approximation methods often comes at the expense of optimality. In a few special cases, LDRs are shown to be optimal for the MSARO. Bertsimas et al. (2010) showed that in a one-dimensional MSARO problem with no first-stage decisions, linear second-stage costs, and a hyperrectangle uncertainty set, affine policies are optimal. In a series of studies on 2ARO problems it was shown that the result still holds if we have only random right-hand side and/or technology matrix and a simplex uncertainty set (Bertsimas and Goyal 2012, Iancu et al. 2013, Zhen et al. 2018). Nevertheless, for the majority of the MSARO problems, decision rules yield suboptimal policies.

In order to assess the quality of the solutions of an LDR-restricted problem, Kuhn et al. (2011) applied LDRs on the dual variables associated with the linear programming dual of the original problem with no integer decisions, obtaining a dual bound and thus an optimality gap. From
from the developments in the highly active 2ARO literature in the future, e.g., the incorporation
of machine learning for computational enhancements (Julien et al. 2022, Dumouchelle et al. 2023).

While the aforementioned ideas are aimed at providing feasible solutions for MSARO problems,
an important question arises as to the quality of the obtained solution. In order to evaluate the
quality of a feasible policy, one could use a dual bound on the optimal value of the MSARO
problem. Unfortunately, obtaining dual bounds for MSARO problems is a largely unexplored topic
in the literature, especially in the case of mixed-integer recourse. To fill this gap, we propose to
develop dual approximations for MSARO. In particular, we propose a Lagrangian dual for the
MSARO problem and apply decision rules to the Lagrangian multipliers, leveraging ideas rooted
in the MSP literature, namely, Lagrangian dual decision rules (Daryalal et al. 2024). In deriving
a Lagrangian dual of the MSARO problem, we assign a probability distribution with the support
as the uncertainty set and use the assigned distribution in dualizing a subset of constraints along
with their Lagrangian multipliers. As a result, we obtain a dual approximation of MSARO in
the form of a two-stage stochastic optimization problem, which can be solved with the help of
state-of-the-art methods for two-stage stochastic problems. Since the quality of the resulting dual
bound depends on the assigned distribution used while developing the dual formulation, we define
a distribution optimization problem with the aim of identifying the strongest such dual bound.
We develop appropriate solution methods tailored to the nature of the recourse variables for the
resulting distribution optimization problem.

The contributions of our work are summarized as follows:

* We develop a solution framework for MSARO problems that returns adaptable policies as well as a dual bound measuring the quality of these policies. We do this by employing two-stage and Lagrangian dual decision rules, leading to novel techniques that reveal new theoretical and practical avenues.

* We propose to approximate MSARO problems via 2ARO problems to leverage existing solution methodologies and future developments for the latter in designing adaptable policies for the former. To this end, we present two-stage decision rules for MSARO, the first adaptation of a generalization of two-stage linear decision rules from the MSP literature in robust optimization, which can be applied to a broad range of problems. We employ, for an illustration of our approach, a tailored constraint-and-column generation algorithm to solve the resulting 2ARO approximation. The optimal solution of this approximation not only provides a primal policy but can also contribute to the calculation of a dual bound through identification of critical realizations in the uncertainty set.

* With a similar motivation, we derive a dual approximation of MSARO in the form of a two-stage stochastic optimization problem, which we show to be a strong dual in certain cases.
another perspective, Hadjiyiannis et al. (2011) solved an MSARO problem for a finite set of scenarios sampled from the uncertainty set to reach a dual bound. Inspired by these studies, Bertsimas and de Ruiter (2016) proposed a reformulation of MSARO models with polyhedral uncertainty, by dualizing over the recourse decisions and then uncertain parameters, leading to an equivalent problem that achieves a dual bound (potentially different than the one obtained from the original nested formulation) for a finite set of scenarios, that is the binding scenarios in an LDR-restricted approximation of the problem.

## 1.2. Contributions

In this work, we present decision rule-based approaches to obtain primal and dual bounds for an MSARO problem. On the primal side, we apply two-stage decision rules that approximate the multistage problem with a 2ARO, and leverage existing solution methods from the 2ARO literature to solve the obtained approximate model. This results in feasible adjustable policies. On the dual side, we use Lagrangian duality and *Lagrangian dual decision rules* for deriving a dual problem for an MSARO problem, for which we develop reformulations and solution methods. Accordingly, we obtain primal and dual bounds and an optimality gap for measuring the quality of the proposed policies. Figure 3 provides an overview of our solution framework.

<figure><img src="image_4.png" /><figcaption>**Figure 3** Proposed bounding framework for MSARO</figcaption></figure>

More specifically, our new policies are inspired by two numerically promising decision rules recently developed in the stochastic programming community for MSP: two-stage linear decision rules (Bodur and Luedtke 2018) and Lagrangian dual decision rules (Daryalal et al. 2021), that approximate the MSP with a two-stage problem in the primal and dual space, respectively. On
Moreover, in order to obtain the strongest possible such dual bound, we study the numerical
solution of a distribution optimization problem. More specifically, we apply decision rules to
dual variables, and design a cutting-plane algorithm to solve the obtained restricted dual prob-
lem. We also show that in the special case of continuous recourse, the restricted dual problem
can be reformulated as a monolithic bilinear program. Additionally, we present an alternative
decomposable dual problem which offers an improved numerical performance. These novel
techniques contribute to the scarce literature for obtaining dual bounds for MSARO problems
with mixed-integer recourse.

* We evaluate the performance of our solution framework over multistage versions of three classical problems from the MSARO literature: (i) the newsvendor problem, (ii) the location-transportation problem, and (iii) the capital budgeting problem. Each of these problem classes is suitable for a different solution method developed in this work and our analysis over various instances attests to the quality of the returned primal and dual bounds.

The remainder of the paper is organized as follows. In Section 2 we review the literature relevant to our work. In Section 3 we introduce two-stage decision rules for obtaining primal adaptable policies. In Section 4 we present our approach to deriving dual bounds. This is followed by numerical experiments in Section 5 and concluding remarks. We remark that all proofs are deferred to Appendix B.

## 2. Literature Review

Figure 2 presents a summary of existing solution methods for obtaining exact/approximate solutions and dual bounds for MSARO, with methods developed specifically for 2ARO separately categorized. In the following, we briefly discuss each method and the specific problem structure it can address.

Exact solution methods are scarce in the MSARO literature and the existing studies mostly focus on 2ARO. For 2ARO problems with fixed recourse and finite or polyhedral uncertainty set, Zeng and Zhao (2013) developed a constraint-and-column generation algorithm. For the same type of problems restricted to continuous recourse, Bertsimas et al. (2012) designed a Benders decomposition-type algorithm and applied it to a unit commitment problem, whereas Georghiou et al. (2020) proposed a convergent method based on enumeration of the extreme points of the uncertainty set combined with affine decision rules to provide gradually improving primal and dual bounds. For 2ARO problems with continuous fixed recourse, Zhen et al. (2018) used Fourier-Motzkin elimination iteratively to remove the second-stage decisions, eventually forming an equivalent static robust optimization problem. This computationally expensive approach is also extended to multi-stage problems. In the case of mixed-binary recourse and only objective uncertainty, Arslan and
the primal side we have generalized the idea of two-stage LDRs and designed a new decision
rule with a broader application, including some classes of MSARO with mixed-integer recourse.
Bringing this idea to the context of robust optimization makes it possible to benefit from the
rich literature of 2ARO in (approximately) solving multistage problems. On the dual side, the
literature of stochastic programming and distributionally robust optimization (DRO) prove to be
useful in the development of a new dual problem for MSARO, building a bridge between these
close communities. As in our work in Daryalal et al. (2021), dual solutions can also be potentially
used to derive primal policies, which we leave for future research. A summary of our contributions
is given below:

* For a broad class of multistage adaptive robust optimization problems, we propose a solution framework that returns adaptable policies as well as a dual bound measuring the quality of these policies. We do this by employing two-stage and Lagrangian dual decision rules, leading to novel techniques that reveal new theoretical and practical avenues.

* We present the class of two-stage decision rules that impose a predefined structure on a subset of the adaptable/adjustable decision variables, with the goal of building an approximate model that returns feasible first-stage solutions. The result is a 2ARO approximation of the original multistage problem. We show that in some special cases, constraint-and-column generation algorithm converges to the exact solution of the 2ARO approximation, otherwise the K-adaptability approach is applicable.

* For a general MSARO, we present a dual problem and prove that, under certain conditions, strong duality holds. This is one of the few approaches for obtaining dual bounds for MSARO with integer recourse. To the best of our knowledge, the only other approach that applies to mixed-integer recourse is discretization/sampling. For this dual, we introduce Lagrangian dual decision rules, resulting in a two-stage problem. We design a cutting-plane algorithm that solves the restricted dual problem for MSARO with mixed-integer recourse. We also show that in the case of continuous recourse, the problem can be reformulated as a monolithic bilinear program. Additionally, we present an alternative dual problem that is decomposable by scenario and offers computational efficacy.

* We have evaluated the performance of our solution framework over multistage versions of three classical problems in the literature of MSARO: (i) the newsvendor problem, (ii) the location-transportation problem, and (iii) the capital budgeting problem. Each of these problem classes are suitable for a different solution method developed in this work and our analysis over various instances attests to the quality of the returned primal and dual bounds.

The remainder of this article is organized as follows. In Section 2 we introduce the notion of a state
variable and design two-stage decision rules for retrieving feasible solutions and adaptable policies.
<figure><img src="image_4.png" /><figcaption>**Figure 2** Solution methods for MSARO</figcaption></figure>

Detienne (2022) proposed an exact method based on a Dantzig-Wolfe reformulation of the recourse problem based on a technical assumption on the structure of the linking constraints. Similarly, using Dantzig-Wolfe reformulation, for a subclass of 2ARO problems with fixed and mixed-integer recourse, block diagonal recourse matrix and a finite uncertainty set, Hashemi Doulabi et al. (2021) derived a static formulation which is amenable to Benders decomposition. For continuous MSARO problems with a stage-wise rectangular uncertainty set, Georghiou et al. (2019) developed robust dual dynamic programming (RDDP) and proved finite/asymptomatic convergence for various problem sub-classes. RDDP is an adaptation of the stochastic dual dynamic programming algorithm from the MSP literature (Pereira and Pinto 1991) to MSARO.

Approximate solution methods are more common in the MSARO literature, with the central idea of restricting adaptable/adjustable decisions to follow a certain functional form, known as *decision rules*. Ben-Tal et al. (2004) proposed the first decision rule for MSARO problems with *continuous recourse*, LDRs, where recourse decisions are expressed as affine functions of uncertain parameters where the parameters of this function are to be optimized. The resulting LDR-restricted problem being a static optimization problem, it can be reformulated as a linear optimization problem in certain cases. Nonlinear decision rules were also proposed, such as deflected and segregated
In Section 3 we reformulate the problem by employing the nonanticipativity constraints, and then propose a dual problem through their Lagrangian relaxation. To solve this dual problem, we apply linear decision rules on the Lagrangian multipliers of the relaxed nonanticipativity constraints. We formalize the problem of obtaining the best such dual bound, and develop solution methods depending on the nature of the recourse problem. This is followed by numerical experiments in Section 4 and concluding remarks.

*Notation.* Throughout, we use $[a] := \{1, 2, \dots, a\}$ and $[a, b] := \{a, a+1, \dots, b\}$ for positive integers $a$ and $b$ (with $a \le b$), and $(\cdot)^{\top}$ for the transpose operator.

## 2. Primal Bounding

In this work, we study MSARO problems of the following form:

$$ \nu^* = \min z \qquad (1a) $$

$$ \text{s.t.} \quad \sum_{t \in [T]} c_t (\xi^t)^{\top} x_t (\xi^t) \leq z \qquad \xi^T \in \Xi \qquad (1b) $$

$$ A_t(\xi^t)x_t(\xi^t) + B_t(\xi^t)x_{t-1}(\xi^{t-1}) \leq b_t(\xi^t) \qquad t \in [2, T], \xi^t \in \Xi^t \qquad (1c) $$

$$ D_t(\xi^t)x_t(\xi^t) \leq d_t(\xi^t) \qquad t \in [T], \xi^t \in \Xi^t \qquad (1d) $$

$$ x_t(\xi^t) \in \mathbb{R}^{n_t - n_t^i} \times \mathbb{Z}^{n_t^i} \qquad t \in [T], \xi^t \in \Xi^t, \qquad (1e) $$

where $T$ denotes the number of decision stages, $c_t: \mathbb{R}^{\ell^t} \to \mathbb{R}^{n_t}$, $b_t: \mathbb{R}^{\ell^t} \to \mathbb{R}^{m_t^s}$, $d_t: \mathbb{R}^{\ell^t} \to \mathbb{R}^{m_t^R}$, $A_t: \mathbb{R}^{\ell^t} \to \mathbb{R}^{m_t^s \times n_t}$, $B_t: \mathbb{R}^{\ell^t} \to \mathbb{R}^{m_t^s \times n_{t-1}}$, $D_t: \mathbb{R}^{\ell^t} \to \mathbb{R}^{m_t^R \times n_t}$. The uncertainty is represented by $\{\xi_t\}_{t=1}^T$ with support $\Xi$ where $\xi_t$ is an uncertain vector of random variables at stage $t$ with outcomes $\xi_t \in \Xi_t = \text{proj}_{\xi_t}(\Xi) \subseteq \mathbb{R}^{\ell^t}$ (where by convention $\xi_1 = 1$). A realization of uncertainty is referred to as a *scenario*, $\{\xi_t\}_{t=1}^T$, which is a set of outcomes for all random variables at every stage such that $\xi_t \in \Xi_t$. History of uncertainty at stage $t$ is characterized by $\xi^t = (\xi_1, \dots, \xi_t)$. Decision variables are denoted by $x_t(\xi^t) \in X_t(\xi^t) = \{x_t \in \mathbb{R}^{n_t - n_t^i} \times \mathbb{Z}^{n_t^i} : D_t(\xi^t)x_t \leq d_t(\xi^t)\}$, i.e., nonanticipative policies that only depend on the history of uncertainty. Throughout, we make the following assumptions:

**ASSUMPTION 1.** The problem has relatively complete recourse, meaning that for all $t \in [1, T]$, $\xi^t$ and a history of feasible decisions $x_t(\xi^t)$ made up to $t$, there always exists a feasible decision at stage $t$.

**ASSUMPTION 2.** For $t \in [T]$ and $\xi^t \in \Xi$, the maximum diameter of $X_t(\xi^t)$ is finite, i.e., the feasibility sets are bounded.

When the set $X_t(\xi^t)$, $t \in [T]$ does not have integrality restrictions ($n_t^i = 0$), we have continuous recourse, otherwise we are dealing with an MSARO problem with mixed-integer recourse. Together with constraints (1b), the objective function (1a) minimizes the worst outcome. Constraints (1c)
affine (Chen et al. 2008), extended affine (Chen and Zhang 2009), piecewise affine (Goh and Sim 2010), truncated linear (See and Sim 2010), piecewise affine with exponentially many pieces (Ben-Tal et al. 2020), quadratic (Xu and Hanasusanto 2021) and polynomial (Bertsimas et al. 2011b) decision rules. However, the resulting reformulations when using non-linear decision rules are often nonlinear, e.g. semidefinite or copositive programs. For a comprehensive list of nonlinear decision rules, interested reader may refer to the survey by Yanıkoğlu et al. (2019).

In the case of mixed-integer recourse, LDRs and most of its aforementioned extensions lead to non-adjustable decisions for the integer variables. Thus, alternative approaches have been proposed, with the key idea of (implicitly or explicitly) partitioning the uncertainty set and determining a constant recourse solution corresponding to each subset. A popular approach for 2ARO problems uses the notion of finite adaptability, first introduced by Bertsimas and Caramanis (2010). In finite or *K*-adaptability, the decision-maker a priori commits to *K* recourse decisions (while making the first-stage decisions), and then chooses among them after observing the uncertainty realization which leads to an implicit *K*-partition of the uncertainty set. While Bertsimas and Caramanis (2010) presented an exact formulation for the 2-adaptability problem, for the general *K*-adaptability case, Hanasusanto et al. (2015) proposed a monolithic formulation for problems with binary recourse, and Subramanyam et al. (2020) developed a branch-and-bound algorithm for problems with mixed-integer recourse. For MSARO problems on the other hand, explicit uncertainty set partitioning is considered in an iterative heuristic framework, with the aim of obtaining a sequence of improving approximations (Bertsimas and Dunning 2016, Postek and Hertog 2016, Romeijnders and Postek 2020). Lastly, for MSARO problems with pure-binary recourse, Bertsimas and Georghiou (2018) introduced binary decision rules, whereas for the mixed-binary recourse case, Bertsimas and Georghiou (2015) implicitly designed piecewise linear/constant decision rules.

While the aforementioned primal approximations can be shown to be exact in some special cases (Bertsimas et al. 2010, Bertsimas and Goyal 2012, Iancu et al. 2013, Hanasusanto et al. 2015, Zhen et al. 2018), in general they do not provide optimal solutions. In order to assess the quality of their feasible solutions, dual bounds can be used. To this end, Kuhn et al. (2011) presented the idea of deriving a dual problem for MSARO with only continuous variables and applying LDRs on the dual variables. Since their approach was originally derived for stochastic programs its application to MSAROs requires assigning a probability distribution to the uncertainty set. The impact of the chosen distribution on the quality of the obtained dual bound was observed by Kuhn et al. (2011), as such a distribution optimization problem was mentioned. This problem was later formalized by Hadjiyiannis et al. (2011) for a 2ARO with continuous variables and shown to be of the same theoretical difficulty as the original problem. Hadjiyiannis et al. (2011) proposed to solve instead a 2ARO problem for a finite set of scenarios from the uncertainty set, selected based on a primal
and (1d) are state and recourse constraints, respectively: while the former link different stages, the latter are local restrictions for a specific stage.

For an MSARO with continuous recourse, it is quite common to restrict the decisions $x_t(\xi^t)$ to an affine form, which, by breaking the temporal dependencies between stages, approximates problem (1) with a static robust optimization problem. Let $\nu^{LDR}$ be the bound obtained by applying such linear decision rules (LDRs) to the set of all decision variables. Our goal in this section is to employ a new paradigm where a specific subset of the decision variables $x_t(\xi^t)$ are enforced to follow a structured decision rule, leading to an approximation in the form of a two-stage adaptive robust optimization problem.

## 2.1. Two-stage Decision Rules

In a similar manner to constraints (1c) and (1d), let us partition the decision variables $x_t(\xi^t)$, $t \in [2, T]$ into $x_t^s(\xi^t) \in \mathbb{R}^{qt}$ and $x_t^r(\xi^t) \in \mathbb{R}^{pt}$, state and recourse variables, respectively, with $q_t + p_t = n_t$, with integrality restrictions, if any, embedded in the set $X_t(\xi^t)$. Then the MSARO can be reformulated as follows:

$$ \nu^* = \min z \qquad (2a) $$

$$ \text{s.t.} \quad \sum_{t \in [T]} c_t^s(\xi^t)^\top x_t^s(\xi^t) + c_t^r(\xi^t)^\top x_t^r(\xi^t) \leq z \qquad \xi^T \in \Xi \qquad (2b) $$

$$ A_t^s(\xi^t)x_t^s(\xi^t) + B_t^s(\xi^t)x_{t-1}^s(\xi^{t-1}) + A_t^r(\xi^t)x_t^r(\xi^t) \leq b_t(\xi^t) \quad t \in [2, T], \xi^t \in \Xi^t \qquad (2c) $$

$$ (x_t^s(\xi^t), x_t^r(\xi^t)) \in X_t(\xi^t) \qquad t \in [2, T], \xi^t \in \Xi^t \qquad (2d) $$

$$ x_1 \in X_1, \qquad (2e) $$

where $x_1^s(\xi^1) = x_1$, and $c_t^s(\xi^t)$, $A_t^s(\xi^t)$, $B_t^s(\xi^t)$ are sub-arrays/sub-matrices of $c_t(\xi^t)$, $A_t(\xi^t)$, $B_t(\xi^t)$ associated with state variables with appropriate dimensions, while $c_t^r(\xi^t)$, $A_t^r(\xi^t)$ have the same role for the recourse variables. Let $x_t^s(\xi^t)$ be approximated by a decision rule, i.e., $x_t^s(\xi^t) = \Theta_t(\xi^t, \beta_t)$, where $\Theta_t: \mathbb{R}^{\ell^t} \times \mathbb{R}^{K_t} \rightarrow \mathbb{R}^{qt}$ is a given decision rule, and $\beta_t \in \mathbb{R}^{K_t}$ is its vector of design parameters. By substituting the decision rule in problem (2), we obtain an approximation that can be reformulated as:

$$ \nu^{2s} = \min c_1^\top x_1 + \max_{\xi^T \in \Xi} \min_{x^r \in \mathcal{X}(x_1, \beta, \xi^T)} \sum_{t \in [2, T]} c_t^s(\xi^t)^\top \Theta_t(\xi^t, \beta_t) + c_t^r(\xi^t)^\top x_t^r \qquad (3a) $$

$$ \text{s.t.} \quad x_1 \in X_1 \qquad (3b) $$

$$ \beta_t \in \mathbb{R}^{K_t} \quad t \in [2, T], \qquad (3c) $$

where:

$$ \mathcal{X}(x_1, \beta, \xi^T) = \{ (x_t^r(\xi^t))_{t \in [2, T]} \in \mathbb{R}^{p_2} \times \mathbb{R}^{p_3} \times \cdots \times \mathbb{R}^{p_T} : $$
decision rule restriction, to reach a dual bound. This procedure can also be extended to obtain dual bounds for general MSAROs. Finally, for special cases of 2ARO problems with continuous recourse, Georghiou et al. (2020) proposed a framework to derive progressive dual bounds, by considering the linear programming dual of their primal extreme point reformulation.

### 3. Primal Bounding

The MSARO problem given in a nested form in (1) can be reformulated as a monolithic optimization problem by explicitly introducing the functional form of the decision variables, $x_t(\xi^t)$: $\Xi^t \to \mathbb{R}^{n_t - n_t^1} \times \mathbb{Z}^{n_t^1}$ for all $t \in [T]$, along with a deterministic variable, $z$, representing the worst-case objective value:

$$ \nu^* := \min z \qquad (2a) $$

$$ \text{s.t.} \quad \sum_{t \in [T]} c_t(\xi^t)^\top x_t(\xi^t) \le z \qquad \xi^T \in \Xi^T \qquad (2b) $$

$$ A_t(\xi^t)x_t(\xi^t) + B_t(\xi^t)x_{t-1}(\xi^{t-1}) \le b_t(\xi^t) \qquad t \in [2, T], \xi^t \in \Xi^t \qquad (2c) $$

$$ x_t(\xi^t) \in X_t(\xi^t) \qquad t \in [T], \xi^t \in \Xi^t. \qquad (2d) $$

Together with constraints (2b), the objective function (2a) minimizes the worst outcome. Constraints (2c) and (2d) are state and recourse constraints, respectively: while the former link different stages, the latter are local restrictions for a specific stage.

Throughout the paper, we make the following assumptions:

**ASSUMPTION 1.** $X_1(\xi^1)$ is non-empty.

**ASSUMPTION 2.** The problem has relatively complete recourse, i.e., for all $t \in [2, T]$, $\xi^t \in \Xi^t$ and a history of feasible decisions made up to $t$, $\{x_{t'}(\xi^{t'})\}_{t' \in [t-1]}$, there always exists a feasible decision at stage $t$, $x_t(\xi^t)$.

**ASSUMPTION 3.** For $t \in [T]$ and $\xi^t \in \Xi^t$, the maximum diameter of $X_t(\xi^t)$ is finite, i.e., the feasibility sets are bounded.

We note that combined with the compactness assumption of the uncertainty set and its projections, these assumptions imply that the studied MSARO problem has a finite optimal objective value.

To derive feasible policies to the MSARO problem, it is quite common in the literature to restrict all the decisions $x_t(\xi^t)$ to follow a simple functional form, such as an affine or piecewise constant decision rule. By breaking the temporal dependencies between stages, this approach approximates problem (2) with a static robust optimization problem. Our goal in this section is to employ a new paradigm where a specific subset of the decision variables $x_t(\xi^t)$ are enforced to follow a structured decision rule, leading to a restriction in the form of a 2ARO problem. We introduce this approach, two-stage decision rules for MSARO, in Section 3.1, then present its specific instantiation and possible solution methodologies in Sections 3.2 and 3.3.
$$ 
\left. 
\begin{array}{l}
A_t^r(\xi^t)x_t^r \le b_t(\xi^t) - (A_t^s(\xi^t)\Theta_t(\xi^t, \beta_t) + B_t^s(\xi^t)\Theta_{t-1}(\xi^{t-1}, \beta_{t-1})) \quad t \in [2, T] \\
(\Theta_t(\xi^t, \beta_t), x_t^r) \in X_t(\xi^t) \qquad t \in [2, T]
\end{array}
\right\}.
$$

Note that the decision rules are solely applied to the state variables, whereas recourse variables remain fully adjustable to the uncertain parameters (see Figure 4). Problem (3) is a 2ARO and the temporal dependency between stages is removed thanks to the application of the two-stage decision rule. A practical result of such an approximation is that the resulting problem can be solved using the existing solution methods for the 2ARO. In the following sections, we study two special cases that directly benefit from this literature. The first one addresses MSARO problems with continuous state variables, and the second can be applied to those problems that have only bounded integer state variables. Together, they permit the approximation of a (bounded) mixed-integer MSARO with a 2ARO.

<figure><img src="image_5.png" /><figcaption>Figure 4 Two-stage decision rules</figcaption></figure>

## 2.2. Two-stage Linear Decision Rules

If the state variables are continuous, we can approximate them with an affine function. Let $\Phi_t(\xi^t) = (\Phi_{t1}(\xi^t), \dots, \Phi_{tK_t}(\xi^t)) : \mathbb{R}^{\ell^t} \to \mathbb{R}^{K_t}$ be a vector of basis functions. The two-stage LDR is enforced by using $\Theta_t(\xi^t, \beta_t) = \Phi_t(\xi^t)^\top \beta_t$ in (3). The result is a two-stage mixed-integer linear program:

$$ \nu^{2S-LDR} = \min \ c_1^\top x_1 + \max_{\xi^T \in \Xi} \ \min_{x^r \in \mathcal{X}(x_1, \beta, \xi^T)} \ \sum_{t \in [2, T]} c_t^s(\xi^t)^\top \Phi_t(\xi^t)^\top \beta_t + c_t^r(\xi^t)^\top x_t^r \quad (4a) $$

$$ \text{s.t.} \quad x_1 \in X_1 \qquad (4b) $$

$$ \beta_t \in \mathbb{R}^{K_t} \quad t \in [2, T]. \qquad (4c) $$
## 3.1. Two-stage Decision Rules

In a similar manner to constraints (2c) and (2d), we partition the decision variables $x_t(\xi^t)$, $t \in [T]$ into $x_t^s(\xi^t) \in \mathbb{R}^{qt}$ and $x_t^r(\xi^t) \in \mathbb{R}^{pt}$, state and recourse variables, as those that appear in the state constraints of subsequent stages and the others, respectively. We have $q_t + p_t = n_t$, with integrality restrictions on the variables, if any, embedded in the set $X_t(\xi^t)$. Then the MSARO (2) can be written more explicitly as follows:

$$ \nu^* = \min z \tag{3a} $$

$$ \text{s.t.} \quad \sum_{t \in [T]} c_t^s(\xi^t)^\top x_t^s(\xi^t) + c_t^r(\xi^t)^\top x_t^r(\xi^t) \le z \qquad \xi^T \in \Xi^T \tag{3b} $$

$$ A_t^s(\xi^t)x_t^s(\xi^t) + B_t^s(\xi^t)x_{t-1}^s(\xi^{t-1}) + A_t^r(\xi^t)x_t^r(\xi^t) \le b_t(\xi^t) \quad t \in [2, T], \xi^t \in \Xi^t \tag{3c} $$

$$ (x_t^s(\xi^t), x_t^r(\xi^t)) \in X_t(\xi^t) \qquad t \in [T], \xi^t \in \Xi^t \tag{3d} $$

where $c_t^s(\xi^t)$, $A_t^s(\xi^t)$, $B_t^s(\xi^t)$, $D_t^s(\xi^t)$ are sub-arrays/sub-matrices of $c_t(\xi^t)$, $A_t(\xi^t)$, $B_t(\xi^t)$, $D_t(\xi^t)$ associated with the state variables with appropriate dimensions, while $c_t^r(\xi^t)$, $A_t^r(\xi^t)$, $D_t^r(\xi^t)$ have the same role for the recourse variables, and $X_t(\xi^t) = \{x_t^s \in \mathbb{R}^{q_t - q_t^i} \times \mathbb{Z}^{q_t^i}, x_t^r \in \mathbb{R}^{p_t - p_t^i} \times \mathbb{Z}^{p_t^i} : D_t^s(\xi^t)x_t^s + D_t^r(\xi^t)x_t^r \le d_t(\xi^t)\}$. For notational convenience, we drop the parametrization for the first-stage variables as well as their feasible set and the objective vector, i.e., use $x_1 = (x_1^s, x_1^r)$, $X_1$, and $c_1 = (c_1^s, c_1^r)$.

For $t \in [2, T]$, let $x_t^s(\xi^t)$ be approximated by a decision rule, i.e., $x_t^s(\xi^t) = \Theta_t(\xi^t, \beta_t)$, where $\Theta_t: \mathbb{R}^{l^t} \times \mathbb{R}^{K_t} \to \mathbb{R}^{qt}$ represents the rule, and $\beta_t \in \mathbb{R}^{K_t}$ represents its vector of design parameters. By substituting this rule in problem (3), we obtain an approximation that can be reformulated as:

$$ \nu^{2S} := \min c_1^\top x_1 + \max_{\xi^T \in \Xi^T} \min_{x^r \in \mathcal{X}(x_1^s, \beta, \xi^T)} \sum_{t \in [2, T]} c_t^s(\xi^t)^\top \Theta_t(\xi^t, \beta_t) + c_t^r(\xi^t)^\top x_t^r \tag{4a} $$

$$ \text{s.t.} \quad x_1 \in X_1 \tag{4b} $$

$$ \beta_t \in \mathbb{R}^{K_t} \quad t \in [2, T], \tag{4c} $$

where:

$$ \begin{aligned} \mathcal{X}(x_1^s, \beta, \xi^T) := & \{ (x_t^r)_{t \in [2, T]} \in \mathbb{R}^{p_2} \times \mathbb{R}^{p_3} \times \cdots \times \mathbb{R}^{p_T} : \\ & A_t^r(\xi^t)x_t^r \le b_t(\xi^t) - (A_t^s(\xi^t)\Theta_t(\xi^t, \beta_t) + B_t^s(\xi^t)x_1^s) \quad t=2 \\ & A_t^r(\xi^t)x_t^r \le b_t(\xi^t) - (A_t^s(\xi^t)\Theta_t(\xi^t, \beta_t) + B_t^s(\xi^t)\Theta_{t-1}(\xi^{t-1}, \beta_{t-1})) \quad t \in [3, T] \\ & (\Theta_t(\xi^t, \beta_t), x_t^r) \in X_t(\xi^t) \quad t \in [2, T] \}. \end{aligned} $$

Note that the decision rules are solely applied to the state variables, whereas the recourse variables remain fully adjustable to the uncertain parameters (see Figure 3). Problem (4) is a 2ARO since the temporal dependency between stages is removed thanks to the application of the two-stage
Since recourse variables $x_t^r(\xi^t)$ in problem (4) are fully adjustable, as opposed to the bound $\nu^{LDR}$ obtained by restricting them to LDRs, it immediately follows for MSARO problems with continuous state variables, that

$$ \nu^* \le \nu^{2S-LDR} \le \nu^{LDR}. $$

In general, we can approximate the 2ARO problem (4) with the K-adaptability approach (see Section 1.1 for an overview of the related literature), while under certain conditions it can be solved exactly. The focus of the following is on the latter case, where we assume that the original MSARO problem has fixed recourse, i.e., $c_t(\xi^t) = c_t, A_t(\xi^t) = A_t, B_t(\xi^t) = B_t, D_t(\xi^t) = D_t, t \in [2, T]$.

For solving a 2ARO problem with $\Xi$ either a finite discrete set or a polyhedron, and only right-hand side uncertainty, Zeng and Zhao (2013) proposed the constraint-and-column generation (C&CG) method, by drawing on the fact that not all realizations in $\Xi$ contribute to the optimal objective value. The C&CG strives to identify the most significant scenarios by starting from a smaller uncertainty set and gradually expanding it, which leads to the generation of new columns and constraints, respectively corresponding to recourse variables and second-stage constraints for the newly identified scenario. The C&CG algorithm of Zeng and Zhao (2013) assumes fixed technology matrix and first-stage costs, i.e., that the coefficients of the first-stage variables in the objective function and constraints are deterministic. However, in problem (4), basis functions $\Phi_t(\xi^t)$ appear as coefficients of $\beta_t$, both in the objective function (4a) (in the term $c_t^s(\xi^t)^\top \Phi_t(\xi^t)$) and the second-stage constraints (terms $A_t^s \Phi_t(\xi^t)$, $B_t^s \Phi_{t-1}(\xi^{t-1})$ and $D_t^s \Phi_t(\xi^t)$). Next, we present the adaptation of the C&CG algorithm for solving problem (4) when the technology matrix and first-stage costs of the MSARO are not fixed, followed by a sufficient condition under which the C&CG algorithm is finitely convergent.

**Constraint-and-column Generation for Two-stage LDRs**

Consider a relaxation of problem (4) where instead of considering the uncertainty set $\Xi$, a potentially empty subset of it, $\hat{\Xi} \subseteq \Xi$, is used. With an attempt to build a stronger relaxation, the C&CG gradually expands $\hat{\Xi}$ with what we call significant scenarios. The optimal objective value of the following problem is a lower bound on $\nu^{2S-LDR}$:

$$ \min_{c_1^\top x_1 + \eta} \qquad (5a) $$

$$ \text{s.t.} \quad \eta \ge \sum_{t \in [2, T]} (c_t^{s\top} \Phi_t(\xi^t)^\top \beta_t + c_t^{r\top} x_{t, \xi^t}^r) \qquad \xi^T \in \Xi^{\text{FEAS}} \qquad (5b) $$

$$ A_t^r x_{t, \xi^t}^r + A_t^s \Phi_t(\xi^t)^\top \beta_t + B_t^s \Phi_{t-1}(\xi^{t-1})^\top \beta_{t-1} \le b_t(\xi^t) \qquad t \in [2, T], \xi^t \in \hat{\Xi}^t \qquad (5c) $$

$$ (\Phi_t(\xi^t)^\top \beta_t, x_{t, \xi^t}^r) \in X_t(\xi^t) \qquad t \in [2, T], \xi^t \in \hat{\Xi}^t \qquad (5d) $$
<figure><img src="image_5.png" /><figcaption>**Figure 3** Two-stage decision rules</figcaption></figure>

decision rules (see Appendix A for a detailed proof). We remark that the relatively complete recourse assumption stated in Assumption 2 for the MSARO problem, does not guarantee that the 2ARO model has relatively complete recourse, but it can be ensured, for instance, by following the techniques mentioned in (Bodur and Luedtke 2018).

A practical result of such an approximation is that the resulting problem can be solved using the existing solution methods for 2ARO. In the following sections, we present two possible choices for the decision rule $\Theta_t(\xi^t, \beta_t)$, respectively applicable to continuous and integer state variables. Together, they permit the approximation of an MSARO problem with mixed-integer state variables via a 2ARO model. We remark that the nature of the recourse variables does not impact the reduction of an MSARO to a 2ARO, but plays an important role in the choice of an appropriate solution method for the resulting 2ARO model. In what follows, we illustrate the application of these decision rules and the algorithmic solution of ensuing 2ARO models. For ease of exposition, we present MSARO problems with only continuous and only integer state variables separately.

## 3.2. Two-stage Linear Decision Rules for MSAROs with Continuous State Variables

If the state variables are continuous, we can approximate them via a decision rule with an affine form. For $t \in [2, T]$, letting $\Phi_t(\xi^t) = (\Phi_{t1}(\xi^t), \dots, \Phi_{tK_t}(\xi^t)) : \mathbb{R}^{\ell^t} \to \mathbb{R}^{q_t \times K_t}$ be a vector of chosen basis functions, the two-stage LDR is enforced by using

$$ \Theta_t(\xi^t, \beta_t) = \Phi_t(\xi^t)\beta_t \quad (5) $$

in (4) where we use a compact matrix representation<fnref n="1" /> for notational convenience. The resulting 2ARO is written as:

$$ \nu^{2\text{S-LDR}} := \min c_1^\top x_1 + \max_{\xi^T \in \Xi^T} \min_{x^r \in \mathcal{X}(x_1^s, \beta, \xi^T)} \sum_{t \in [2,T]} c_t^s(\xi^t)^\top \Phi_t(\xi^t)\beta_t + c_t^r(\xi^t)^\top x_t^r \quad (6a) $$

<footnote n="1">¹ This representation is obtained, without loss of generality, by concatenating individual LDR restrictions applied to each state variable $x_{ti}^s(\xi^t)$ for $i \in [q_t]$. For example, consider an instance where there are two state variables at stage $t=2$ and the history consists of two components $\xi^2 = (\xi_1, \xi_2) \in \mathbb{R}^2$ (with the convention that $\xi_1 = 1$). Then the decisions rules $x_{21}(\xi^2) = \beta_{21}^1\xi_1 + \beta_{22}^1\xi_2$ and $x_{22}(\xi^2) = \beta_{21}^2\xi_1 + \beta_{22}^2\xi_2$ can be represented in the more compact matrix form with $K_2=4$ using the concatenated decision vector $\beta \in \mathbb{R}^4$ and the basis function matrix $\Phi_2(\xi^2) = \begin{bmatrix} \xi_1 & \xi_2 & 0 & 0 \\ 0 & 0 & \xi_1 & \xi_2 \end{bmatrix}$.</footnote>
$$x_1 \in X_1, \quad \eta \in \mathbb{R} \qquad (5e)$$

$$\beta_t \in \mathbb{R}^{K_t} \qquad t \in [2, T], \qquad (5f)$$

where $\hat{\Xi}^t = \text{proj}_{\xi^t}(\hat{\Xi})$, and $\Xi^{\text{FEAS}} \subseteq \hat{\Xi}$ is the set of identified significant scenarios $\xi^T$ for which there exists a feasible first-stage solution $\hat{\beta}$ such that the feasibility space $\mathcal{X}(x_1, \hat{\beta}, \xi^T)$ is not empty. In order to have an exact solution to problem (4), we need to expand $\hat{\Xi}$ with significant scenarios (and consequently the set of recourse variable copies $x_{t,\xi^t}^r$), until it can be proven that adding another scenario does not improve the optimal objective value of $\mathcal{MP}(\hat{\Xi})$, at which point its optimal objective value would be equal to $\nu^{2S-LDR}$. This implies that, at convergence, solving $\mathcal{MP}(\hat{\Xi})$ returns an optimal solution $(\hat{x}_1, \hat{\beta}, \hat{\eta})$ such that $\hat{\eta}$ accurately measures the worst-case second-stage cost over the complete uncertainty set $\Xi$. To check whether the convergence criterion is satisfied, we solve the adversarial problem for a given first-stage solution $\hat{\beta}$, which results in the following subproblem:

$$\mathcal{SP}(\hat{\beta}) = \max_{\xi^T \in \Xi} \left\{ \sum_{t \in [2,T]} c_t^{\top} \Phi_t (\xi^t)^{\top} \hat{\beta}_t + \min_{x^r \in \mathcal{X}(x_1, \hat{\beta}, \xi^T)} \sum_{t \in [2,T]} c_t^{\top} x_t^r \right\}. \qquad (6)$$

If $\hat{\eta} = \mathcal{SP}(\hat{\beta})$, then $\hat{\eta}$ exactly measures the worst-case cost of the second-stage problem and $(\hat{x}_1, \hat{\beta}, \hat{\eta})$ is an optimal solution to problem (4), i.e., $\nu^{2S-LDR} = c_1^\top \hat{x}_1 + \hat{\eta}$. Otherwise, let $\hat{\xi}^T$ be the optimal solution of subproblem (6) if it is feasible, or a scenario such that $\mathcal{X}(x_1, \hat{\beta}, \hat{\xi}^T)$ is an empty set. We create new variables $x_{t,\hat{\xi}^t}^r, t \in [2, T]$ and update the set of significant scenarios with $\hat{\Xi} = \hat{\Xi} \cup \{\hat{\xi}^T\}$, and accordingly update (5c) and (5d). If subproblem (6) is feasible, then we update $\Xi^{\text{FEAS}} = \Xi^{\text{FEAS}} \cup \{\hat{\xi}^T\}$ and (5b) as well. In this case, constraints (5b)-(5d) make up the optimality cuts. If subproblem (6) is infeasible, then (5c)-(5d) act as feasibility cuts. We repeat solving the master problem and the subproblem in a cutting-plane fashion and add the appropriate cuts, until $\hat{\Xi}$ includes all the significant scenarios and $\hat{\eta} = \mathcal{SP}(\hat{\beta})$. The following proposition provides a sufficient condition for finite convergence of the C&CG algorithm.

**PROPOSITION 1.** Assume that the uncertainty set $\Xi$ is either a finite discrete set or a compact polyhedron, and the basis functions $\Phi_t(\xi^t)$ are affine in $\xi^t$. Then the C&CG algorithm converges to $\nu^{2S-LDR}$ of an MSARO with continuous and fixed recourse in a finite number of iterations.

*Proof* Denote by $f(x_1, \xi^T, \beta)$, the objective function of the inner maximization in the two-stage problem (4). Then, we can rewrite problem (4) as:

$$\nu^{2S-LDR} = \min \left\{ c_1^\top x_1 + \max_{\xi^T \in \Xi} f(x_1, \xi^T, \beta) \mid x_1 \in X_1, \beta_t \in \mathbb{R}^{K_t}, t \in [2, T] \right\}.$$

First, we show that $f(x_1, \xi^T, \beta)$ is convex in $\xi^T$ for given $x_1 \in X_1$. For $\xi_1^T, \xi_2^T \in \Xi$ and $\lambda \in [0, 1]$, the following (in)equalities hold:

$$\lambda f(x_1, \xi_1^T, \beta) + (1 - \lambda) f(x_1, \xi_2^T, \beta) =$$
$$ \text{s.t.} \quad x_1 \in X_1 \qquad (6b) $$

$$ \beta_t \in \mathbb{R}^{K_t} \quad t \in [2, T]. \qquad (6c) $$

As reviewed in Section 2, the 2ARO problem (6) can either be solved by means of approximation (e.g., K-adaptability (Subramanyam et al. 2020), uncertainty set partitioning (Postek and Hertog 2016, Bertsimas and Dunning 2016), Neur2RO (Dumouchelle et al. 2023)) or exactly, most notably via the commonly used constraint-and-column generation (C&CG) method, initially proposed by Zeng and Zhao (2013), which we detail next.

The C&CG method draws on the fact that not all realizations in $\Xi^T$ contribute to the worst-case objective value. It then strives to identify necessary realizations by starting from a smaller uncertainty set and gradually expanding it. This leads to the generation of new columns and constraints, respectively corresponding to recourse variables and second-stage constraints for the newly identified realization. More specifically, consider a relaxation of problem (6) where instead of the uncertainty set $\Xi^T$, a potentially empty subset $\hat{\Xi} \subseteq \Xi^T$ is used to obtain the following master problem, which we denote by $MP(\hat{\Xi})$:

$$ \min \quad c_1^\top x_1 + \eta \qquad (7a) $$

$$ \text{s.t.} \quad \eta \ge \sum_{t \in [2, T]} \left( c_t^s (\xi^t)^\top \Phi_t (\xi^t) \beta_t + c_t^r (\xi^t)^\top x_{t, \xi^T}^r \right) \qquad \xi^T \in \hat{\Xi}^{\text{FEAS}} \qquad (7b) $$

$$ A_t^r(\xi^t)x_{t,\xi^T}^r + A_t^s(\xi^t)\Phi_t(\xi^t)\beta_t + B_t^s(\xi^t)x_1^s \le b_t(\xi^t) \qquad t=2, \xi^T \in \hat{\Xi} \qquad (7c) $$

$$ A_t^r(\xi^t)x_{t,\xi^T}^r + A_t^s(\xi^t)\Phi_t(\xi^t)\beta_t + B_t^s(\xi^t)\Phi_{t-1}(\xi^{t-1})\beta_{t-1} \le b_t(\xi^t) \qquad t \in [3, T], \xi^T \in \hat{\Xi} \qquad (7d) $$

$$ (\Phi_t(\xi^t)\beta_t, x_{t,\xi^T}^r) \in X_t(\xi^t) \qquad t \in [2, T], \xi^T \in \hat{\Xi} \qquad (7e) $$

$$ x_1 \in X_1, \quad \eta \in \mathbb{R} \qquad (7f) $$

$$ \beta_t \in \mathbb{R}^{K_t} \qquad t \in [2, T], \qquad (7g) $$

where $\hat{\Xi}^{\text{FEAS}} \subseteq \hat{\Xi}$ is the set of identified necessary realizations $\xi^T$ for which there exists a feasible first-stage solution $\hat{\beta}$ such that the feasibility space $\mathcal{X}(x_1^s, \hat{\beta}, \xi^T)$ is nonempty. If the recourse variables are all continuous (i.e., $p_t^j = 0, t \in [2, T]$), then (7) is a linear program, otherwise it is a mixed-integer linear program. We remark that if $\hat{\Xi}^{\text{FEAS}}$ is empty at initialization, a valid lower bound on the $\eta$ variable can be added to the model to avoid unboundedness. If $MP(\hat{\Xi})$ is infeasible for any $\hat{\Xi} \subseteq \Xi^T$, the 2ARO model (6) is proven to be infeasible. Next, we consider the cases where $MP(\hat{\Xi})$ is feasible and bounded.

The optimal objective value of $MP(\hat{\Xi})$ is a lower bound on $\nu^{2S-LDR}$, the optimal objective value of (6). To obtain an exact solution to problem (6), we may need to gradually expand $\hat{\Xi}$ with necessary realizations (and consequently the set of recourse variable copies $x_{t,\xi^T}^r$). At convergence,
$$ \sum_{t \in [2,T]} c_t^{\mathrm{s}^{\top}} (\lambda \Phi(\xi_1^t)^{\top} + (1-\lambda) \Phi(\xi_2^t)^{\top}) \beta_t + \\ \lambda \left( \min_{x^r \in \mathcal{X}(x_1, \beta, \xi_1^T)} \sum_{t \in [2,T]} c_t^{x^r} x_t^r + (1-\lambda) \left( \min_{x^r \in \mathcal{X}(x_1, \beta, \xi_2^T)} \sum_{t \in [2,T]} c_t^{x^r} x_t^r \right) \right) = \qquad (7a) $$

$$ \sum_{t \in [2,T]} c_t^{\mathrm{s}^{\top}} \Phi(\lambda \xi_1^t + (1-\lambda) \xi_2^t)^{\top} \beta_t + \\ \lambda \left( \min_{x^r \in \mathcal{X}(x_1, \beta, \xi_1^T)} \sum_{t \in [2,T]} c_t^{x^r} x_t^r + (1-\lambda) \left( \min_{x^r \in \mathcal{X}(x_1, \beta, \xi_2^T)} \sum_{t \in [2,T]} c_t^{x^r} x_t^r \right) \right) \geq \qquad (7b) $$

$$ \sum_{t \in [2,T]} c_t^{\mathrm{s}^{\top}} \Phi(\lambda \xi_1^t + (1-\lambda) \xi_2^t)^{\top} \beta_t + \min_{x^r \in \mathcal{X}(x_1, \beta, \lambda \xi_1^T (1-\lambda) \xi_2^T)} \sum_{t \in [2,T]} c_t^{x^r} x_t^r = f(x_1, \lambda \xi_1^T + (1-\lambda) \xi_2^T, \beta). $$

Equality (7a) holds because $\Phi(\xi^t)$ is affine. Inequality (7b) follows from the convexity of the linear minimization problem (by the fixed recourse assumption), since $\xi$ appears only on the right-hand-sides of constraints and $x_t^r$ are continuous (by the continuous recourse assumption).

Now, let $\Xi^{\mathrm{EXT}}$ be the set of extreme points of $\Xi$ if it is a polyhedron, or that of its convex-hull $\mathrm{conv}(\Xi)$ if it is a finite discrete set. In maximization of a convex function over a compact polyhedral set, there is an optimal solution that is an extreme point (Hendrix et al. 2010). Then, the two-stage problem becomes:

$$ \nu^{2\mathrm{S}-\mathrm{LDR}} = \min \{c_1^\top x_1 + \eta \mid \eta \ge f(x_1, \xi^T, \beta), \xi^T \in \Xi^{\mathrm{EXT}}, x_1 \in X_1, \beta_t \in \mathbb{R}^{K_t}, t \in [2,T]\}. $$

C&CG is then the process of gradually adding constraints for each extreme point. Because $|\Xi^{\mathrm{EXT}}| < +\infty$, iterating over all extreme points takes finitely many steps, which concludes the proof. $\square$

In case the conditions of Proposition 1 are not satisfied, the constraint-and-column generation algorithm can still converge asymptotically to the optimal solution of problem (4). This is because the lower bound (optimal objective value of $\mathcal{MP}(\Xi)$) improves with each iteration, and under the Assumptions 1 and 2, the optimal value of problem (4) is finite. Consequently, the algorithm will asymptotically approach this optimal value.

**REMARK 1.** With continuous recourse (i.e., $n_t^i = 0, t \in [2, T]$), the subproblem (6) can be reformulated as a bilinear monolithic problem, using LP duality or Karush-Kuhn-Tucker conditions (Zeng and Zhao 2013). Let $\pi^S$ and $\pi^R$ be the LP dual variables associated with the state and recourse constraints in $\mathcal{X}(x_1, \hat{\beta}_t, \xi^T)$, in the inner minimization problem of (6), respectively. Then, subproblem (6) with continuous recourse can be modeled as follows:

$$ \max \sum_{t \in [2,T]} c_t^{\mathrm{s}^{\top}} \Phi_t (\xi^t)^{\top} \hat{\beta}_t + c_t^{\mathrm{x}^{\top}} x_t^{\mathrm{x}} \qquad (8a) $$

$$ \text{s.t.} \quad A_t^{\mathrm{x}} x_t^{\mathrm{x}} + A_t^{\mathrm{s}} \Phi_t (\xi^t)^{\top} \hat{\beta}_t + B_t^{\mathrm{s}} \Phi_{t-1} (\xi^{t-1})^{\top} \hat{\beta}_{t-1} - b_t (\xi^t) \leq 0 \quad t \in [2, T] \qquad (8b) $$

$$ D_t^{\mathrm{s}} \Phi_t (\xi^t)^{\top} \hat{\beta}_t + D_t^{\mathrm{x}} x_t^{\mathrm{x}} - d_t (\xi^t) \leq 0 \quad t \in [2, T] \qquad (8c) $$
solving $\mathcal{MP}(\Xi)$ should return an optimal solution $(\hat{x}_1, \hat{\beta}, \hat{\eta})$ such that $\hat{\eta}$ accurately measures the worst-case second-stage cost over the complete uncertainty set $\Xi^T$ (or conclude infeasibility of (6)). To check whether this convergence criterion is satisfied, we solve the adversarial problem for a given first-stage solution $(\hat{\beta}, \hat{x}_1^s)$, which results in the following subproblem:

$$ \mathcal{SP}(\hat{\beta}, \hat{x}_1^s) := \max_{\xi^T \in \Xi^T} \left\{ \sum_{t \in [2,T]} c_t^s(\xi^t)^\top \Phi_t(\xi^t) \hat{\beta}_t + \min_{x^r \in \mathcal{X}(\hat{x}_1^s, \hat{\beta}, \xi^T)} \sum_{t \in [2,T]} c_t^r(\xi^t)^\top x_t^r \right\}. \quad (8) $$

If $\hat{\eta} = \mathcal{SP}(\hat{\beta}, \hat{x}_1^s)$, then $\hat{\eta}$ exactly measures the worst-case cost of the second-stage problem and $(\hat{x}_1, \hat{\beta}, \hat{\eta})$ is an optimal solution to problem (6), i.e., $\nu^{2S-LDR} = c_1^\top \hat{x}_1 + \hat{\eta}$. Otherwise, let $\hat{\xi}^T$ be the optimal solution of subproblem (8) if it is feasible, or a scenario such that $\mathcal{X}(\hat{x}_1^s, \hat{\beta}, \hat{\xi}^T)$ is an empty set. We create new variables $x_t^r, t \in [2, T]$ and update the set of necessary realizations with $\hat{\Xi} = \hat{\Xi} \cup \{\hat{\xi}^T\}$, and accordingly update (7c)-(7e). If subproblem (8) is feasible, then we update $\hat{\Xi}^{FEAS} = \hat{\Xi}^{FEAS} \cup \{\hat{\xi}^T\}$ and (7b) as well. In this case, constraints (7b)-(7e) make up the optimality cuts. If subproblem (8) is infeasible, then (7c)-(7e) act as feasibility cuts. We repeat solving the master problem and the subproblem in a cutting-plane fashion and add the appropriate cuts until $\hat{\Xi}$ includes all the necessary realizations and $\hat{\eta} = \mathcal{SP}(\hat{\beta}, \hat{x}_1^s)$, or $\mathcal{MP}(\Xi)$ becomes infeasible.

In general, subproblem (8) can be numerically challenging to solve especially if the inner minimization problem contains integer variables and/or non-linearities in $\xi^T$. However, there are certain practical cases in which (8) can be reformulated as a mixed-integer linear program and then directly be given to an optimization solver. In the context of our subproblem, one notable example of this is provided in Remark 1. Further, extending this case to allow integer recourse variables ($p_t^i \neq 0, t \in [2, T]$), while keeping the other assumptions in Remark 1, Zhao and Zeng (2012) proposed a nested constraint-and-column generation algorithm. In addition, we note that it is not necessary to solve the subproblem exactly at each iteration; the solution process can be stopped as soon as a violated cut for the master problem is identified. Recent studies have also considered the use of neural networks in order to approximate the inner minimization problem (Dumouchelle et al. 2023), especially to handle integer recourse variables and non-linearities.

**REMARK 1.** Consider an MSARO where the uncertainty set is a polytope, the basis functions $\Phi_t(\xi^t)$ are chosen to be affine in $\xi^t$ for all $t \in [2, T]$, all the recourse variables are continuous (i.e., $p_t^i = 0, t \in [2, T]$) and we have fixed parameters associated with the state variables, i.e., $c_t^s(\xi^t) = c_t^s, A_t^s(\xi^t) = A_t^s, B_t^s(\xi^t) = B_t^s, D_t^s(\xi^t) = D_t^s$ for all $t \in [2, T]$ and $\xi^t \in \Xi^t$. In this case, using linear programming duality or the KKT optimality conditions (under a further relatively complete recourse assumption for the 2ARO problem (6)) yields a monolithic bilinear continuous optimization problem, which can further be linearized using big-M constraints to obtain a mixed-integer linear program (see Appendix C.1 for details).
$$A_t^{\mathrm{r}^\top} \pi_t^\mathrm{s} + D_t^{\mathrm{r}^\top} \pi_t^\mathrm{R} = c_t^\mathrm{r} \qquad t \in [2, T] \qquad (8d)$$

$$(A_t^\mathrm{r} x_t^\mathrm{r} + A_t^\mathrm{s} \Phi_t (\xi^t)^\top \hat{\beta}_t + B_t^\mathrm{s} \Phi_{t-1} (\xi^{t-1})^\top \hat{\beta}_{t-1} - b_t(\xi^t))^{\top} \pi_t^\mathrm{s} = 0 \qquad t \in [2, T] \qquad (8e)$$

$$(D_t^\mathrm{s} \Phi_t (\xi^t)^\top \hat{\beta}_t + D_t^\mathrm{r} x_t^\mathrm{r} - d_t(\xi^t))^{\top} \pi_t^\mathrm{R} = 0 \qquad t \in [2, T] \qquad (8f)$$

$$\xi^T \in \Xi, \quad x_t^\mathrm{r} \in \mathbb{R}^{p_t}, \quad \pi_t^\mathrm{s} \in \mathbb{R}_-^{m_t^s}, \quad \pi_t^\mathrm{R} \in \mathbb{R}_-^{m_t^R} \qquad t \in [2, T]. \qquad (8g)$$

Inequalities (8b)-(8c) are the primal feasibility constraints at $\hat{\beta}_t$, (8d) the dual feasibility constraints, and (8e)-(8f) are the complementary slackness constraints. The latter include bilinear terms, as basis functions $\Phi_t(\xi^t)$ are functions of $\xi^T$ which are decision variables in (8). They can be linearized with the addition of binary decision variables via the so-called big-M constraints. A detailed example of building a mixed-integer linear model for a multistage location-allocation problem is provided in Appendix A.

## 2.3. Two-stage Piecewise-constant Decision Rules

In this section we study the application of two-stage decision rules in another special case of MSAROs, where each state variable $x_{ti}^s(\xi^t) \in \mathbb{Z}, i \in [q_t], t \in [2, T]$ is a bounded integer with a given domain $[\underline{\kappa}_{ti}, \overline{\kappa}_{ti}]$. We define the two-stage piecewise-constant decision rule (PCDR) as follows:

$$\Theta_{ti}(\xi^t, \beta_t) = \begin{cases} \underline{\kappa}_{ti} & \Upsilon_{ti}(\xi^t, \beta_t^i) \in \mathcal{K}_{t1} \\ \underline{\kappa}_{ti} + 1 & \Upsilon_{ti}(\xi^t, \beta_t^i) \in \mathcal{K}_{t2} \\ \vdots & \vdots \\ \overline{\kappa}_{ti} & \Upsilon_{ti}(\xi^t, \beta_t^i) \in \mathcal{K}_{tJ}, \end{cases}$$

where $\mathcal{K}_{tj} \subset \mathbb{R}, j \in [J]$ are disjoint sets with $\cup_{j \in [J]} \mathcal{K}_{tj} = \mathbb{R}$, and $\Upsilon_{ti}(\xi^t, \beta_t^i) : \mathbb{R}^{l_t} \times \mathbb{R}^{K_t} \to \mathbb{R}$ are functions defining the policy. We remark that the PCDR is presented for each individual variable $x_{ti}^s(\xi^t)$, and has its own decision variables $\beta_t^i$. Semantically, the PCDR partitions the set $\mathbb{R}$ into $J$ subsets, and then assigns an integer value to each partition.

A special case of PCDRs can be defined by restricting the form of all mappings $\Upsilon_{ti}(\xi^t, \beta_t^i)$ to be a linear function of the history. In the following we show that such a decision rule results in a model that is structurally very similar to (4), thus it can be solved via the C&CG method. Let $\Upsilon_{ti}(\xi^t, \beta_t^i) = \Upsilon_{ti}^A(\xi^t)^\top \beta_t^i$ be an affine function of basis functions $\Upsilon_{ti}^A(\xi^t)$, and $\mathcal{K}_{tj} = (a_{tj}, b_{tj}], j \in [J]$ be intervals (in practice $\mathbb{R}$ is replaced by $[-M, M]$ for $M$ sufficiently large). Then problem (3) becomes:

$$\nu^{2S-\mathrm{PCDR}} = \min_{x_1, \beta} \left\{ c_1^\top x_1 + \mathcal{SP}^{\mathrm{PCDR}}(x_1, \beta) \mid x_1 \in X_1, \ \beta_t \in \mathbb{R}^{q_t \times K_t}, t \in [2, T] \right\},$$

where:

$$\mathcal{SP}^{\mathrm{PCDR}}(x_1, \beta) = \max_{\xi^T \in \Xi} \sum_{t \in [T]} \sum_{i \in [q_t]} c_{ti}^{\mathrm{s}} x_{ti}^{\mathrm{s}} + \min_{x^\mathrm{r} \in \mathcal{X}(x_1, \beta, \xi^T)} \sum_{t \in [2, T]} c_t^{\mathrm{r}^\top} x_i^{\mathrm{r}} \qquad (9a)$$
We further remark that the presence of decision-rule design variables $\beta$ may make the master problem $MP(\hat{\Xi})$ numerically challenging to solve, compared to traditional use cases of C&CG. In that regard, recently proposed enhancement ideas can be employed, such as the inexact C&CG algorithm proposed by Tsang et al. (2023). In their framework, master problems are solved to a given relative optimality gap which gradually reduces to zero over the course of the algorithm. Additionally, it is possible to design the basis functions $\Phi_t(\cdot)$ to use a small information basis rather than the entire history $\xi^t$ so that the number of design variables is reduced. This would be especially helpful for problems with larger number of decision stages.

Another consideration concerning the exact solution of (6) using the C&CG algorithm is its convergence. Indeed, the C&CG may not have finite convergence in general. If $\Xi^T$ is a finite set, then the finite convergence is straightforward, otherwise more conditions are needed to ensure this property. Zeng and Zhao (2013) proved finite convergence, considering a (bounded) polyhedral uncertainty set, for problems with only right-hand-side uncertainty represented as an affine function of uncertain parameters. However, in problem (6), basis functions $\Phi_t(\xi^t)$ appear as coefficients of $\beta_t$, both in the objective function (6a) (in the term $c_t^s(\xi^t)^\top \Phi_t(\xi^t)$) and the second-stage constraints (terms $A_t^s(\xi^t)\Phi_t(\xi^t)$, $B_t^s(\xi^t)\Phi_{t-1}(\xi^{t-1})$ and $D_t^s(\xi^t)\Phi_t(\xi^t)$). The following proposition provides a sufficient condition for finite convergence of the C&CG algorithm (proof is given in Appendix B) in this more general context.

**PROPOSITION 1.** Consider an MSARO with only right-hand-side uncertainty, continuous recourse, and (bounded) polyhedral uncertainty set. If the basis functions $\Phi_t(\xi^t)$ are chosen to be affine in $\xi^t$ for all $t \in [2, T]$, the C&CG algorithm converges to $\nu^{2S-LDR}$ in a finite number of iterations.

In case the conditions of Proposition 1 are not satisfied, the C&CG algorithm still converges but asymptotically to an optimal solution of problem (6) if it is feasible (since it is bounded under the boundedness assumption imposed on the original MSARO problem).

Lastly, we note that in the case of continuous recourse, a linear decision rule can be applied to the recourse decision variables as well resulting in the LDR approach proposed by Ben-Tal et al. (2004). Since recourse variables $x_t^s(\xi^t)$ in (6) are fully adjustable, it immediately follows that

$$ \nu^* \le \nu^{2S-LDR} \le \nu^{LDR}, $$

where $\nu^{LDR}$ refers to the bound obtained from the commonly used LDR approach.

### 3.3. Two-stage Piecewise-constant Decision Rules for MSAROs with Integer State Variables

In this section, we study the application of two-stage decision rules in another special case of MSAROs, where each state variable $x_{ti}^s(\xi^t) \in \mathbb{Z}, i \in [q_t], t \in [2, T]$ is a bounded integer with a given
$$ \text{s.t.} \quad \sum_{j \in [J]} (\kappa_{ti} + j - 1) v_{tij} = x_{ti}^s \qquad t \in [2, T], i \in [q_t] \qquad (9b) $$

$$ \sum_{j \in [J]} \omega_{tij} = \Upsilon_{ti}^{\mathrm{A}} (\boldsymbol{\xi}^t)^{\top} \beta_t^i \qquad t \in [2, T], i \in [q_t] \qquad (9c) $$

$$ (a_{tj} + \epsilon) v_{tij} \le \omega_{tij} \le b_{tj} v_{tij} \qquad t \in [2, T], i \in [q_t], j \in [J] \qquad (9d) $$

$$ \sum_{j \in [J]} v_{tij} = 1 \qquad t \in [2, T], i \in [q_t] \qquad (9e) $$

$$ v_{tij} \in \{0, 1\} \qquad t \in [2, T], i \in [q_t], j \in [J]. \qquad (9f) $$

In the outer maximization problem, the PCDR is modeled using the auxiliary variables $v_{tij}, \omega_{tij}$ and constraints (9b)-(9f) ($\epsilon$ is added to lower bound on $\omega_{tij}$ to ensure that partitions $K_{tj}$ are disjoint). Note that, similar to Remark 1, if recourse decision variables $x_t^r, t \in [2, T]$ are continuous, we can reformulate problem (9) as a monolithic mixed-integer linear program. Furthermore, applying the same arguments discussed in Section 2.2 shows that the C&CG algorithm will converge asymptotically to the optimal solution of $\nu^{2S-PCDR}$.

### 3. Dual Bounding

In this section, our goal is to introduce a new dual problem that provides a lower bound for MSARO problems with mixed-integer recourse. Note that, the linear programming (LP) dual of continuous MSARO problems has been studied before (Kuhn et al. 2011), showing that even then, establishing strong duality is not straightforward. Recall that a random variable $\xi \in \mathbb{R}^n$ is characterized by its probability space $(\mathbb{R}^n, \mathcal{F}, \mathbb{P})$, in which $\mathcal{F}$ is a collection of subsets of $\mathbb{R}^n$ and $\mathbb{P}: \mathcal{F} \to [0, 1]$ is a probability measure that assigns a value between 0 and 1 to each element of $F \in \mathcal{F}$. Then, we define the support of $\xi$ as the smallest closed set $\Xi$ in the Borel $\sigma$-algebra such that $\mathbb{P}(\Xi) = 1$. Obtaining a strong LP dual is challenging due to a fundamental aspect of MSARO: the lack of knowledge about the true probability distribution and therefore $\mathbb{P}$. To see the importance of $\mathbb{P}$ and its impact on the dual, we formulate the LP dual of model (1), with no integer variables (only to demonstrate the concept) for an *assigned* probability distribution. Denote by $\pi^{(1b)}(\xi^T)$, $\pi_t^{(1c)}(\xi^t)$, $\pi_t^{(1d)}(\xi^t)$, the LP dual variables associated with constraints (1b), (1c) and (1d), respectively. Let $p: \Xi^t \to \mathbb{R}^+$ be the probability density function of $\xi^t$ with respect to $\mathbb{P}$. First we multiply constraints (1b)-(1d) by their associated $p(\xi^t)$, and then apply the standard linear programming duality theory. The LP dual of continuous MSARO for $\mathbb{P}$ is:

$$ \begin{aligned} \mathcal{L}^{\mathrm{LP}}(\mathbb{P}) = & \\ \max \quad & \mathbb{E}_{\xi^T \sim \mathbb{P}} \sum_{t \in [T]} [b_t(\xi^t)^{\top} \pi_t^{(1c)}(\xi^t) + d_t(\xi^t)^{\top} \pi_t^{(1d)}(\xi^t)] \\ \text{s.t.} \quad & \mathbb{E}_{\xi^T \sim \mathbb{P}} [\pi^{(1b)}(\xi^T)] = 1 \end{aligned} $$
domain $[\underline{\kappa}_{ti}, \bar{\kappa}_{ti}]$ where boundedness follows from Assumption 3. We enforce the two-stage piecewise-constant decision rule (PCDR) by using $x_{ti}^s(\xi^t) = \Theta_{ti}(\xi^t, \beta_t^i)$ where

$$ \Theta_{ti}(\xi^t, \beta_t^i) = \begin{cases} \underline{\kappa}_{ti} & \Upsilon_{ti}(\xi^t, \beta_t^i) \in \mathcal{K}_{t1}^i \\ \underline{\kappa}_{ti} + 1 & \Upsilon_{ti}(\xi^t, \beta_t^i) \in \mathcal{K}_{t2}^i \\ \vdots & \vdots \\ \bar{\kappa}_{ti} & \Upsilon_{ti}(\xi^t, \beta_t^i) \in \mathcal{K}_{tJ_i}^i, \end{cases} \quad (9) $$

where $\mathcal{K}_{tj}^i \subset [-1, 1], j \in [J_i]$ are disjoint sets with $\cup_{j \in [J_i]} \mathcal{K}_{tj}^i = [-1, 1]$, and $\Upsilon_{ti}(\xi^t, \beta_t^i) : \mathbb{R}^{\ell_t} \times \mathbb{R}^{K_t} \to [-1, 1]$ are functions defining the policy for $t \in [2, T]$ and $i \in [q_t]$. Semantically, the PCDR partitions the interval $[-1, 1]$ into subsets, and then assigns an integer value to each partition.

A special case of PCDRs can be defined by restricting the form of all mappings $\Upsilon_{ti}(\xi^t, \beta_t^i)$ to be a linear function of the decision rule design variables $\beta$. In the following, we show that such a decision rule results in a model that is structurally very similar to (6), thus it is amenable to the C&CG method. Let, without loss of generality, $\mathcal{K}_{t1}^i = [a_{t1}^i, b_{t1}^i]$, and $\mathcal{K}_{tj}^i = (a_{tj}^i, b_{tj}^i], j \in [J_i] \setminus \{1\}$ be intervals with $a_{t1}^i = -1$, $b_{t1}^i = 1$, and $a_{tj}^i = b_{t,j-1}^i$ for $j \in [2, J_i]$. Let further, for $t \in [2, T]$, $\Upsilon_{ti}(\xi^t, \beta_t^i) = \hat{\Upsilon}_{ti}(\xi^t)^\top \beta_t^i$ be chosen as an affine function of basis functions $\hat{\Upsilon}_{ti}(\xi^t)$. This implies, in particular, that $-1 \le \hat{\Upsilon}_{ti}(\xi^t)^\top \beta_t^i \le 1$ which is enforced through robust constraints. Then problem (4) becomes:

$$ \nu^{2S-PCDR} := \min c_1^\top x_1 + SP^{PCDR}(\beta, x_1^s) \quad (10a) $$

$$ \text{s.t.} \quad x_1 \in X_1 \quad (10b) $$

$$ \beta_t \in \mathbb{R}^{K_t} \qquad t \in [2, T] \quad (10c) $$

$$ -1 \le \hat{\Upsilon}_{ti}(\xi^t)^\top \beta_t^i \le 1 \qquad t \in [2, T], i \in [q_t], \xi^T \in \Xi^T \quad (10d) $$

where:

$$ SP^{PCDR}(\beta, x_1^s) := \max_{\xi^T \in \Xi^T} \sum_{t \in [2, T]} c_t^s(\xi^t)^\top x_t^s + \min_{x^r \in \mathcal{X}(x_1^s, \beta, \xi^T)} \sum_{t \in [2, T]} c_t^r(\xi^t)^\top x^r \quad (11a) $$

$$ \text{s.t.} \quad \sum_{j \in [J_i]} (\underline{\kappa}_{ti} + j - 1) v_{tij} = x_{ti}^s \qquad t \in [2, T], i \in [q_t] \quad (11b) $$

$$ \sum_{j \in [J_i]} \omega_{tij} = \hat{\Upsilon}_{ti}(\xi^t)^\top \beta_t^i \qquad t \in [2, T], i \in [q_t] \quad (11c) $$

$$ (a_{tj}^i + \epsilon_j) v_{tij} \le \omega_{tij} \le b_{tj}^i v_{tij} \qquad t \in [2, T], i \in [q_t], j \in [J_i] \quad (11d) $$

$$ \sum_{j \in [J_i]} v_{tij} = 1 \qquad t \in [2, T], i \in [q_t] \quad (11e) $$

$$ v_{tij} \in \{0, 1\} \qquad t \in [2, T], i \in [q_t], j \in [J_i]. \quad (11f) $$

with $\epsilon_1 = 0$. The PCDR is modeled using the auxiliary variables $v_{tij}, \omega_{tij}$ and constraints (11b)-(11f). Variables $v$ determine in which interval the quantity $\hat{\Upsilon}_{ti}(\xi^t)^\top \beta_t^i$ falls, in accordance with the
$$ \mathbb{E}_{\xi^T \sim \mathbb{P}} [-c_1^\top \pi^{(1b)}(\xi^T) + B_2(\xi^2)^\top \pi_2^{(1c)}(\xi^2)] + D_1 \pi_1^{(1d)} = \mathbf{0} $$

$$ \mathbb{E}_{\xi^T \sim \mathbb{P}} [-c_t(\xi^t)^\top \pi^{(1b)}(\xi^T) | \xi^t] + p(\xi^t)A_t(\xi^t)^\top \pi_t^{(1c)}(\xi^t) + p(\xi^{t+1})B_{t+1}(\xi^{t+1})^\top \pi_{t+1}^{(1c)}(\xi^{t+1}) + p(\xi^t)D_t(\xi^t)^\top \pi_t^{(1d)}(\xi^t) = \mathbf{0} \quad t \in [2, T-1], \xi^t \in \Xi^t $$

$$ -p(\xi^T)c_T(\xi^T)^\top \pi^{(1b)}(\xi^T) + p(\xi^T)A_T(\xi^T)^\top \pi_T^{(1c)}(\xi^T) + p(\xi^T)D_T(\xi^T)^\top \pi_T^{(1d)}(\xi^T) = \mathbf{0} $$

$$ \xi^T \in \Xi $$

$$ \pi^{(1b)}(\xi^T) \geq \mathbf{0} $$

$$ \xi^T \in \Xi $$

$$ \pi_t^{(1c)}(\xi^t) \leq \mathbf{0} $$

$$ t \in [2, T], \xi^t \in \Xi^t $$

$$ \pi_t^{(1d)}(\xi^t) \leq \mathbf{0} $$

$$ t \in [T], \xi^t \in \Xi^t. $$

Although we have used an assigned probability distribution $\mathbb{P}$ in the derivation of the LP dual bound, as mentioned, the underlying assumption in robust optimization is that we do not know the true probability distribution. As discussed by Kuhn et al. (2011), when this dual bound is used for robust optimization problems, the quality of the bound can be drastically different based on the choice of the assigned distribution $\mathbb{P}$. To obtain the best LP dual bound, denote by $\mathcal{P} = \{\mathbb{P} | p^\mathbb{P}(\xi^T) > 0, \xi^T \in \Xi\}$, an ambiguity set including all probability distributions with support $\Xi$, i.e., every $\mathbb{P} \in \mathcal{P}$ has a density function that assigns a strictly positive value to all $\xi^T \in \Xi$. The following distributionally robust optimization model returns the best lower bound over $\mathcal{P}$:

$$ \sup_{\mathbb{P} \in \mathcal{P}} \mathcal{L}^{\mathbb{P}}(\mathbb{P}). $$

This shows that, even in the continuous setting, there exists a trade-off between the quality of the LP bound from an assigned probability distribution and the computational effort of searching over the entire ambiguity set. To the best of our knowledge, such a lower bounding problem has not been studied before. In the following sections, we introduce the *nonanticipative* (NA) dual problem of an MSARO with *mixed-integer* recourse, which relies on the use of linear decision rules in a Lagrangian dual of problem (1).

## 3.1. Nonanticipative Dual of the MSARO

The NA dual is based on a reformulation of the MSARO problem where we create a copy of every decision variable for every realization and explicitly enforce nonanticipativity. To this end, we introduce the copy variables $y(\xi^T) = (y_1(\xi^T), \dots, y_T(\xi^T))$ as perfect information variables depending on the entire realization $\xi^T = (\xi_1, \dots, \xi_T)$. We can reformulate the MSARO problem (1) as:

$$ \min z \qquad (10a) $$
variables $\omega$, and the integer values assigned to the state variables. Here, $\epsilon_j$ is added to the lower bound on $\omega_{tij}$ to ensure that partitions $K_{tj}$ are disjoint. We remark that, it is possible to choose $\epsilon_j = 0$ for $j \in [J_i]$, in which case the intervals would intersect at their boundaries. In this case, whenever the quantity $\hat{\Upsilon}_{ti}(\xi^t)^\top \beta_t^i$ is a boundary point, the model allows assigning either one of the corresponding integer values to the associated state variable. The objective function then dictates that the solution leading to the worst objective value is chosen.

Similar to Remark 1, in certain special cases, we can reformulate problem (11) as a monolithic mixed-integer linear program.

REMARK 2. Consider an MSARO where the uncertainty set is a polytope, the basis functions $\hat{\Upsilon}_{ti}(\xi^t)$ are chosen to be affine in $\xi^t$ for all $t \in [2, T]$ and $i \in [q_t]$, all the recourse variables are continuous, and we have fixed parameters associated with the state variables, i.e., $c_t^s(\xi^t) = c_t^s, A_t^s(\xi^t) = A_t^s, B_t^s(\xi^t) = B_t^s, D_t^s(\xi^t) = D_t^s$ for all $t \in [2, T]$ and $\xi^t \in \Xi^t$. In this case, using linear programming duality yields a monolithic bilinear continuous optimization problem whose objective function involves products between variables $x_{it}^s$ and linear programming dual variables of the inner minimization problem. Thanks to constraints (11b), variables $x_{it}^s$ can be substituted for a weighted sum of binary variables $v$, as such these bilinear terms can be linearized using big-M constraints to obtain a mixed-integer linear program.

Further, in applying the C&CG method to model (10), robust constraints (10d) will appear in the master problem. If the uncertainty set is a polytope and the basis functions $\hat{\Upsilon}_{ti}(\xi^t)$ are chosen to be affine in $\xi^t$ for all $t \in [2, T]$ and $i \in [q_t]$ then these semi-infinite constraints can be reformulated as a finite set of linear constraints using classical robust optimization techniques based on linear programming duality. We also note that the arguments presented in Section 3.2 imply similarly that the C&CG algorithm converges asymptotically to the optimal solution of $v^{2S-PCDR}$ if it is feasible.

Lastly, the methods presented in Sections 3.2 and 3.3 can be combined to address MSAROs with mixed-integer state variables.

REMARK 3. For an MSARO with mixed-integer state variables, the application of linear and piecewise constant decision rules, given by equations (5) and (9), to the continuous and integer state variables, respectively, yields a 2ARO approximation. The resulting model is presented in detail in Appendix C.2. This model is similarly amenable to the C&CG method for exact solution but can also benefit from other 2ARO solution methods from the literature.
$$ \text{s.t.} \quad \sum_{t \in [T]} c_t (\boldsymbol{\xi}^t)^\top y_t (\boldsymbol{\xi}^T) \leq z \qquad \boldsymbol{\xi}^T \in \Xi \qquad (10b) $$

$$ A_t(\boldsymbol{\xi}^t)y_t(\boldsymbol{\xi}^T) + B_t(\boldsymbol{\xi}^t)y_{t-1}(\boldsymbol{\xi}^T) \leq b_t(\boldsymbol{\xi}^t) \qquad t \in [T], \boldsymbol{\xi}^T \in \Xi \qquad (10c) $$

$$ D_t(\boldsymbol{\xi}^t)y_t(\boldsymbol{\xi}^T) \leq d_t(\boldsymbol{\xi}^t) \qquad t \in [T], \boldsymbol{\xi}^T \in \Xi \qquad (10d) $$

$$ y_t(\boldsymbol{\xi}^T) = y_t(\boldsymbol{\xi}'^T) \qquad t \in [T], \boldsymbol{\xi}^T, \boldsymbol{\xi}'^T \in \Xi \text{ with } \boldsymbol{\xi}^t = \boldsymbol{\xi}'^t \qquad (10e) $$

$$ y_t(\boldsymbol{\xi}^T) \in \mathbb{R}^{n_t - n_t^i} \times \mathbb{Z}^{n_t^i} \qquad t \in [T], \boldsymbol{\xi}^T \in \Xi. \qquad (10f) $$

Constraints (10e) are nonanticipativity constraints which ensure that at stage t for every partial realization of $\boldsymbol{\xi}^T$, the decisions made are consistent (i.e., the decisions made in all realizations sharing the history $\boldsymbol{\xi}^t$ are the same). Before introducing the NA dual derived from (10), we present the following lemma, which we use to reformulate the nonanticipativity constraints.

LEMMA 1. For given $\mathbb{P} \in \mathcal{P}$, constraints (10e) are equivalent to the following:

$$ y_t(\boldsymbol{\xi}^T) = \mathbb{E}_{\boldsymbol{\xi}'^T \sim \mathbb{P}} [y_t(\boldsymbol{\xi}'^T) \mid \boldsymbol{\xi}'^t = \boldsymbol{\xi}^t], \quad t \in [T], \boldsymbol{\xi}^T \in \Xi. \qquad (11) $$

Proof For a fixed $t$ and $\boldsymbol{\xi}^t$, and every $\boldsymbol{\xi}'^T$ that shares the same history with $\boldsymbol{\xi}^T$ up to $t$ (including $\boldsymbol{\xi}^T$), multiply both sides of (10e) by $p(\boldsymbol{\xi}'^T)$. The corresponding constraint (11) is derived by integrating over all such $\boldsymbol{\xi}'^T$ scenarios, thus it is a valid equality for problem (10). Moreover, this constraint makes sure that the decision variables under all the scenarios sharing the same history take the same value, i.e., $\mathbb{E}_{\boldsymbol{\xi}'^T \sim \mathbb{P}} [y_t(\boldsymbol{\xi}'^T) \mid \boldsymbol{\xi}'^t = \boldsymbol{\xi}^t]$, which enforces the semantics of the nonanticipativity constraints (10e). $\square$

Denote by $Y(\boldsymbol{\xi}^T)$ the scenario feasibility space defined by constraints (10b)-(10d) and (10f) for all $t \in [T]$ and some $\boldsymbol{\xi}^T$. For an assigned probability distribution $\mathbb{P}$ and with given dual functionals $\lambda_t(\cdot) : \mathbb{R}^t \to \mathbb{R}^{n_t}$, relaxing the nonanticipativity constraints (11) results in the following nonanticipative Lagrangian relaxation problem:

$$ \mathcal{L}_{\text{LR}}^{\text{NA}}(\mathbb{P}, \lambda_1(\cdot), \dots, \lambda_T(\cdot)) = \min \quad z + \sum_{t \in [T]} \mathbb{E}_{\boldsymbol{\xi}'^T \sim \mathbb{P}} \left[ \lambda_t(\boldsymbol{\xi}^T)^\top \left( y_t(\boldsymbol{\xi}^T) - \mathbb{E}_{\boldsymbol{\xi}'^T \sim \mathbb{P}} [y_t(\boldsymbol{\xi}'^T) \mid \boldsymbol{\xi}'^t = \boldsymbol{\xi}^t] \right) \right] \qquad (12a) $$

$$ \text{s.t.} \quad (z, y_1(\boldsymbol{\xi}^T), \dots, y_T(\boldsymbol{\xi}^T)) \in Y(\boldsymbol{\xi}^T) \qquad \boldsymbol{\xi}^T \in \Xi. \qquad (12b) $$

$\mathcal{L}_{\text{LR}}^{\text{NA}}$ is a lower bound for the original MSARO problem. The NA Lagrangian dual problem aims to find the best bound among all such Lagrangian relaxation bounds:

$$ \mathcal{L}^{\text{NA}}(\mathbb{P}) = \max_{\lambda_1(\cdot), \dots, \lambda_T(\cdot)} \mathcal{L}_{\text{LR}}^{\text{NA}}(\mathbb{P}, \lambda_1(\cdot), \dots, \lambda_T(\cdot)). \qquad (13) $$

The following proposition shows that, regardless of the choice of $\mathbb{P}$, $\mathcal{L}^{\text{NA}}(\mathbb{P})$ is an exact dual bound for continuous MSARO problems.
## 4. Dual Bounding

In this section, we introduce a new dual problem that provides a lower bound for MSARO problems with mixed-integer recourse. Due to the existence of integer variables, we rely on Lagrangian duality techniques, where we create a Lagrangian relaxation and optimize over the Lagrangian dual multipliers. As the MSARO involves constraints corresponding to every realization of uncertainty, Lagrangian multipliers are functions of uncertainty, as such they are usually high (possibly infinite) dimensional. To overcome the difficulty in their optimization, we propose to apply decision rule restrictions to Lagrangian multipliers, leveraging ideas rooted in the MSP literature (Kuhn et al. 2011, Daryalal et al. 2024). In deriving a Lagrangian dual of the MSARO problem, we choose a probability distribution with the support as the uncertainty set, and use the associated density function to scale the constraints to be dualized, resulting in expectation terms in the objective function of the relaxation. Accordingly, we obtain a dual approximation of MSARO in the form of a two-stage stochastic program. This probability distribution-based approach has several benefits, most notably the possibility of leveraging state-of-the-art stochastic programming techniques to solve the dual problem. However, the quality of the resulting dual bound depends on the probability distribution used while developing the dual formulation, as previously observed by Kuhn et al. (2011) for MSAROs with continuous recourse. With the aim of identifying the strongest such dual bound, we formally pose a distribution optimization problem (akin to what was developed in (Hadjiyiannis et al. 2011)) and develop appropriate solution methods (tailored to the nature of the recourse variables) for the resulting distribution optimization problem. To the best of our knowledge, numerical solution of such a bounding problem and the quality of the obtained bounds have not been studied before.

In what follows, in Section 4.1, we introduce the nonanticipative reformulation of the MSARO problem and its Lagrangian dual. In Section 4.2, we present the restricted Lagrangian dual problem and define the associated distribution optimization problem for which we develop solution methods in Section 4.3. Lastly, in Section 4.4, we propose an alternative dual problem, which can be weaker in terms of the quality of the obtained bound but has computational advantages thanks to its decomposable structure.

### 4.1. Nonanticipative Dual of the MSARO

The nonanticipative (NA) dual is based on a reformulation of the MSARO problem where we create a copy of decision variables for every stage and every realization, and explicitly enforce nonanticipativity constraints. To this end, we introduce the copy variables $y(\xi^T) = (y_1(\xi^T), \dots, y_T(\xi^T))$ for all $\xi^T \in \Xi^T$ as perfect information variables depending on the entire realization $\xi^T = (\xi_1, \dots, \xi_T)$.
**PROPOSITION 2.** Let $\mathbb{P}$ be any probability distribution in $\mathcal{P}$. For continuous MSARO problems, $\mathcal{L}^{\text{NA}}(\mathbb{P})$ is a strong dual.

*Proof* The strength of a Lagrangian dual problem can be studied by its primal characterization, derived by Geoffrion (1974) for a mixed-integer linear optimization problem. Then, using standard Lagrangian duality theory (see, for instance, Wolsey and Nemhauser (1999)), the primal characterization of (13) is:

$$ \min z \qquad (14a) $$

$$ \text{s.t.} \quad (z, y_1(\xi^T), \dots, y_T(\xi^T)) \in \text{conv}(Y(\xi^T)) \qquad \xi^T \in \Xi \qquad (14b) $$

$$ y_t(\xi^T) = \mathbb{E}_{\xi'^T \sim \mathbb{P}} \left[ y_t(\xi'^T) \mid \xi'^t = \xi^t \right] \qquad t \in [T], \xi^T \in \Xi. \qquad (14c) $$

The result follows from the equivalence of (11) and (10e), as well as the fact that for continuous MSARO problems, we have that $\text{conv}(Y(\xi^T)) = Y(\xi^T)$ for all $\xi^T \in \Xi$. $\square$

We remark that our construction of the dual postulates that we multiply the constraints with a density function whose support is $\Xi$. Therefore, the strictly positive density property of $\mathbb{P} \in \mathcal{P}$ is necessary for the exactness of our formulation.

Proposition 2 suggests that the solution of problem (13) gives an exact dual bound for (10) (hence for (1)) if all decision variables are continuous. Although this bound is not exact in the case of mixed-integer recourse, the potential of leveraging the literature of multistage stochastic programming in achieving a dual bound for MSARO is quite appealing.

The objective function of the Lagrangian relaxation problem (12) contains (conditional) expectations of decision variables $y_t(\xi^T)$, which is computationally challenging. Because of our initial assumption of relatively complete recourse (Assumption 1), and the fact that $X_t(\xi^t)$ are bounded (Assumption 2), so is $Y(\xi^T)$ and thanks to Lemma 1 of Daryalal et al. (2021) we can replace the expectation term for $t$ in the objective function (12a) with:

$$ \mathbb{E}_{\xi^T \sim \mathbb{P}} \left[ \left( \lambda_t(\xi^T) - \mathbb{E}_{\xi'^T \sim \mathbb{P}} \left[ \lambda_t(\xi'^T) \mid \xi'^t = \xi^t \right] \right)^\top y_t(\xi^T) \right]. $$

For given $\lambda_t(\xi^T)$, this exchange allows us to compute the coefficients of $y_t(\xi^T)$ in the Lagrangian relaxation problem. Still, the optimal form of the dual functionals $\lambda_t(\cdot)$ need to be determined, making the problem (13) intractable. In the next section, we restrict these Lagrangian multipliers to follow LDRs and obtain a restricted dual problem with decision variables of finite dimension. Furthermore, this new dual problem is amenable to well-known solution techniques from the literature of two-stage stochastic programming, such as sample average approximation, which are designed to approximately solve a problem with expectation in the objective function.
We denote the decision variables by $y$ instead of $x$ used in previous sections to emphasize the fact that they are perfect information variables. We also note that we do not need to distinguish state and recourse variables in this section, thus vectors $y$ involve all the decisions. We can then obtain the NA reformulation of the MSARO problem (2) as:

$$ \min z \qquad (12a) $$

$$ \text{s.t.} \quad \sum_{t \in [T]} c_t (\xi^t)^\top y_t (\xi^T) \leq z \qquad \xi^T \in \Xi^T \qquad (12b) $$

$$ A_t(\xi^t)y_t(\xi^T) + B_t(\xi^t)y_{t-1}(\xi^T) \leq b_t(\xi^t) \qquad t \in [2, T], \xi^T \in \Xi^T \qquad (12c) $$

$$ D_t(\xi^t)y_t(\xi^T) \leq d_t(\xi^t) \qquad t \in [T], \xi^T \in \Xi^T \qquad (12d) $$

$$ y_t(\xi^T) = y_t(\xi'^T) \qquad t \in [T], \xi^T, \xi'^T \in \Xi^T \text{ with } \xi^T = \xi'^T \qquad (12e) $$

$$ y_t(\xi^T) \in \mathbb{R}^{n_t - n_t^1} \times \mathbb{Z}^{n_t^1} \qquad t \in [T], \xi^T \in \Xi^T. \qquad (12f) $$

Constraints (12e) are nonanticipativity constraints which ensure that at stage $t$ for every partial realization of $\xi^T$, the decisions made are consistent (i.e., the decisions made in all realizations sharing the history $\xi^t$ are the same).

**REMARK 4.** The nonanticipativity constraints are redundant for stage $T$, however, we include them in model (12) for notational convenience. In our implementation for the numerical results presented in Section 5, we exclude those redundant constraints.

In deriving a Lagrangian relaxation, we will first scale the constraints to be relaxed. As the constraints correspond to uncertainty realizations, we choose the scaling factors in such a way that they induce a probability measure over the support $\Xi^T$. To this end, we let $\mathbb{P}$ denote this probability measure such that $\mathbb{P}(\Xi^T) = 1$, which we interchangeably refer to as the probability distribution. We define $p^\mathbb{P}: \Xi^T \to \mathbb{R}^+$ as the associated density function. Lastly, we let $\mathcal{P}^> := \{\mathbb{P} \mid p^\mathbb{P}(\xi^T) > 0, \xi^T \in \Xi^T\}$, i.e., every $\mathbb{P} \in \mathcal{P}^>$ has a density function assigning a strictly positive value to all $\xi^T \in \Xi^T$.

Before introducing the NA dual derived from (12), we present the following lemma, which we use to reformulate the nonanticipativity constraints (proof in Appendix B).

**LEMMA 1.** For any $\mathbb{P} \in \mathcal{P}^>$, constraints (12e) are equivalent to the following:

$$ y_t(\xi^T) = \mathbb{E}_{\xi'^T \sim \mathbb{P}} [y_t(\xi'^T) \mid \xi'^t = \xi^t], \quad t \in [T], \xi^T \in \Xi^T. \qquad (13) $$

One advantage of this reformulation is the reduction in the number of constraints and in turn in the number of dual multipliers to be introduced. Let, to this end, $\lambda_t(\cdot): \mathbb{R}^{\ell^t} \to \mathbb{R}^{n_t}$ for $t \in [T]$, where $\mathbb{E}_{\xi^T \sim \mathbb{P}}[\lambda_t(\xi^T)] < +\infty$, be the dual functionals to be used in relaxing the nonanticipativity constraints (13). Let further the feasibility space for $\xi^T \in \Xi^T$ be $Y(\xi^T) := \{(z, y_1(\xi^T), \dots, y_T(\xi^T)) : $
## 3.2. Lagrangian Dual Decision Rules

The Lagrangian multipliers $\lambda_t(\xi^T)$ are infinite-dimensional. One way to obtain finite-dimensional decision variables is to enforce LDRs on the Lagrangian multipliers of problem (13) with a fixed $\mathbb{P}$. This is referred to as Lagrangian dual decision rules (LDDRs). For a set of pre-determined basis functions $\Psi_t: \Xi \to \mathbb{R}^{n_t \times K_t}$ and LDR decision variables $\alpha_t \in \mathbb{R}^{K_t}$, we restrict the form of $\lambda_t(\xi^T)$ at stage t as follows:

$$ \lambda_t(\xi^T) = \Psi_t(\xi^T)\alpha_t, $$

which gives us a restricted NA Lagrangian dual problem with respect to $\mathbb{P}$:

$$ \mathcal{L}_{R}^{\mathrm{NA}}(\mathbb{P}) = \max_{\alpha_1, \dots, \alpha_T} \mathcal{L}_{\mathrm{LR}}^{\mathrm{NA}}(\mathbb{P}, \Psi_1(\xi^T)\alpha_1, \dots, \Psi_T(\xi^T)\alpha_T). \quad (15) $$

Since problem (15) is a restriction of (13), we have $\mathcal{L}_{R}^{\mathrm{NA}}(\mathbb{P}) \le \mathcal{L}^{\mathrm{NA}}(\mathbb{P})$. An immediate question is whether by this restriction the bound gets weaker. More importantly, whether the choice of $\mathbb{P}$ affects the quality of the bound $\mathcal{L}_{R}^{\mathrm{NA}}(\mathbb{P})$. Using Lemma 2 in Daryalal et al. (2021), the primal characterization of $\mathcal{L}_{R}^{\mathrm{NA}}(\mathbb{P})$ is:

$$ \min z \tag{16a} $$

$$ \text{s.t.} \quad (z, y_1(\xi^T), \dots, y_T(\xi^T)) \in \text{conv}(Y(\xi^T)) \qquad \xi^T \in \Xi \tag{16b} $$

$$ \mathbb{E}_{\xi^T \sim \mathbb{P}} \left[ \Psi_t(\xi^T)^\top \left( y_t(\xi^T) - \mathbb{E}_{\xi'^T \sim \mathbb{P}} [y_t(\xi'^T) | \xi'^t = \xi^t] \right) \right] = 0 \qquad t \in [T] \tag{16c} $$

Comparing the problems (16) and (14), it is clear that for an assigned $\mathbb{P}$, (16) is a relaxation of (14), as constraints (16c) are an aggregation of (14c). Consequently, unlike $\mathcal{L}^{\mathrm{NA}}(\mathbb{P})$, even for continuous MSARO, $\mathcal{L}_{R}^{\mathrm{NA}}(\mathbb{P})$ is not a strong bound. Furthermore, as a result of constraints (16c), the strength of the restricted NA Lagrangian dual bound depends on the choice of probability distribution $\mathbb{P}$. Similar to the discussion on the LP dual at the beginning of the section, this observation motivates us to optimize over a set of probability distributions and find the best bound among all possible $\mathbb{P}$ and LDR variables $\alpha$. As such, we propose to solve the following DRO problem with ambiguity set $\mathcal{P}$:

$$ \nu_{R}^{\mathrm{NA}} = \sup_{\mathbb{P} \in \mathcal{P}} \mathcal{L}_{R}^{\mathrm{NA}}(\mathbb{P}). \tag{17} $$

We remark that (17) is a sup-max-min problem and is not in the familiar form of min-max-min that DRO problems are often represented with. However, as the notion of optimization over a family of probability distributions fits the field of DRO, we still refer to (17) as a DRO problem, with a slight abuse of terminology.

In linear and mixed-integer programming, strict inequalities such as the ones required for $\mathbb{P} \in \mathcal{P}$ (the strictly positive density property for all $\xi^T \in \Xi$) often cause numerical and theoretical
(12b) – (12d) and (12f)}. Then, after scaling (13) with the probability densities associated with the distribution $\mathbb{P}$, we obtain the following NA Lagrangian relaxation problem:

$$ \mathcal{L}_{\text{LR}}^{\text{NA}}(\mathbb{P}, \lambda_1(\cdot), \dots, \lambda_T(\cdot)) := \min z + \sum_{t \in [T]} \mathbb{E}_{\xi^T \sim \mathbb{P}} \left[ \lambda_t(\xi^T)^\top (y_t(\xi^T) - \mathbb{E}_{\xi^{T'} \sim \mathbb{P}} [y_t(\xi^{T'}) | \xi^{t'} = \xi^t]) \right] \quad (14\text{a}) $$

$$ \text{s.t.} \quad (z, y_1(\xi^T), \dots, y_T(\xi^T)) \in Y(\xi^T) \quad \xi^T \in \Xi^T. \quad (14\text{b}) $$

whose optimal objective value yields a lower bound for the original MSARO problem for any $\mathbb{P} \in \mathcal{P}^>$. The NA Lagrangian dual problem aims to find the best bound among all such Lagrangian relaxation bounds:

$$ \mathcal{L}^{\text{NA}}(\mathbb{P}) := \max_{\lambda_1(\cdot), \dots, \lambda_T(\cdot)} \mathcal{L}_{\text{LR}}^{\text{NA}}(\mathbb{P}, \lambda_1(\cdot), \dots, \lambda_T(\cdot)). \quad (15) $$

The following proposition shows that regardless of the choice of $\mathbb{P}$, $\mathcal{L}^{\text{NA}}(\mathbb{P})$ is an exact dual bound for MSARO problems with continuous recourse (proof in Appendix B).

**PROPOSITION 2.** Let $\mathbb{P}$ be any probability measure in $\mathcal{P}^>$. For MSARO problems with continuous recourse, (15) is a strong dual of (2), i.e., $\mathcal{L}^{\text{NA}}(\mathbb{P}) = \nu^*$.

We remark that our construction of the dual postulates that we multiply the constraints with a density function whose support is $\Xi^T$. Therefore, the strictly positive density property of $\mathcal{P}^>$ is necessary for the exactness of our formulation. Proposition 2 suggests that the solution of problem (15) gives an exact dual bound for (2) (hence for (12)) if all decision variables are continuous. Although this bound is not necessarily exact in the case of mixed-integer recourse, the potential of leveraging the literature of multistage stochastic programming in achieving a dual bound for MSARO is quite appealing.

The objective function of the NA Lagrangian relaxation problem (14) contains (conditional) expectations of decision variables $y_t(\xi^T)$, which is computationally challenging. Because of our initial assumption that $X_t(\xi^t)$ are bounded (Assumption 3) we have, by letting $\tilde{Y}(\xi^T) := \text{proj}_{y_1(\cdot), \dots, y_T(\cdot)} Y(\xi^T)$, that $\mathbb{E}_{\xi^T \sim \mathbb{P}}[\text{diam}(\tilde{Y}(\xi^T))] < +\infty$. Since we further have that $\mathbb{E}_{\xi^T \sim \mathbb{P}}[\lambda_t(\xi^T)] < +\infty$, we can apply Lemma 1 of (Daryalal et al. 2024) and replace the expectation term for $t$ in the objective function (14a) with:

$$ \mathbb{E}_{\xi^T \sim \mathbb{P}} \left[ \left( \lambda_t(\xi^T) - \mathbb{E}_{\xi^{T'} \sim \mathbb{P}} [\lambda_t(\xi^{T'}) | \xi^{t'} = \xi^t] \right)^\top y_t(\xi^T) \right]. $$

For given $\lambda_t(\xi^T)$, this exchange allows us to compute the coefficients of $y_t(\xi^T)$ in the Lagrangian relaxation problem. Still, the optimal form of the dual functionals $\lambda_t(\cdot)$ need to be determined, making the problem (15) computationally intractable. In the next section, we restrict these Lagrangian multipliers to follow LDRs and obtain a restricted dual problem with decision variables of smaller (finite) dimension. Furthermore, this new dual problem is amenable to well-known solution techniques from the literature of two-stage stochastic programming which are designed to approximately solve a problem with expectation in the objective function.
difficulties and in general are not desirable. To avoid these inequalities, in the following discussions we modify the DRO model to improve its numerical behavior. Denote by $\bar{\mathcal{P}}$ a superset of $\mathcal{P}$ that also admits distributions with support $\bar{\Xi} \subset \Xi$, i.e., allows $p(\xi^T) = 0$ for some $\xi^T \in \Xi$. Consider the problem:

$$ \bar{\nu}_R^{\text{NA}} = \max_{\mathbb{P} \in \bar{\mathcal{P}}} \mathcal{L}_R^{\text{NA}}(\mathbb{P}). \quad (18) $$

Because (18) is a relaxation of (17), it yields an upper bound for $\nu_R^{\text{NA}}$. Thus, it does not immediately follow that such a bound is a valid lower bound for $\nu^*$. In the following proposition, we prove that problem (18) indeed leads to a valid dual (lower) bound for the MSARO problem (1).

**PROPOSITION 3.** $\bar{\nu}_R^{\text{NA}}$ is a valid lower bound for $\nu^*$.

*Proof* We have shown that $\mathcal{L}_R^{\text{NA}}(\mathbb{P}) \le \nu^*$ for any $\mathbb{P} \in \mathcal{P}$. Now we claim that, for any $\hat{\mathbb{P}} \in \bar{\mathcal{P}} \setminus \mathcal{P} = \{\mathbb{P} \in \bar{\mathcal{P}} \mid \exists \xi^T \in \Xi : p^\mathbb{P}(\xi^T) = 0\}$ where $p^\mathbb{P}: \Xi \to \mathbb{R}_+$ is the density function of $\hat{\mathbb{P}}$, the inequality $\mathcal{L}_R^{\text{NA}}(\hat{\mathbb{P}}) \le \nu^*$ also holds. Consider the nonanticipativity constraints (10e) in the NA reformulation of the MSARO problem (1). For a fixed $\hat{\mathbb{P}}$ and $\xi^T$ in the uncertainty set, and all $\xi'^T \in \Xi$ sharing the same history with $\xi^T$ up to $t$ ($\xi^t = \xi'^t$), multiply both sides of their associated constraint (10e) by $p^\mathbb{P}(\xi'^T)$ and integrate over all such $\xi'^T$. Note that, any $p^\mathbb{P}(\xi'^T) = 0$ removes a nonanticipativity constraint from the problem. Thus, constraints (11) obtained with respect to $\hat{\mathbb{P}}$ are not equivalent to constraints (10e) and are simply an aggregation hence a valid equality. Therefore, replacing the nonanticipativity constraints with (11) obtained with respect to $\hat{\mathbb{P}}$, results in a relaxation of (1) and not an exact reformulation. Subsequently, relaxation of (11) and construction of the Lagrangian dual problem (13) with respect to $\hat{\mathbb{P}}$ yields a relaxation of a relaxation of a minimization problem, where for each fixed $\lambda_t(\cdot)$ we have $\mathcal{L}_{\text{LR}}^{\text{NA}}(\hat{\mathbb{P}}, \lambda_1(\cdot), \dots, \lambda_T(\cdot)) \le \nu^*$, and consequently $\mathcal{L}^{\text{NA}}(\hat{\mathbb{P}}) \le \nu^*$. From the primal characterizations (14) and (16), the inequality $\mathcal{L}_R^{\text{NA}}(\hat{\mathbb{P}}) \le \mathcal{L}^{\text{NA}}(\hat{\mathbb{P}})$ follows, which completes the proof. □

We refer to the DRO problem (18) as the *restricted NA dual*. It can be written explicitly as:

$$ \bar{\nu}_R^{\text{NA}} = \max_{\mathbb{P}, \alpha} Q(\mathbb{P}, \alpha) \quad (19a) $$

$$ \text{s.t.} \quad \alpha_t \in \mathbb{R}^{K_t} \quad t \in [T] \quad (19b) $$

$$ \mathbb{P} \in \bar{\mathcal{P}}, \quad (19c) $$

where

$$ Q(\mathbb{P}, \alpha) = \min z + \sum_{t \in [T]} \mathbb{E}_{\xi^T \sim \mathbb{P}} \left[ \left( \left( \Psi_t(\xi^T) - \mathbb{E}_{\xi'^T \sim \mathbb{P}} \left[ \Psi_t(\xi'^T) \mid \xi'^t = \xi^t \right] \right) \alpha_t \right)^\top y_t(\xi^T) \right] \quad (20a) $$

$$ \text{s.t.} \quad (z, y_1(\xi^T), \dots, y_T(\xi^T)) \in \text{conv}(Y(\xi^T)) \quad \xi^T \in \Xi. \quad (20b) $$

Problem (19) is a two-stage stochastic program (2SP) and can benefit from its rich literature. In the next section, we build on well-known stochastic programming techniques to design a decomposition method for the restricted NA dual (19).
## 4.2. Lagrangian Dual Decision Rules

We restrict the NA Lagrangian dual problem (15) for a given $\mathbb{P} \in \mathcal{P}^>$ by enforcing LDRs on the Lagrangian multipliers, referred to as Lagrangian dual decision rules (LDDRs). For a set of pre-determined basis functions $\Psi_t: \Xi^T \to \mathbb{R}^{n_t \times K_t}$ and LDR decision variables $\alpha_t \in \mathbb{R}^{K_t}$, we restrict the form of $\lambda_t(\xi^T)$ at stage $t \in [T]$ as follows:

$$ \lambda_t(\xi^T) = \Psi_t(\xi^T)\alpha_t, $$

which gives us a restricted NA Lagrangian dual problem with respect to $\mathbb{P}$:

$$ \mathcal{L}_R^{\text{NA}}(\mathbb{P}) := \max_{\alpha_1, \dots, \alpha_T} \mathcal{L}_{\text{LR}}^{\text{NA}}(\mathbb{P}, \Psi_1(\xi^T)\alpha_1, \dots, \Psi_T(\xi^T)\alpha_T). \quad (16) $$

Since problem (16) is a restriction of (15), we have $\mathcal{L}_R^{\text{NA}}(\mathbb{P}) \le \mathcal{L}^{\text{NA}}(\mathbb{P})$.

Using Lemma 2 in (Daryalal et al. 2024), the primal characterization of $\mathcal{L}_R^{\text{NA}}(\mathbb{P})$ is:

$$ \min z \tag{17a} $$

$$ \text{s.t.} \quad (z, y_1(\xi^T), \dots, y_T(\xi^T)) \in \text{conv}(Y(\xi^T)) \qquad \xi^T \in \Xi^T \tag{17b} $$

$$ \mathbb{E}_{\xi^T \sim \mathbb{P}} \left[ \Psi_t(\xi^T)^\top (y_t(\xi^T) - \mathbb{E}_{\xi'^T \sim \mathbb{P}} [y_t(\xi'^T) \mid \xi'^t = \xi^t]) \right] = 0 \qquad t \in [T] \tag{17c} $$

For a given $\mathbb{P}$, comparing (17) to the primal characterization of the (unrestricted) NA Lagrangian dual problem (15) (provided in (42) of Appendix B), it is clear that, the former is a relaxation of the latter since constraints (17c) are an aggregation of their counterpart (42c). Consequently, unlike $\mathcal{L}^{\text{NA}}(\mathbb{P})$, even for MSARO with continuous recourse, $\mathcal{L}_R^{\text{NA}}(\mathbb{P})$ is not necessarily a strong bound. Furthermore, due to constraints (17c), the strength of the restricted NA Lagrangian dual bound depends on the choice of probability measure $\mathbb{P}$. A similar observation was made by Kuhn et al. (2011) and Hadjiyiannis et al. (2011) concerning their dual bound for MSARO problems with only continuous variables. This observation motivates us to optimize over the probability distribution $\mathbb{P}$ and LDR variables $\alpha$ to find the best such dual bound. As such, we propose to solve a distribution optimization (DO) problem over the set of probability distributions $\mathcal{P}^>$, defined as follows:

$$ \nu_R^{\text{NA-DO}} := \sup_{\mathbb{P} \in \mathcal{P}^>} \mathcal{L}_R^{\text{NA}}(\mathbb{P}) \tag{18} $$

where, as before, $\mathcal{P}^> = \{\mathbb{P} \mid p^>(\xi^T) > 0, \xi^T \in \Xi^T\}$.

In linear and mixed-integer programming, strict inequalities such as the ones required for $\mathbb{P} \in \mathcal{P}^>$ (the strictly positive density property for all $\xi^T \in \Xi^T$) often cause numerical and theoretical difficulties, thus are not desirable. To avoid these inequalities, in the following discussions, we
### 3.3. Solving the Restricted NA Dual Problem

For a discrete set of scenarios and given $\alpha$, the expectation terms in the objective function of $Q(\mathbb{P}, \alpha)$ in (19), i.e., the coefficients of the $y$ variables, can be calculated. In the following sections we use this fact to solve or approximate the restricted NA dual problem. We first present a solution approach designed for finite discrete uncertainty sets. Then we show how we can use this method for a general uncertainty set.

#### 3.3.1. Discrete Uncertainty Set
In a discrete uncertainty set with realizations $\xi^T \in \Xi$, $\bar{\mathcal{P}}$ can be modeled by a set of (in)equalities, that is $\bar{\mathcal{P}} = \{\rho \in \mathbb{R}_+^{|\Xi|} | \mathbf{1}^\top \rho = 1\}$, where a vector $\rho^\mathbb{P} \in \bar{\mathcal{P}}$ characterizes the probability measure $\mathbb{P}$ such that for all $\xi^T \in \Xi$, $\rho_{\xi^T}^\mathbb{P}$ is the probability of scenario $\xi^T$ with respect to $\mathbb{P}$. Given that $\Xi$ is discrete, the expectation term in the objective function (20a) can be written as:

$$ \sum_{t \in [T]} \sum_{\xi^T \in \Xi} \rho_{\xi^T}^{\mathbb{P}} \left( \left( \Psi_t(\xi^T) - \sum_{\substack{\xi'^T \in \Xi: \\ \xi'^t = \xi^t}} \rho_{\xi'^T}^{\mathbb{P}} \Psi_t(\xi'^T) \right) \alpha_t \right)^{\top} y_t(\xi^T). $$

Let $\Psi_{tk}(\xi^T)$ be the $k^{th}$ column of the matrix $\Psi_t(\xi^T)$. Then $Q(\mathbb{P}, \alpha)$ can be expressed as:

$$ \begin{array}{ll} \min & z + \sum_{t \in [T]} \sum_{\xi^T \in \Xi} \rho_{\xi^T}^{\mathbb{P}} \left( \sum_{k \in [K_t]} \left( \Psi_{tk}(\xi^T) - \sum_{\substack{\xi'^T \in \Xi: \\ \xi'^t = \xi^t}} \rho_{\xi'^T}^{\mathbb{P}} \Psi_{tk}(\xi'^T) \right) \alpha_{tk} \right)^{\top} y_t(\xi^T) \\ \text{s.t.} & (z, y_1(\xi^T), \dots, y_T(\xi^T)) \in Y(\xi^T) \quad \xi^T \in \Xi. \end{array} $$

Let $\gamma_{tk\xi^T} = \rho_{\xi^T}^\mathbb{P} \alpha_{tk}$ and $\beta_{t\xi^T} = \sum_{k \in [K_t]} (\gamma_{tk\xi^T} \Psi_{tk}(\xi^T) - \sum_{\substack{\xi'^T \in \Xi: \\ \xi'^t = \xi^t}} \rho_{\xi'^T}^\mathbb{P} \gamma_{tk\xi^T} \Psi_{tk}(\xi'^T))$. Then with a change of variables in the restricted NA dual (19) we have:

$$ \bar{\nu}_R^{\text{NA}} = \max_{\rho, \alpha, \gamma, \beta} Q(\beta) \qquad (21a) $$

$$ \text{s.t.} \quad \sum_{\xi^T \in \Xi} \rho_{\xi^T} = 1 \qquad (21b) $$

$$ \gamma_{tk\xi^T} = \rho_{\xi^T} \alpha_{tk} \qquad t \in [T], k \in [K_t], \xi^T \in \Xi \qquad (21c) $$

$$ \beta_{t\xi^T} = \sum_{k \in [K_t]} \left( \gamma_{tk\xi^T} \Psi_{tk}(\xi^T) - \sum_{\substack{\xi'^T \in \Xi: \\ \xi'^t = \xi^t}} \rho_{\xi'^T} \gamma_{tk\xi^T} \Psi_{tk}(\xi'^T) \right) \quad t \in [T], \xi^T \in \Xi \qquad (21d) $$

$$ \rho \in \mathbb{R}_+^{|\Xi|} \qquad (21e) $$

$$ \alpha_t \in \mathbb{R}^{K_t}, \quad \gamma \in \mathbb{R}^{T \times K_t \times |\Xi|}, \quad \beta \in \mathbb{R}^{n_t \times T \times |\Xi|} \qquad t \in [T] \qquad (21f) $$

where

$$ Q(\beta) = \min \left\{ z + \sum_{t \in [T]} \sum_{\xi^T \in \Xi} \beta_{t\xi^T}^{\top} y_t(\xi^T) \mid (z, y_1(\xi^T), \dots, y_T(\xi^T)) \in Y(\xi^T), \xi^T \in \Xi \right\}. \qquad (22) $$
modify the DO model to improve its numerical behaviour. Denote by $\mathcal{P}^\ge$ a superset of $\mathcal{P}^>$ that also admits distributions that allow $p(\xi^T) = 0$ for some $\xi^T \in \Xi^T$. Consider the problem:

$$ \bar{\nu}_R^{\text{NA-DO}} := \max_{\mathbb{P} \in \mathcal{P}^\ge} \mathcal{L}_R^{\text{NA}}(\mathbb{P}). \quad (19) $$

Because (19) is a relaxation of (18), it yields an upper bound for $\nu_R^{\text{NA-DO}}$. Thus, it does not immediately follow that such a bound is a valid lower bound for the optimal value of the MSARO problem, $\nu^*$. The following proposition shows that (19) indeed leads to a valid dual (lower) bound (see Appendix B for proof).

**PROPOSITION 3.** $\bar{\nu}_R^{\text{NA-DO}}$ is a lower bound for $\nu^*$.

Hereafter, we refer to (19) or its equivalent explicit form

$$ \bar{\nu}_R^{\text{NA-DO}} = \max_{\mathbb{P}, \alpha} Q(\mathbb{P}, \alpha) \quad (20a) $$

$$ \text{s.t.} \quad \alpha_t \in \mathbb{R}^{K_t}, \quad t \in [T] \quad (20b) $$

$$ \mathbb{P} \in \mathcal{P}^\ge, \quad (20c) $$

as the DO problem, where

$$ Q(\mathbb{P}, \alpha) := \min z + \sum_{t \in [T]} \mathbb{E}_{\xi^T \sim \mathbb{P}} \left[ \left( \left( \Psi_t(\xi^T) - \mathbb{E}_{\xi^T \sim \mathbb{P}} [\Psi_t(\xi^T^t) | \xi^t = \xi^t] \right) \alpha_t \right)^{\top} y_t(\xi^T) \right] \quad (21a) $$

$$ \text{s.t.} \quad (z, y_1(\xi^T), \dots, y_T(\xi^T)) \in Y(\xi^T) \quad \xi^T \in \Xi^T. \quad (21b) $$

Problem (21) is a two-stage stochastic program (2SP) and can benefit from its rich literature. In the next section, we build on well-known stochastic programming techniques to design a decomposition method to solve the DO problem (20).

### 4.3. Solving the DO Problem

There are two main challenges associated with the solution of the DO problem: the expectation terms in the objective of (21a) and the max-min structure in (20).

If the uncertainty set of the MSARO problem, $\Xi^T$, is not discrete, objective function (21a) includes the expectation of a nonsmooth concave function. Further, even when $\Xi^T$ is discrete calculating the expectation term exactly can be prohibitive from a computational point of view. The literature of two-stage stochastic programming addresses such a difficulty by means of sampling-based approaches that replace the expectation in the objective function with the average of a sample drawn from the underlying distribution and has favourable theoretical convergence results (see e.g., Shapiro et al. (2009)). We follow the sample average approximation (SAA) approach in
Model (21) can be solved via a cutting-plane method in which $Q(\beta)$ is approximated by a set of linear inequalities. At each iteration, we solve the following bilinear program as the master problem:

$$ \max_{p,\alpha,\gamma,\beta} \left\{ \eta \mid (21b) - (21f), (\eta, \beta) \in \mathcal{H} \right\}, \qquad (23) $$

where $\eta$ is an auxiliary variable representing $Q(\beta)$, and $\mathcal{H}$ is a set described by optimality cuts approximating $Q(\beta)$. Note that, as $\beta$ only parameterizes the objective function of $Q(\beta)$, i.e., it does not impact the feasibility space, there is no need for feasibility cuts. With $(\hat{\eta}, \hat{\beta})$ returned from solving the master problem (23), we solve the subproblem (22) to compute $Q(\hat{\beta})$, resulting in $\hat{y}_t(\xi^T)$ as the optimal solution. If $\hat{\eta} \le Q(\hat{\beta})$, we have found the optimal solution of the restricted NA dual. Otherwise we add the following optimality cut to the master problem:

$$ \eta \le Q(\hat{\beta}) + \sum_{t \in [T]} \sum_{\xi^T \in \Xi} (\beta_{t\xi^T} - \hat{\beta}_{t\xi^T})^\top \hat{y}_t(\xi^T). \qquad (24) $$

This procedure continues until no more optimality cuts are found. The objective function of the subproblem, $Q(\beta)$, is a concave function in $\beta$ (pointwise minimum of linear functions with respect to $\beta$). Each cut (24) is a hyperplane approximating the subproblem from above. The cutting-plane method iteratively finds improving approximations of $Q(\beta)$. Enforcing large lower and upper bounds on LDR variables $\alpha$ can ensure that the algorithm converges to $\epsilon$-optimality.

## Discrete Uncertainty in Continuous MSARO

As a special case, we study continuous MSARO with discrete uncertainty. Our goal here is to use this particular structure and derive a monolithic formulation as an alternative to the cutting-plane algorithm, to leverage off-the-shelf solvers. Denote by $u_{\xi^T}$, $v_{t\xi^T}$ and $w_{t\xi^T}$, the dual variables associated with the set of constraints described by $Y(\xi^T)$ for given $\xi^T \in \Xi$ (corresponding to (10b), (10c) and (10d), respectively). The LP dual of the inner minimization problem (22), i.e., the subproblem of the cutting-plane algorithm, is:

$$ Q^D(\beta) = \max \sum_{t \in [T]} \sum_{\xi^T \in \Xi} b_t(\xi^T)^\top v_{t\xi^T} + \sum_{t \in [T]} \sum_{\xi^T \in \Xi} d_t(\xi^T)^\top w_{t\xi^T} \qquad (25a) $$

$$ \text{s.t.} \quad \sum_{\xi^T \in \Xi} u_{\xi^T} = 1 \qquad (25b) $$

$$ -c_T(\xi^T)u_{\xi^T} + A_T(\xi^T)^\top v_{t\xi^T} + D_{T\xi^T}^\top w_{t\xi^T} = \beta_{T\xi^T} \qquad \xi^T \in \Xi \qquad (25c) $$

$$ -c_t(\xi^t)u_{\xi^T} + A_t(\xi^t)^\top v_{t\xi^T} + D_t(\xi^t)^\top w_{t\xi^T} + B_{t+1}(\xi^{t+1})^\top v_{t+1,\xi^T} = \beta_{t\xi^T} \qquad t \in [T-1], \xi^T \in \Xi \qquad (25d) $$

$$ u_{\xi^T} \ge 0 \qquad \xi^T \in \Xi \qquad (25e) $$
overcoming the first challenge and show in Section 4.3.1 that it leads to a valid dual bound for the MSARO problem.

With regards to the second challenge, we propose a cutting plane algorithm that iteratively constructs improving approximations of $Q(\mathbb{P}, \alpha)$ through its supporting hyperplanes obtained by solving the SAA approximation of (21). We present this general algorithm in Section 4.3.2 and propose an alternative monolithic formulation in the special case of MSARO with continuous recourse in Section 4.3.3. Figure 4 summarizes the methods presented in this section for solving the DO model.

<figure><img src="image_6.png" /><figcaption>**Figure 4** Solution methods for the DO problem</figcaption></figure>

**4.3.1. Sample average approximation (SAA)** Let $\Omega \subseteq \Xi^T$ be a finite subset of the uncertainty set $\Xi^T$. We define the set of probability measures $\mathcal{P}_{\Omega}^{\geq}$ such that $\mathbb{P}_{\Omega} \in \mathcal{P}_{\Omega}^{\geq}$ implies that $\mathbb{P}_{\Omega}(\Omega) = 1$ and the associated density function has value zero for any realization not in $\Omega$:

$$p_{\mathbb{P}_{\Omega}}(\xi^T) = 0, \xi^T \in \Xi^T \setminus \Omega, \text{ and } p_{\mathbb{P}_{\Omega}}(\xi^T) \geq 0, \xi^T \in \Omega.$$

Since $\mathcal{P}_{\Omega}^{\geq} \subseteq \mathcal{P}^{\geq}$ and $\bar{\nu}_{R}^{\text{NA-DO}}$ is obtained by maximizing $\mathcal{L}_{R}^{\text{NA}}(\mathbb{P})$ over $\mathbb{P} \in \mathcal{P}^{\geq}$, we have, for any $\mathbb{P}_{\Omega} \in \mathcal{P}_{\Omega}^{\geq}$, that $\mathcal{L}_{R}^{\text{NA}}(\mathbb{P}_{\Omega}) \leq \bar{\nu}_{R}^{\text{NA-DO}} \leq \nu^*$. It then follows that

$$\max_{\mathbb{P}_{\Omega} \in \mathcal{P}_{\Omega}} \mathcal{L}_{R}^{\text{NA}}(\mathbb{P}_{\Omega}) \leq \bar{\nu}_{R}^{\text{NA-DO}} \leq \nu^*. \qquad (22)$$

We finally have that

$$\max_{\mathbb{P}_{\Omega} \in \mathcal{P}_{\Omega}^{\geq}} \mathcal{L}_{R}^{\text{NA}}(\mathbb{P}_{\Omega}) = \max_{\mathbb{P}_{\Omega}, \alpha} Q(\mathbb{P}_{\Omega}, \alpha) \qquad (23)$$

$$\text{s.t.} \quad \alpha_t \in \mathbb{R}^{K_t} \quad t \in [T]$$

$$\mathbb{P}_{\Omega} \in \mathcal{P}_{\Omega}^{\geq},$$
$$v_{t\xi^T}, w_{t\xi^T} \le \mathbf{0} \qquad t \in [T], \xi^T \in \Xi. \qquad (25f)$$

Merging the two maximization problems in (21), we get the following monolithic bilinear program:

$$\bar{\nu}_R^{\text{NA}} = \max \sum_{t \in [T]} b_t (\xi^T)^\top v_{t\xi^T} + \sum_{t \in [T]} \sum_{\xi^T \in \Xi} d_t (\xi^t)^\top w_{t\xi^T} \qquad (26a)$$

$$\text{s.t.} \quad (21b) - (21f) \qquad (26b)$$

$$(25b) - (25f). \qquad (26c)$$

Recall that, in adapting the C&CG algorithm for two-stage LDRs and Remark 1, we use KKT conditions to obtain a monolithic model of the subproblem and the bilinear terms in the constraints can be linearized by adding big-$M$ constraints. The same technique would not lead to a linear reformulation of (21). Indeed, the source of nonlinearity in constraint (8e) and (8f) is the multiplication of a linear function and a dual variable. Thus, we are able to deal with them separately (and make sure that at least one of the two terms is evaluated to zero) and derive a linear model. However, in (21) the nonlinearity is rooted in the constraints of the outer maximization problem itself and enforcing the complementary slackness through separate constraints does not resolve it. Nevertheless, there is a large body of research on solution methods for bilinear problems that can be used in solving model (26). Furthermore, many nonlinear optimization solvers such as MOSEK (ApS 2022) include algorithms specifically tailored for large-scale quadratic problems and even a commercial MIP solver, Gurobi, has recently added a bilinear feature that solves non-convex quadratic optimization models to global optimality (Gurobi Optimization, LLC 2022).

### 3.3.2. Continuous Uncertainty Set

If the uncertainty set of MSARO, $\Xi$, is not discrete, objective function (20a) includes the expectation of a nonsmooth concave function, making (19) unsuitable for many existing solution methods for 2SPs. The literature of two-stage stochastic programming addresses this difficulty by means of sampling-based approaches that approximate the underlying problem. Sample average approximation (SAA) is a well-known method that replaces the expectations in the objective function of a 2SP with the average of a given sample and has favorable theoretical convergence results (see e.g., Shapiro et al. (2009)). In this section we show that (approximately) solving the 2SP (19) over a sample returns a valid dual bound for the original MSARO problem (1).

From the discussions in Section 3, for any $\mathbb{P} \in \bar{\mathcal{P}}$, we have:

$$\mathcal{L}_R^{\text{NA}}(\mathbb{P}) \le \bar{\nu}_R^{\text{NA}} \le \nu^*. \qquad (27)$$
where for any $\mathbb{P}_{\Omega} \in \mathcal{P}_{\Omega}^{\ge}$ the expectation terms in the objective function (21a), used in calculating $\mathcal{Q}(\mathbb{P}_{\Omega}, \alpha)$, are replaced by their sample average over $\Omega$.

In the remainder of this section, we omit the notation $\Omega$ and write our models with a finite discrete uncertainty set $\Xi^T$ which can either be the full uncertainty set of the MSARO problem or a set of realizations sampled from it.

### 4.3.2. MSARO with mixed-integer recourse
In a discrete uncertainty set with realizations $\xi^T \in \Xi^T$, $\mathcal{P}^\ge$ can be modeled by a set of (in)equalities, that is $\mathcal{P}^\ge = \{\rho \in \mathbb{R}_+^{\Xi^T} | 1^\top \rho = 1\}$, where a vector $\rho^\mathbb{P} \in \mathcal{P}^\ge$ characterizes the probability measure $\mathbb{P}$ such that for all $\xi^T \in \Xi^T$, $\rho^\mathbb{P}_{\xi^T}$ is the probability of realization $\xi^T$ with respect to $\mathbb{P}$. Define $\rho^\mathbb{P}_{|\xi^T} := \sum_{\xi^{T'} \in \Xi^T: \xi^{T'}=\xi^T} \rho^\mathbb{P}_{\xi^{T'}|\xi^T}$ as the sum of the probabilities of all realizations sharing the same history $\xi^t$ up to stage $t$ and let $\rho^\mathbb{P}_{\xi^{T'}|\xi^t}$ be the conditional probability of realization $\xi^{T'} \in \Xi^T$ given history $\xi^t$. More precisely, given $\xi^t$, for $\xi^{T'} \in \Xi^T$ with $\xi^{t'} = \xi^t$ we have that $\rho^\mathbb{P}_{\xi^{T'}|\xi^t} = \frac{\rho^\mathbb{P}_{\xi^{T'}|\xi^t}}{\rho^\mathbb{P}_{\xi^t}}$ if $\rho^\mathbb{P}_{|\xi^t} > 0$ and that $\rho^\mathbb{P}_{\xi^{T'}|\xi^t} = 0$ otherwise. Then, the expectation term in the objective function (21a) can be written as:

$$ \sum_{t \in [T]} \sum_{\xi^T \in \Xi^T} \rho^\mathbb{P}_{\xi^T} \left( \left( \Psi_t(\xi^T) - \sum_{\substack{\xi^{T'} \in \Xi^T: \\ \xi^{t'} = \xi^t}} \rho^\mathbb{P}_{\xi^{T'}|\xi^t} \Psi_t(\xi^{T'}) \right) \alpha_t \right)^{\top} y_t(\xi^T). $$

Let $\Psi_{tk}(\xi^T)$ be the $k^{th}$ column of the matrix $\Psi_t(\xi^T)$. Then $\mathcal{Q}(\mathbb{P}, \alpha)$ can be expressed as:

$$ \begin{array}{ll} \min & z + \sum_{t \in [T]} \sum_{\xi^T \in \Xi^T} \rho^\mathbb{P}_{\xi^T} \left( \sum_{k \in [K_t]} \left( \Psi_{tk}(\xi^T) - \sum_{\substack{\xi^{T'} \in \Xi^T: \\ \xi^{t'} = \xi^t}} \rho^\mathbb{P}_{\xi^{T'}|\xi^t} \Psi_{tk}(\xi^{T'}) \right) \alpha_{tk} \right)^{\top} y_t(\xi^T) \\ \text{s.t.} & (z, y_1(\xi^T), \dots, y_T(\xi^T)) \in Y(\xi^T) \quad \xi^T \in \Xi^T. \end{array} $$

Let $\gamma_{tk\xi^T} := \rho^\mathbb{P}_{\xi^T} \alpha_{tk}$ and define $\beta_{t\xi^T} := \sum_{k \in [K_t]} \left( \gamma_{tk\xi^T} \Psi_{tk}(\xi^T) - \sum_{\substack{\xi^{T'} \in \Xi^T: \\ \xi^{t'} = \xi^t}} \rho^\mathbb{P}_{\xi^{T'}|\xi^t} \gamma_{tk\xi^T} \Psi_{tk}(\xi^{T'}) \right)$ as the coefficient vector of variables $y_t(\xi^T)$. Then with a change of variables in the DO problem (20) we have:

$$ \bar{\nu}_R^{\text{NA-DO}} = \max_{\rho, \alpha, \gamma, \beta} \mathcal{Q}(\beta) \qquad (24a) $$

$$ \text{s.t.} \quad \sum_{\xi^T \in \Xi^T} \rho_{\xi^T} = 1 \qquad (24b) $$

$$ \gamma_{tk\xi^T} = \rho_{\xi^T} \alpha_{tk} \qquad t \in [T], k \in [K_t], \xi^T \in \Xi^T \qquad (24c) $$

$$ \sum_{\substack{\xi^{T'} \in \Xi^T: \\ \xi^{t'} = \xi^t}} \rho^\mathbb{P}_{\xi^{T'}|\xi^T} \beta_{t\xi^T} - \sum_{k \in [K_t]} \left( \sum_{\substack{\xi^{T'} \in \Xi^T: \\ \xi^{t'} = \xi^t}} \rho^\mathbb{P}_{\xi^{T'}|\xi^T} \gamma_{tk\xi^T} \Psi_{tk}(\xi^T) - \sum_{\substack{\xi^{T'} \in \Xi^T: \\ \xi^{t'} = \xi^t}} \rho_{\xi^{T'}|\xi^T} \gamma_{tk\xi^T} \Psi_{tk}(\xi^{T'}) \right) = 0 \qquad (24d) $$

$$ -M \rho_{\xi^T} \le \beta_{t\xi^T} \le M \rho_{\xi^T} \qquad t \in [T], \xi^T \in \Xi^T \qquad (24e) $$
First inequality † is valid because the bound $\bar{\nu}_R^{NA}$ is the result of maximization over all $\mathbb{P} \in \bar{\mathcal{P}}$, and the second inequality ‡ is from Proposition 3. For a given finite sample $\Omega \subseteq \Xi$, we define a probability measure $\mathbb{P}_\Omega$ with $\mathbb{P}_\Omega(\Omega) = 1$ such that the density function assigns zero to any scenario not in $\Omega$:

$$p^{\mathbb{P}_{\Omega}}(\boldsymbol{\xi}^{T}) = 0, \quad \boldsymbol{\xi}^{T} \in \Xi \setminus \Omega, \text{ and } p^{\mathbb{P}_{\Omega}}(\boldsymbol{\xi}^{T}) \geq 0, \quad \boldsymbol{\xi}^{T} \in \Omega.$$

Clearly $\mathbb{P}_\Omega \in \bar{\mathcal{P}}$, thus inequalities (27) hold for $\mathbb{P}_\Omega$, i.e., $\mathcal{L}_R^{NA}(\mathbb{P}_\Omega) \le \bar{\nu}_R^{NA} \le \nu^*$, and $\mathcal{L}_R^{NA}(\mathbb{P}_\Omega)$ is a valid dual bound on $\nu^*$. As $\Omega$ is a finite set, we can represent $\mathbb{P}_\Omega$ with $\rho^{\mathbb{P}_\Omega} \in \mathbb{R}_+^{|\Omega|}$ such that $\mathbf{1}^\top \rho^{\mathbb{P}_\Omega} = 1$. Denote by $\mathcal{P}^\Omega$, the set of all such vectors $\rho$. Consequently, solving the DRO problem (21) with discrete uncertainty set $\Omega$ finds the best bound over all probability measures $\mathbb{P} \in \mathcal{P}^\Omega$, which is a valid dual bound on $\nu^*$. Calling the optimal value of this new DRO problem $\bar{\nu}_R^{NA}(\Omega)$, the following proposition establishes the validity of the bound $\bar{\nu}_R^{NA}(\Omega)$.

**PROPOSITION 4.** For any given sample $\Omega \subseteq \Xi$, weak duality holds between $\bar{\nu}_R^{NA}(\Omega)$ and $\nu^*$, i.e.,

$$\bar{\nu}_R^{NA}(\Omega) \le \bar{\nu}_R^{NA} \le \nu^*.$$

So, for a continuous uncertainty set, we can get a discrete sample $\Omega$ which turns the problem into a case of discrete uncertainty. As such, we can use the arguments in Section 3.3.1 to obtain $\bar{\nu}_R^{NA}(\Omega)$.

Figure 5 provides a summary of the methods presented in this section for solving the restricted NA dual problem. For continuous MSARO with a discrete uncertainty set, the restricted NA dual bound can be computed by solving a bilinear program, while for mixed-integer MSARO, the devised cutting-plane method is applicable. For MSARO with a continuous uncertainty set, we use SAA, i.e., we first take a sample and then approximately solve the restricted NA dual problem using the methods developed in the discrete setting.

<figure><img src="image_6.png" /><figcaption>**Figure 5** Solution methods for the restricted NA dual problem</figcaption></figure>
$$ \rho \in \mathbb{R}_+^{\Xi^T} \qquad (24f) $$

$$ \alpha_t \in \mathbb{R}^{K_t}, \quad \gamma_t \in \mathbb{R}^{K_t \times |\Xi^T|}, \quad \beta_t \in \mathbb{R}^{|\Xi^T| \times n_t} \qquad t \in [T] \qquad (24g) $$

where

$$ Q(\beta) = \min \left\{ z + \sum_{t \in [T]} \sum_{\xi^T \in \Xi^T} \beta_{t\xi^T}^{\top} y_t(\xi^T) \mid (z, y_1(\xi^T), \dots, y_T(\xi^T)) \in Y(\xi^T), \xi^T \in \Xi^T \right\}, \quad (25) $$

and $M \ge 0$ is a vector of sufficiently large numbers.

In model (24), constraints (24b) along with the nonnegativity restrictions imposed on variables $\rho$ induce a probability distribution over the realizations $\Xi^T$ with some realizations potentially assigned zero probability. For given $\xi^T$, if its associated probability $\rho_{\xi^T}$ is zero, then constraints (24c) and (24e) imply, respectively, that all variables $\gamma$ and $\beta$ indexed by $\xi^T$ are zero, as such constraints (24d) trivially hold. Otherwise, since probability $\rho_{\xi^T} > 0$, we have that $\sum_{\substack{\xi'^T \in \Xi^T \\ \xi'^t = \xi^t}} \rho_{\xi'^T}^{\mathbb{P}} > 0$, as such constraints (24d) correctly impose the definition of $\beta$ by dividing all terms by $\sum_{\substack{\xi'^T \in \Xi^T \\ \xi'^t = \xi^t}} \rho_{\xi'^T}^{\mathbb{P}}$ and plugging in the definition of $\gamma$ from constraints (24c). Then given coefficients $\beta$, (25) evaluates the expected value of the optimal decisions $y(\cdot)$.

**REMARK 5.** The choice of values for the vector $M$ impacts the quality of the bound obtained from model (24). In particular, in the limit case where $M = 0$, one obtains the perfect information bound, that is, all nonanticipativity constraints are relaxed from the NA reformulation of the MSARO problem. Otherwise, for larger values of $M$, the optimal value of (24) is lower bounded by the perfect information bound since choosing $\beta = 0$ is always feasible. One can therefore expect to obtain a better bound from (24).

Model (24) can be solved via a cutting-plane method in which $Q(\beta)$ is approximated by a set of linear inequalities. At each iteration, we solve the following bilinear program as the master problem:

$$ \max_{\rho, \alpha, \gamma, \beta} \left\{ \eta \mid (24b) - (24g), (\eta, \beta) \in \mathcal{H} \right\}, \qquad (26) $$

where $\eta$ is an auxiliary variable representing $Q(\beta)$, and $\mathcal{H}$ is a set described by optimality cuts approximating $Q(\beta)$. Note that, as $\beta$ only parameterizes the objective function of $Q(\beta)$, i.e., it does not impact the feasibility space, there is no need for feasibility cuts. With $(\hat{\eta}, \hat{\beta})$ returned from solving the master problem (26), we solve the subproblem (25) to compute $Q(\hat{\beta})$, resulting in $\hat{y}_t(\xi^T)$ as the optimal solution. If $\hat{\eta} \le Q(\hat{\beta})$, we have found the optimal solution of the DO problem. Otherwise we add the following optimality cut to the master problem:

$$ \eta \le Q(\hat{\beta}) + \sum_{t \in [T]} \sum_{\xi^T \in \Xi^T} (\beta_{t\xi^T} - \hat{\beta}_{t\xi^T})^{\top} \hat{y}_t(\xi^T). \qquad (27) $$
### 3.4. Restricted Decomposable NA Dual Bound

In solving problem (21) using the cutting-plane method, we frequently optimize (22) to compute $Q(\beta)$. This can become computationally demanding when there is a large number of scenarios (i.e., large size of a discrete uncertainty set or a sample). Furthermore, we cannot decompose $Q(\beta)$ by the scenarios, as they are linked through the $z$ variables. In this section, we present an alternative NA dual problem that can yield a potentially weaker bound than the one provided by (19), but is decomposable by scenarios, thus offering a computational advantage. With decision variables $y_t(\xi^T)$ defined as before, we introduce copy variables $z(\xi^T)$ and explicitly enforce them all to be equal. For an assigned probability distribution $\mathbb{P}$ which has a strictly positive density function for every realization of the uncertainty set, the alternative NA reformulation for the MSARO problem (1) is:

$$ \min \quad \mathbb{E}_{\xi^T \sim \mathbb{P}} [z(\xi^T)] \qquad (28a) $$

$$ \text{s.t.} \quad \sum_{t \in [T]} c_t (\xi^t)^\top y_t (\xi^T) \le z(\xi^T) \qquad \xi^T \in \Xi \qquad (28b) $$

$$ z(\xi^T) = \mathbb{E}_{\xi^T \sim \mathbb{P}} [z(\xi^T)] \qquad \xi^T \in \Xi \qquad (28c) $$

$$ (10c), (10d), (10f), (11). \qquad (28d) $$

Together, constraints (28b) and (28c) capture the semantics of the worst-case outcome, which is minimized in the objective function (28a). Relaxation of the nonanticipativity constraints (11) and (28c) with the Lagrangian multipliers $\lambda_t^y(\cdot)$ and $\lambda^z(\cdot)$, leads to the decomposable NA Lagrangian dual problem $\mathcal{L}^{\text{DNA}}(\mathbb{P}) = \max_{\mathbb{P}, \lambda^y, \lambda^z} \mathcal{L}_{LR}^{\text{DNA}}(\mathbb{P}, \lambda_1^y, \dots, \lambda_T^y, \lambda^z)$, where

$$ \mathcal{L}_{LR}^{\text{DNA}}(\mathbb{P}, \lambda_1^y, \dots, \lambda_T^y, \lambda^z) = \min \quad \mathbb{E}_{\xi^T \sim \mathbb{P}} \left[ \left(1 + \lambda^z(\xi^T) - \mathbb{E}_{\xi^T \sim \mathbb{P}} \left[ \lambda^z(\xi^T) \right] \right) z(\xi^T) \right] + \\ \sum_{t \in [T]} \mathbb{E}_{\xi^T \sim \mathbb{P}} \left[ \left( \lambda_t^y(\xi^T) - \mathbb{E}_{\xi^T \sim \mathbb{P}} \left[ \lambda_t^y(\xi^T) \mid \xi^{t'} = \xi^t \right] \right)^\top y_t(\xi^T) \right] \\ \text{s.t.} \quad (z(\xi^T), y_1(\xi^T), \dots, y_T(\xi^T)) \in Y(\xi^T) \qquad \xi^T \in \Xi, $$

with $Y(\xi^T)$ the scenario feasibility space described by constraints (28b), (10c), (10d), and (10f). We remark that, in the first expectation term in the objective function the inner expectation is not conditional, since constraints (28c) (and accordingly their associated dual functions) are not defined at every stage $t \in [T]$ as decision variable $z(\xi^T) \in \mathbb{R}$ captures the cost of an entire scenario at $T$. After substituting the decision rules $\lambda_t^y(\xi^T) = \Psi_t(\xi^T)\alpha_t^y, t \in [T]$ and $\lambda^z(\xi^T) = \Psi_T(\xi^T)\alpha^z$ in $\mathcal{L}^{\text{DNA}}(\mathbb{P})$, where $\alpha^z \in \mathbb{R}^{K_T}$, the restricted decomposable NA dual problem is:

$$ \nu_R^{\text{DNA}} = \max_{\mathbb{P}, \alpha^y, \alpha^z} \mathcal{L}_{LR}^{\text{DNA}}(\mathbb{P}, \Psi_1(\cdot)\alpha_1^y, \dots, \Psi_T(\cdot)\alpha_T^y, \Psi_T(\cdot)\alpha^z). $$
This procedure continues until no more optimality cuts are found. The objective function of the subproblem, $Q(\beta)$, is a concave function in $\beta$ (pointwise minimum of linear functions with respect to $\beta$). Each cut (27) is a hyperplane approximating the subproblem from above. The cutting-plane method iteratively finds improving approximations of $Q(\beta)$. We note that it is not necessary to execute the cutting-plane procedure until convergence in order to obtain a valid dual bound. Indeed at each iteration of the algorithm the value $Q(\hat{\beta})$ provides a valid dual bound for $\nu^*$.

REMARK 6. In our implementation for the numerical results presented in Section 5, rather than creating the copies of first-stage variables $y_1(\cdot)$ and relaxing their nonanticipativity constraints, we keep them as static variables in (25), same as variable $z$.

### 4.3.3. MSARO with continuous recourse
As a special case, we study MSARO with continuous recourse. Our goal here is to use this particular structure and derive a monolithic formulation as an alternative to the cutting-plane algorithm, to leverage off-the-shelf solvers. Denote by $u_{\xi^T}$, $v_{t\xi^T}$ and $w_{t\xi^T}$, the dual variables associated with the set of constraints described by $Y(\xi^T)$ for given $\xi^T \in \Xi^T$ (corresponding to (12b), (12c) and (12d), respectively). The linear programming dual of the inner minimization problem (25), i.e., the subproblem of the cutting-plane algorithm, is:

$$Q^D(\beta) = \max \sum_{t \in [T]} \sum_{\xi^T \in \Xi^T} b_t(\xi^T)^\top v_{t\xi^T} + \sum_{t \in [T]} \sum_{\xi^T \in \Xi^T} d_t(\xi^T)^\top w_{t\xi^T} \quad (28a)$$

$$\text{s.t.} \quad \sum_{\xi^T \in \Xi^T} u_{\xi^T} = 1 \quad (28b)$$

$$-c_T(\xi^T)u_{\xi^T} + A_T(\xi^T)^\top v_{T\xi^T} + D_{T\xi^T}^\top w_{T\xi^T} - \beta_{T\xi^T} = \mathbf{0} \quad \xi^T \in \Xi^T \quad (28c)$$

$$-c_t(\xi^t)u_{\xi^T} + A_t(\xi^t)^\top v_{t\xi^T} + D_t(\xi^t)^\top w_{t\xi^T} +$$

$$B_{t+1}(\xi^{t+1})^\top v_{t+1,\xi^T} - \beta_{t\xi^T} = \mathbf{0} \quad t \in [T-1], \xi^T \in \Xi^T \quad (28d)$$

$$u_{\xi^T} \ge 0 \quad \xi^T \in \Xi^T \quad (28e)$$

$$v_{t\xi^T}, w_{t\xi^T} \le \mathbf{0} \quad t \in [T], \xi^T \in \Xi^T. \quad (28f)$$

Merging the two maximization problems in (24), we get the monolithic bilinear program:

$$\bar{\nu}_R^{\text{NA-DO}} = \max \sum_{t \in [T]} b_t(\xi^T)^\top v_{t\xi^T} + \sum_{t \in [T]} \sum_{\xi^T \in \Xi^T} d_t(\xi^T)^\top w_{t\xi^T} \quad (29a)$$

$$\text{s.t.} \quad (24b) - (24g) \quad (29b)$$

$$\quad (28b) - (28f). \quad (29c)$$

There is a large body of research on solution methods for bilinear problems that can be used in solving model (29). Furthermore, many optimization solvers, such as MOSEK (ApS 2022) and
The restricted decomposable NA dual problem can be solved using the same methods developed in the previous sections for the non-decomposable problem. Due to the relaxation of the nonan-ticipativity constraints on z variables, for given $\mathbb{P}$ and $\lambda^y(\cdot)$, subproblem $\mathcal{L}_{LR}^{DNA}(\cdot)$ is a relaxation of $\mathcal{L}_{LR}^{NA}(\cdot)$. Therefore $\nu_R^{DNA} \le \nu_R^{NA}$, i.e., $\nu_R^{NA}$ is a potentially stronger bound. However, the fact that $\mathcal{L}_{LR}^{DNA}(\cdot)$ is decomposable is highly desirable. Particularly for continuous uncertainty sets where we rely on sampling, and the quality of the bound varies based on the selected sample. Since the decomposable model can afford samples of larger sizes it can potentially yield better bounds com-pared to the bound obtained from non-decomposable model over a smaller sample. We explore this possibility in our numerical section.

# 4. Numerical Experiments

In this section we evaluate the performance of the proposed bounding framework over multistage versions of three classical decision-making problems under uncertainty: (i) the newsvendor problem, (ii) the location-transportation problem, and (iii) the capital budgeting problem. Depending on their structure, each problem is solved by using the appropriate models and methods described in Sections 2 and 3, illustrating the applicability of the developed concepts to a large array of problem classes.

## 4.1. Benchmarks and Implementation Details

To assess the quality of the primal and dual bounds, we measure the relative distance of the bounds from the true optimal value when an exact solution of MSARO is available (in small-size instances). Otherwise, we report the optimality gap between the bounds obtained from the proposed methods, and compare it against a gap obtained via traditional bounding methods if one exists. In continuous problems, we consider LDRs (i.e., their application to all decision variables) as the benchmark for the primal decision rules. On the dual side, we use the perfect information (PI) bound (denoted by $\nu^{PI}$) for comparison, which often can be conveniently evaluated for a general MSARO problem. The PI bound corresponds to the optimal objective value of the MSARO problem reformulated as in (10) without the nonanticipativity constraints (10e), i.e., it finds the cost of every scenario in the uncertainty set individually, and then selects the one with the worst-case cost.

The algorithms are implemented in Python and use the Gurobi Optimizer 9.5.1 (Gurobi Optimization, LLC 2022) as the MIP/bilinear program solver. The computational experiments are carried out on the Niagara supercomputer servers (Loken et al. 2010, Ponce et al. 2019). The programs stop if the running time reaches 10 hours, at which point we report the best valid observed lower/upper bound on the instance. As a common design choice for the basis functions of the LDRs, we use the uncertain parameters themselves, i.e., the standard basis functions. Any implementation nuances and enhancements used for improving the performance of the algorithms are discussed for each problem class in a dedicated section, along with the characteristics of the studied instances.
Gurobi (Gurobi Optimization, LLC 2022), offer off-the-shelf alternatives to solve problems of type (29). Further, (29) can be solved heuristically to obtain a valid dual bound. For instance, one could alternate between optimizing over variables $\rho$ and $(\alpha, \beta, \gamma)$ (given that all bilinear terms involve variables $\rho$) to converge towards a local optimal solution.

### 4.4. Restricted Decomposable NA Dual

In solving problem (24) using the cutting-plane method, we frequently optimize (25) to compute $Q(\beta)$. This can become computationally demanding when there is a large number of realizations. Further, we cannot decompose $Q(\beta)$ by realizations, as they are linked through the $z$ variables. In this section, we present an alternative NA reformulation of the MSARO problem that can yield, through its associated DO problem, a potentially weaker bound than the one provided by (24). However, this alternative reformulation offers a computational advantage since in the framework of the cutting-plane method it leads to decomposable subproblems when calculating $Q(\beta)$. To this end, in addition to the decision variable copies $y_t(\xi^T)$, we introduce copy variables $z(\xi^T)$ and explicitly enforce them to be equal via nonanticipativity constraints. For an assigned probability measure $\mathbb{P} \in \mathcal{P}^>$, the alternative NA reformulation of the MSARO problem (2) is:

$$ \min \quad \mathbb{E}_{\xi^T \sim \mathbb{P}} [z(\xi^T)] \qquad (30a) $$

$$ \text{s.t.} \quad \sum_{t \in [T]} c_t (\xi^t)^\top y_t (\xi^T) \le z(\xi^T) \qquad \xi^T \in \Xi^T \qquad (30b) $$

$$ z(\xi^T) = \mathbb{E}_{\xi^T \sim \mathbb{P}} [z(\xi^T)] \qquad \xi^T \in \Xi^T \qquad (30c) $$

$$ (12c), (12d), (12f), (13). \qquad (30d) $$

Together, (30b) and (30c) capture the semantics of the worst-case outcome, which is minimized in the objective function (30a). We remark that constraints (30c) are obtained from individual nonanticipativity constraints similarly to the derivation of constraints (13) provided in Lemma 1.

Relaxation of the nonanticipativity constraints (13) and (30c) with the Lagrangian multipliers $\lambda_t^y(\cdot)$ and $\lambda_t^z(\cdot)$ such that $\mathbb{E}_{\xi^T \sim \mathbb{P}}[\lambda_t^y(\xi^T)] < +\infty$ and $\mathbb{E}_{\xi^T \sim \mathbb{P}}[\lambda_t^z(\xi^T)] < +\infty$, leads to the decomposable NA Lagrangian dual problem $\mathcal{L}^{\text{DNA}}(\mathbb{P}) = \max_{\lambda^y(\cdot), \lambda^z(\cdot)} \mathcal{L}_{\text{LR}}^{\text{DNA}}(\mathbb{P}, \lambda_1^y(\cdot), \dots, \lambda_T^y(\cdot), \lambda^z(\cdot))$, where

$$ \mathcal{L}_{\text{LR}}^{\text{DNA}}(\mathbb{P}, \lambda_1^y(\cdot), \dots, \lambda_T^y(\cdot), \lambda^z(\cdot)) = \min \quad \mathbb{E}_{\xi^T \sim \mathbb{P}} \left[ \left(1 + \lambda^z(\xi^T) - \mathbb{E}_{\xi^T \sim \mathbb{P}} \left[ \lambda^z(\xi^T) \right] \right) z(\xi^T) \right] + \\ \sum_{t \in [T]} \mathbb{E}_{\xi^T \sim \mathbb{P}} \left[ \left( \lambda_t^y(\xi^T) - \mathbb{E}_{\xi^T \sim \mathbb{P}} \left[ \lambda_t^y(\xi^T) \middle| \xi^t = \xi^z \right] \right) \right]^\top y_t(\xi^T) \\ \text{s.t.} \quad (z(\xi^T), y_1(\xi^T), \dots, y_T(\xi^T)) \in Y(\xi^T) \qquad \xi^T \in \Xi^T, $$

with $Y(\xi^T)$ the scenario feasibility space described by constraints (30b), (12c), (12d), and (12f). We note that, in the first expectation term in the objective function the inner expectation is not
## 4.2. Robust Budgeted Newsvendor Problem

In this section, we extend the two-stage newsvendor problem studied in Xu and Hanasusanto (2021) to the multistage setting. In this problem, a decision-maker (the newsvendor) needs to order from a set of items to be sold (only) at the next decision stage, with the objective of maximizing the worst-case profit over the planning horizon. Let $d_{it}(\xi^t)$ be the uncertain demand of item $i \in [I]$ at stage $t \in [2, T]$, $c_i$ and $s_i$ the purchase and shortage costs of item $i$, respectively, and $r_i$ its sale price. To meet the customers' demands of stage $t$, at stage $t-1$ the decision-maker decides on the amounts to be ordered from each item, such that the total spending over the $T$ stages does not exceed a predetermined budget of $B$. Denote by $x_{it}(\xi^t)$ the decision variable for the amount of item $i$ ordered at stage $t \in [T-1]$. The multistage multi-item budgeted newsvendor problem is:

$$ \max z \qquad (30a) $$

$$ \text{s.t.} \quad z \le \sum_{i \in [I]} \sum_{t \in [2,T]} y_{it}(\xi^t) \qquad \xi^T \in \Xi \qquad (30b) $$

$$ y_{it}(\xi^t) \le (r_i - c_i)x_{i,t-1}(\xi^{t-1}) - r_i(x_{i,t-1}(\xi^{t-1}) - d_i(\xi^t)) \qquad i \in [I], t \in [2,T], \xi^t \in \Xi^t \qquad (30c) $$

$$ y_{it}(\xi^t) \le (r_i - c_i)x_{i,t-1}(\xi^{t-1}) - s_i(d_i(\xi^t) - x_{i,t-1}(\xi^{t-1})) \qquad i \in [I], t \in [2,T], \xi^t \in \Xi^t \qquad (30d) $$

$$ \sum_{i \in [I]} \sum_{t \in [T-1]} x_{it}(\xi^t) \le B \qquad \xi^{T-1} \in \Xi^{T-1} \qquad (30e) $$

$$ x_t(\xi^t) \in \mathbb{R}_+^I \qquad t \in [T-1], \xi^t \in \Xi^t. \qquad (30f) $$

where $y_{it}(\xi^t)$ is an auxiliary decision variable capturing the profit from item $i$ at stage $t$, by means of the constraints (30c) and (30d). Constraints (30e) impose a budget of $B$ over the order amounts throughout the planning horizon. The objective function is the worst-case profit of the newsvendor, modeled via (30a) and (30b).

### 4.2.1. Problem Instances

Our instance generation loosely follows the procedure described in Ardestani-Jaafari and Delage (2021) for the two-stage robust newsvendor problem. The deterministic parameters $r_i$, $s_i$ and $c_i$ are drawn uniformly from the intervals $[140, 160]$, $[80, 90]$ and $[50, 70]$, respectively. We consider a discrete uncertainty set modeled as a stagewise-dependent scenario tree with branching factor BR (i.e., every node of the tree prior to the leaves has BR many child nodes). The demand realizations $d_{it}(\xi^t)$, $i \in [I]$, $t \in [2, T]$ at a child node are drawn uniformly from $[\mu_{it} - \sigma_{it}, \mu_{it} + \sigma_{it}]$, where $\mu_{it}$ and $\sigma_{it}$ are deterministic parameters uniformly drawn from the intervals $[20, 40]$ and $[10, 20]$. We have generated 26 small-size instances with $T \in [3, 5]$, $I \in [2, 5]$, $B \in \{100, 150, 200, 250, 300\}$, and $BR \in \{2, 3, 4, 5, 10\}$, such that the number of scenarios $|\Xi| = BR^{T-1}$ is less than 150. Additionally, we have generated 18 large-size instances with $T \in [4, 8]$. For $T=4$, the number of items $I$ lies in the set $[3, 5]$, with a budget $B \in \{200, 300\}$ and $BR \in \{10, 15, 20\}$. For $T \in [5, 8]$, our instances have $I \in \{3, 4\}$ items, budget of $B \in \{300, 400, 500, 600\}$, and $BR \in [3, 6]$, restricted to the cases with $|\Xi| \le 3000$.
conditional, since constraints (30c) (and accordingly their associated dual functions) are not defined at every stage $t \in [T]$ as decision variable $z(\xi^T) \in \mathbb{R}$ captures the cost of an entire realization $\xi^T$. Since in $\mathcal{L}_{LR}^{DNA}(\cdot)$, the objective function, constraints and variables are decomposable in $\xi^T$ they can be optimized individually.

In deriving the decomposable NA Lagrangian relaxation problem $\mathcal{L}_{LR}^{DNA}(\cdot)$ we apply Lemma 1 of (Daryalal et al. 2024) to obtain both expectation terms. We remark that the lemma requires the condition $\mathbb{E}_{\xi^T \sim \mathbb{P}}[\text{diam}(\text{proj}_{z(\cdot)}Y(\xi^T))] < +\infty$ which is not naturally satisfied. However, as a result of Assumptions 1-3 and the compactness of the uncertainty set $\Xi^T$ the optimal value of the MSARO problem is bounded. As such the functionals $z(\cdot)$ can be artificially bounded without changing the optimal value of $\mathcal{L}_{LR}^{DNA}(\cdot)$.

After substituting the decision rules $\lambda_t^y(\xi^T) = \Psi_t(\xi^T)\alpha_t^y, t \in [T]$ and $\lambda^z(\xi^T) = \Psi_T(\xi^T)\alpha^z$ in $\mathcal{L}^{DNA}(\mathbb{P})$, where $\alpha^z \in \mathbb{R}^{K_T}$, and merging with the optimization over the probability distributions $\mathbb{P} \in \mathcal{P}^\ge$ the decomposable DO problem is:

$$ \bar{\nu}_R^{DNA-DO} := \max_{\mathbb{P} \in \mathcal{P}^\ge, \alpha_1^y, \dots, \alpha_T^y, \alpha^z} \mathcal{L}_{LR}^{DNA}(\mathbb{P}, \Psi_1(\xi^T)\alpha_1^y, \dots, \Psi_T(\xi^T)\alpha_T^y, \Psi_T(\xi^T)\alpha^z). $$

This DO problem can be solved using the same methods developed in the previous sections for the non-decomposable DO problem. Due to the relaxation of the nonanticipativity constraints on $z(\cdot)$ variables, for given $\mathbb{P}$ and $\lambda^y(\cdot)$, subproblem $\mathcal{L}_{LR}^{DNA}(\cdot)$ is a relaxation of $\mathcal{L}_{LR}^{NA}(\cdot)$. Therefore $\bar{\nu}_R^{DNA-DO} \le \bar{\nu}_R^{NA-DO}$, i.e., $\bar{\nu}_R^{NA-DO}$ is a potentially stronger bound. However, the fact that $\mathcal{L}_{LR}^{DNA}(\cdot)$ is decomposable is highly desirable. Particularly for continuous or large discrete uncertainty sets where we rely on sampling, and the quality of the bound varies based on the selected sample. Since the decomposable model can afford samples of larger sizes it can potentially yield better bounds compared to the bound obtained from non-decomposable model over a smaller sample. We explore this trade-off in our numerical section.

# 5. Numerical Experiments

We evaluate the performance of the proposed bounding framework over multistage versions of three classical decision-making problems under uncertainty: (i) the newsvendor problem, (ii) the location-transportation problem, and (iii) the capital budgeting problem. Depending on their structure, each problem is solved by using the appropriate models and methods described in Sections 3 and 4, illustrating the applicability of the developed concepts to a large array of problem classes.

## 5.1. Benchmarks and Implementation Details

To assess the quality of the primal and dual bounds, we measure the relative distance of the bounds from the true optimal value when an exact solution of MSARO is available (in small-size instances).
**4.2.2. Quality of the Bounds** The small-size instances can be conveniently optimized by solving model (10) over all the scenarios in the uncertainty set. Using this, we can examine the quality of a primal/dual bound by measuring its relative distance to the optimal objective value $\nu^*$. For this problem, all primal and dual problems are solved by the extensive form (i.e., monolithic) of their respective models. Figures 6 and 7 present the gap between the bound and the optimal value of the exact solution, defined as $100(\frac{\nu^* - \nu^{(\cdot)}}{\nu^*})$ and $100(\frac{\nu^{(\cdot)} - \nu^*}{\nu^*})$ for primal and dual bounds, respectively, and presented as a percentage (detailed results are given in the e-companion). In each figure, the solid bars depict the performance of the newly proposed bounds, while the hatched bars represent the benchmarks. The results show that $\nu^{2S-LDR}$ and $\bar{\nu}_R^{NA}$ outperform the benchmark bounds by orders of magnitude. More precisely, $\nu^{2S-LDR}$ on average achieves 84% improvement

<figure><img src="image_7.png" /><figcaption>**Figure 6** Quality of the bounds from LDRs and two-stage LDRs for small-size newsvendor instances</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>**Figure 7** Quality of the bounds from the PI and LDDRs for small-size newsvendor instances</figcaption></figure>

over $\nu^{LDR}$, with reductions in relative distance ranging from 64% to 98%. Interestingly, the quality
Otherwise, we report the optimality gap between the bounds obtained from the proposed methods, and compare it against a gap from traditional bounding methods if one exists. In problems with continuous recourse, we consider LDRs (i.e., their application to all decision variables) as the benchmark for the primal decision rules. On the dual side, we use the perfect information (PI) bound (denoted by $v^{PI}$) for comparison, which often can be conveniently evaluated for a general MSARO problem, as well as the bound obtained by solving model (2) using only the binding realizations identified from the primal decision rule solution. The PI bound corresponds to the optimal objective value of the MSARO problem reformulated as in (12) without the nonanticipativity constraints (12e), i.e., it finds the cost of every realization in the uncertainty set individually, and then selects the one with the worst-case cost.

The algorithms are implemented in Python and use the Gurobi Optimizer 9.5.1 (Gurobi Optimization, LLC 2022) as the mixed-integer/bilinear programming solver. The computational experiments are carried out on the Niagara supercomputer servers (Loken et al. 2010, Ponce et al. 2019). For instances with discrete uncertainty sets, the programs have a time limit of 1 hour. We report the valid lower/upper bound at the point of termination. For instances with continuous uncertainty sets, this time limit is extended to 10 hours. As a common design choice for the basis functions of the LDRs, we use the uncertain parameters themselves, i.e., the standard basis functions. Any implementation nuances and enhancements used for improving the performance of the algorithms are discussed for each problem class in a dedicated section, along with the characteristics of the studied instances.

## 5.2. Robust Newsvendor Problem

In this section, we extend the two-stage newsvendor problem studied in (Xu and Hanasusanto 2021) to the multistage setting. In this problem, a decision-maker (the newsvendor) needs to order from a set of items to be sold (only) at the next decision stage, with the objective of maximizing the worst-case profit over the planning horizon. Let $d_{it}(\xi^t)$ be the uncertain demand of item $i \in [I]$ at stage $t \in [2, T]$, $c_i$ and $s_i$ the purchase and shortage costs of item $i$, respectively, and $r_i$ its sale price. To meet the customers' demands of stage $t$, at stage $t-1$ the decision-maker decides on the amounts to be ordered from each item, such that the total spending over the $T$ stages does not exceed a predetermined budget of $B$. Denote by $x_{it}(\xi^t)$ the decision variable for the amount of item $i$ ordered at stage $t \in [T-1]$. The multistage multi-item budgeted newsvendor problem is:

$$ \max z \qquad (32a) $$

$$ \text{s.t.} \quad z \le \sum_{i \in [I]} \sum_{t \in [2,T]} y_{it}(\xi^t) \qquad (32b) $$
of the bound $\nu^{2S-LDR}$ remains rather stable with changes in the number of items $I$ and budget $B$, compared to the drastic changes of $\nu^{LDR}$ with variations in the inputs. The notable performance of the two-stage LDRs for the newsvendor problem contrasted with the LDRs can be explained by the nature of the recourse variables $y$. In model (30), $y_{it}(\xi^t)$ variables determine the profit of a given realization at stage $t$ for item $i$, which for the newsvendor problem is by definition nonlinear. In fact, $y_{it}(\xi^t)$ is an auxiliary variable, defined to linearize the following net profit at stage $t$ from item $i$:

$$r_i \min \{x_{i,t-1}(\xi^{t-1}), d_{it}(\xi^t)\} - c_i x_{i,t-1}(\xi^{t-1}) - s_i \max \{d_{it}(\xi^t) - x_{i,t-1}(\xi^{t-1}), 0\}.$$

Therefore, for the newsvendor problem LDRs always return suboptimal decisions as they restrict the form of the nonlinear profit function to be affine, while two-stage LDRs allow them to take any form, giving them an immediate advantage over LDRs.

From the dual perspective, $\bar{\nu}_R^{NA}$ achieves an average improvement of 55% compared to $\nu^{PI}$, and in 9 instances fully closes the gap. From the results of Figure 7, a common observation is that LDDRs return a bound of higher quality for smaller values of the ratio $\frac{B}{(T-1) \times I}$, which is an estimate of the average budget available per product at every stage. This trend suggests that the restricted NA dual bound performs better with tighter budget, when there is a higher dependency between stages, in which case the importance of intermediate decisions becomes more pronounced.

### 4.2.3. Optimality Gap
For instances with larger number of scenarios, we compare the optimality gap from the benchmark methods ($OPT^T$) with the gap obtained by applying the newly proposed methods, namely two-stage LDRs and LDDRs ($OPT^N$):

$$OPT^T = 100 \left( \frac{\nu^{PI} - \nu^{LDR}}{\nu^{LDR}} \right), \quad OPT^N = 100 \left( \frac{\bar{\nu}_R^{NA} - \nu^{2S-LDR}}{\nu^{2S-LDR}} \right).$$

Table 1 presents the optimality gaps from the benchmark and proposed models (whose running times are provided in the e-companion). Achieving an average gap reduction of 82%, there is a considerable value in using the two-stage LDRs and LDDRs in devising policies for the multistage newsvendor problem. The optimality gaps $OPT^N$ range from 4% to 46%, with some interesting cases such as instances 3, 5 and 14 where using LDR policies can even lead to a loss of profit in the worst-case. In the majority of the instances, both two-stage LDRs and LDDRs contribute to the improvement of the gap, although the primal side clearly has the larger impact, to the point that in five instances, the dual bound $\bar{\nu}_R^{NA}$ and the benchmark $\nu^{PI}$ are the same. In some instances, the PI bound might already be strong enough so that using LDDRs does not make a tangible difference in strengthening the bound. This is, for instance, the case for instances 8, 12 and 18 that attain a small optimality gap without significantly improving $\nu^{PI}$. As discussed for the small-size instances,
$$y_{it}(\xi^t) \le (r_i - c_i)x_{i,t-1}(\xi^{t-1}) - r_i(x_{i,t-1}(\xi^{t-1}) - d_i(\xi^t)) \quad i \in [I], t \in [2, T], \xi^t \in \Xi^t \quad (32c)$$

$$y_{it}(\xi^t) \le (r_i - c_i)x_{i,t-1}(\xi^{t-1}) - s_i(d_i(\xi^t) - x_{i,t-1}(\xi^{t-1})) \quad i \in [I], t \in [2, T], \xi^t \in \Xi^t \quad (32d)$$

$$\sum_{i \in [I]} \sum_{t \in [T-1]} x_{it}(\xi^t) \le B \qquad \xi^{T-1} \in \Xi^{T-1} \quad (32e)$$

$$x_t(\xi^t) \in \mathbb{R}_+^I \qquad t \in [T-1], \xi^t \in \Xi^t. \quad (32f)$$

where auxiliary variable $y_{it}(\xi^t)$ captures the profit from item $i$ at stage $t$, by means of (32c)-(32d). Constraints (32e) impose a budget of $B$ over the order amounts throughout the planning horizon. The objective function is the worst-case profit of the newsvendor, modeled via (32a) and (32b).

### 5.2.1. Problem Instances
Our instance generation loosely follows the procedure described in (Ardestani-Jaafari and Delage 2021) for the two-stage robust newsvendor problem. Parameters $r_i$, $s_i$ and $c_i$ are drawn uniformly from the intervals $[140, 160]$, $[80, 90]$ and $[50, 70]$, respectively. We consider a discrete uncertainty set modeled as a stagewise-dependent scenario tree with branching factor BR (i.e., every node of the tree prior to the leaves has BR many child nodes). Demand realizations $d_{it}(\xi^t)$, $i \in [I]$, $t \in [2, T]$ at a child node are drawn uniformly from $[\mu_{it} - \sigma_{it}, \mu_{it} + \sigma_{it}]$, where $\mu_{it}$ and $\sigma_{it}$ are uniformly drawn from the intervals $[20, 40]$ and $[10, 20]$. We have generated 26 small-size instances with $T \in [3, 5]$, $I \in [2, 5]$, $B \in \{100, 150, 200, 250, 300\}$, and $BR \in \{2, 3, 4, 5, 10\}$, such that the number of realizations $|\Xi^T| = BR^{T-1}$ is less than 150. Additionally, we have generated 18 large-size instances with $T \in [4, 8]$. For $T=4$, the number of items $I$ lies in the set $[3, 5]$, with a budget $B \in \{200, 300\}$ and $BR \in \{10, 15, 20\}$. For $T \in [5, 8]$, our instances have $I \in \{3, 4\}$ items, budget of $B \in \{300, 400, 500, 600\}$, and $BR \in [3, 6]$, restricted to the cases with $|\Xi^T| \le 3000$.

### 5.2.2. Quality of the Bounds
The small-size instances are easily optimized by solving model (12) over all realizations in the uncertainty set. In our case, this computation takes less than 3 seconds. Using these optimal values, we can examine the quality of a primal/dual bound by measuring its relative distance to the optimal objective value $\nu^*$. For this problem, all primal and dual problems are solved by the extensive form (i.e., monolithic) of their respective models. Figures 5 and 6 present the gap between the bound and the optimal value of the exact solution, defined as $100(\frac{\nu^* - \nu^{(t)}}{\nu^*})$ and $100(\frac{\nu^{(t)} - \nu^*}{\nu^*})$ for primal and dual bounds, respectively, and presented as a percentage (detailed results along with solution times are given in Appendix D.1). In each figure, the solid bars depict the performance of the newly proposed bounds, while the hatched bars represent the benchmarks. The results show that $\nu^{2S-LDR}$ and $\bar{\nu}_R^{NA-DO}$ outperform the benchmark bounds by orders of magnitude. More precisely, $\nu^{2S-LDR}$ on average achieves 84% improvement over $\nu^{LDR}$, with reductions in relative distance ranging from 64% to 98%. Interestingly, the quality of the bound $\nu^{2S-LDR}$ remains rather stable with changes in the number of items $I$ and budget $B$, compared to the drastic changes of $\nu^{LDR}$ with variations in the inputs. The notable performance of
<table><caption>**Table 1** Optimality gaps for larger instances of the newsvendor problem</caption><thead><tr><th rowspan="2">Instance</th><th rowspan="2">T</th><th rowspan="2">BR</th><th rowspan="2">|Ξ|</th><th rowspan="2">I</th><th rowspan="2">B</th><th colspan="2">Primal Bounds</th><th colspan="2">Dual Bounds</th><th colspan="2">Optimality Gap</th><th rowspan="2">Gap reduction</th></tr><tr><th>ν<sup>LDR</sup></th><th>ν<sup>2S</sup></th><th>ν<sup>NA</sup><sub>R</sub></th><th>ν<sup>PI</sup></th><th>OPT<sup>N</sup></th><th>OPT<sup>T</sup></th></tr></thead><tbody><tr><td>1</td><td></td><td></td><td></td><td>3</td><td>200</td><td>5648.8</td><td>8142.9</td><td>8876.4</td><td>9353.0</td><td>9.0%</td><td>65.6%</td><td>86.3%</td></tr><tr><td>2</td><td></td><td></td><td></td><td>3</td><td>300</td><td>9143.5</td><td>13687.6</td><td>17853.0</td><td>17853.0</td><td>30.4%</td><td>95.3%</td><td>68.1%</td></tr><tr><td>3</td><td>4</td><td>10</td><td>1000</td><td>4</td><td>200</td><td>-104.4</td><td>440.0</td><td>642.4</td><td>919.0</td><td>46.0%</td><td>980.4%</td><td>95.3%</td></tr><tr><td>4</td><td></td><td></td><td></td><td>4</td><td>300</td><td>11029.5</td><td>15854.7</td><td>18400.7</td><td>18432.0</td><td>16.1%</td><td>67.1%</td><td>76.1%</td></tr><tr><td>5</td><td></td><td></td><td></td><td>5</td><td>300</td><td>724.9</td><td>6125.5</td><td>6748.6</td><td>7368.0</td><td>10.2%</td><td>916.5%</td><td>98.9%</td></tr><tr><td>6</td><td></td><td>4</td><td>15</td><td>3375</td><td>3</td><td>200</td><td>5111.3</td><td>7072.0</td><td>7646.6</td><td>8222.0</td><td>8.1%</td><td>60.9%</td><td>86.7%</td></tr><tr><td>7</td><td></td><td>4</td><td>20</td><td>8000</td><td>3</td><td>300</td><td>9030.9</td><td>13040.7</td><td>17615.0</td><td>17615.0</td><td>35.1%</td><td>95.1%</td><td>63.1%</td></tr><tr><td>8</td><td></td><td>5</td><td>5</td><td>625</td><td>3</td><td>300</td><td>11134.5</td><td>15053.4</td><td>16191.6</td><td>16680.0</td><td>7.6%</td><td>49.8%</td><td>84.8%</td></tr><tr><td>9</td><td></td><td></td><td></td><td></td><td>4</td><td>300</td><td>3651.5</td><td>9494.8</td><td>10553.8</td><td>10628.0</td><td>11.2%</td><td>191.1%</td><td>94.2%</td></tr><tr><td>10</td><td></td><td>5</td><td>6</td><td>1296</td><td>3</td><td>400</td><td>15013.4</td><td>20271.7</td><td>25157.0</td><td>25157.0</td><td>24.1%</td><td>67.6%</td><td>64.3%</td></tr><tr><td>11</td><td></td><td></td><td></td><td></td><td>3</td><td>400</td><td>15492.3</td><td>21457.9</td><td>28163.0</td><td>28163.0</td><td>31.2%</td><td>81.8%</td><td>61.8%</td></tr><tr><td>12</td><td>6</td><td>4</td><td>1024</td><td>4</td><td>400</td><td>7124.3</td><td>14805.2</td><td>15426.2</td><td>15473.0</td><td>4.2%</td><td>117.2%</td><td>96.4%</td></tr><tr><td>13</td><td></td><td></td><td></td><td></td><td>4</td><td>500</td><td>14445.4</td><td>24887.7</td><td>32946.9</td><td>32973.0</td><td>32.4%</td><td>128.3%</td><td>74.8%</td></tr><tr><td>14</td><td></td><td></td><td></td><td></td><td>3</td><td>300</td><td>24.3</td><td>1495.9</td><td>2137.3</td><td>2383.0</td><td>42.9%</td><td>9692.2%</td><td>99.6%</td></tr><tr><td>15</td><td>7</td><td>3</td><td>729</td><td>3</td><td>400</td><td>12994.7</td><td>17774.1</td><td>19554.5</td><td>19983.0</td><td>10.0%</td><td>53.8%</td><td>81.4%</td></tr><tr><td>16</td><td></td><td></td><td></td><td></td><td>4</td><td>400</td><td>3267.5</td><td>4303.1</td><td>5354.3</td><td>5427.0</td><td>24.4%</td><td>66.1%</td><td>63.0%</td></tr><tr><td>17</td><td>8</td><td>3</td><td>2187</td><td>3</td><td>500</td><td>16149.1</td><td>25321.4</td><td>30259.0</td><td>30259.0</td><td>19.5%</td><td>87.4%</td><td>77.7%</td></tr><tr><td>18</td><td></td><td></td><td></td><td></td><td>4</td><td>600</td><td>16892.3</td><td>27608.1</td><td>28717.7</td><td>28747.0</td><td>4.0%</td><td>70.2%</td><td>94.3%</td></tr></tbody></table>

depending on the number of stages, items, and available budget, in some instances the PI bound can indeed be weak. In such cases $\bar{\nu}_R^{NA}$ can be potentially improved by a better design of the basis functions for LDDRs that better captures the correlation among all items imposed by the budget restriction.

### 4.3. Robust Location-Transportation Problem

The two-stage robust location-transportation problem studied by Zeng and Zhao (2013) is as follows. Given a set of $I$ potential facilities with fixed building cost $f_i$ and unit capacity cost $a_i$, $i \in [I]$, we have to meet the uncertain demand of a set of customers $J$ with unit transportation cost $c_{ij}$, $i \in [I]$, $j \in [J]$. The goal is to decide which facilities to open, as well as their initial capacities, such that the worst-case total cost of facility deployment and future transportation is minimized. We expand the problem to a multistage setting by defining $d_{jt}(\xi^t)$ as the demand of customer $j$ at stage $t$. The MSARO model is then:

$$ \min z \qquad (31a) $$

$$ \text{s.t.} \quad z \ge \sum_{i \in [I]} f_i y_i + \sum_{t \in [T]} \sum_{i \in [I]} a_i s_{it}(\xi^t) + \sum_{t \in [2,T]} \sum_{i \in [I]} \sum_{j \in [J]} c_{ij} x_{ijt}(\xi^t) \quad \xi^T \in \Xi \qquad (31b) $$

$$ s_{i1} \le K_i y_i \qquad i \in [I] \qquad (31c) $$

$$ s_{it}(\xi^t) = s_{i,t-1}(\xi^{t-1}) - \sum_{j \in [J]} x_{ijt}(\xi^t) \qquad i \in [I], t \in [2, T], \xi^t \in \Xi^t \qquad (31d) $$
<figure><img src="image_7.png" /><figcaption>Figure 5 Quality of the primal bounds from LDRs and two-stage LDRs for small-size newsvendor instances</figcaption></figure>

the two-stage LDRs for the newsvendor problem contrasted with the LDRs can be explained by the nature of the recourse variables. In model (32), $y_{it}(\xi^t)$ determines the profit of a given realization at stage $t$ for item $i$, which for the newsvendor problem is by definition nonlinear. In fact, $y_{it}(\xi^t)$ is an auxiliary variable, defined to linearize the following net profit at stage $t$ from item $i$:

$$r_i \min \{x_{i,t-1}(\xi^{t-1}), d_{it}(\xi^t)\} - c_i x_{i,t-1}(\xi^{t-1}) - s_i \max \{d_{it}(\xi^t) - x_{i,t-1}(\xi^{t-1}), 0\}.$$

Therefore, for the newsvendor problem LDRs always return suboptimal decisions as they restrict the form of the nonlinear profit function to be affine, while two-stage LDRs allow them to take any form, giving them an immediate advantage over LDRs.

From the dual perspective, $\bar{\nu}_R^{\text{NA-DO}}$ achieves an average improvement of 55% compared to $\nu^{\text{PI}}$, and in 9 instances fully closes the gap. From the results of Figure 6, a common observation is that LDDRs return a bound of higher quality for smaller values of the ratio $\frac{B}{(T-1) \times I}$, which is an

<figure><img src="image_8.png" /><figcaption>Figure 6 Quality of the dual bounds from the PI and LDDRs for small-size newsvendor instances</figcaption></figure>
$$ \sum_{i \in [I]} x_{ijt}(\xi^t) \geq d_{jt}(\xi^t) \qquad j \in [J], t \in [2, T], \xi^t \in \Xi^t \quad (31e) $$

$$ y \in \{0, 1\}^I, s_1 \in \mathbb{R}_+^I \qquad (31f) $$

$$ s_t(\xi^t) \in \mathbb{R}_+^I, x_t(\xi^t) \in \mathbb{R}_+^{I \times J} \qquad t \in [2, T], \xi^t \in \Xi^t, \quad (31g) $$

where $y_i$ is a binary variable equal to 1 if facility $i$ is built, $s_{i1}$ determines the initial capacity of facility $i$, $s_{it}(\xi^t)$ is the state variable calculating the remaining capacity of facility $i$ at stage $t$, and $x_{ijt}(\xi^t)$ is the amount of goods transported from facility $i$ to customer $j$ at stage $t$. Objective function (31a) together with constraint (31b) measures the worst-case cost. Constraints (31c) bound the initial capacity of the facilities, whereas constraints (31d) are the state equations calculating the remaining capacities. Constraints (31e) ensure that the customer demands are met.

### 4.3.1. Problem Instances

Our instances are generated using the instance parameters described in Zeng and Zhao (2013) for the two-stage problem. For number of stages $T \in [3, 5]$, we have five combinations for the number of facilities and customers, $(I, J) \in \{(5, 5), (5, 7), (5, 10), (10, 10), (20, 20)\}$. Fixed installation, unit capacity, and unit transportation costs are drawn from $f_i \in [100, 1000], a_i \in [10, 100], c_{ij} \in [1, 1000],$ respectively. Maximal capacity is set to $K_i = 2 \times 10^4$ based on preliminary experiments to make sure the instances are feasible. Customer demands are random parameters with support $[\mu_{jt}, (1 + \alpha^d)\mu_{jt}]$, where $\mu_{jt} \in [10, 500]$ and $\alpha^d \in \{0.1, 0.3, 0.5\}$ are given, with $\alpha^d$ a parameter controlling the variation among demand realizations of customer $j$. For this problem, we have generated 32 instances with $\alpha^d = 0.5$ over small scenario trees to compare the bounds with the optimal objective value. To demonstrate the versatility of our solution methods, in addition to the described scenario-tree instances, we have generated 85 instances where the demands belong to the following budgeted uncertainty set:

$$ \Xi = \left\{ d \in \mathbb{R}_+^{J \times T-1} \mid d_{jt} = \mu_{jt} + \delta_{jt}\sigma_{jt}, \delta_{jt} \in [0, 1], j \in [J], t \in [2, T], \sum_{t \in [2, T]} \sum_{j \in [J]} \delta_{jt} \leq \Gamma \right\}. $$

We remark that, $\Gamma$ correlates the demands of all customers and stages together, which results in stagewise (temporal) dependence between the decision stages. In our instances we use $\Gamma = \alpha^u I$, with $\alpha^u \in \{0.1, 0.4, 0.7, 1\}$.

### 4.3.2. Scenario-Tree Instances

In solving the primal (2S-LDR) and exact models of the scenario-tree instances, we have used their respective extensive forms, while for computing the restricted NA dual bound we implemented the cutting-plane method described in Section 3.3. Table 2 presents the gap between the bound and the optimal value of the exact solution, $100(\frac{\nu^{(\cdot)} - \nu^*}{\nu^*})$ and $100(\frac{\nu^* - \nu^{(\cdot)}}{\nu^*})$ for the primal and dual problems, respectively. Results show that, among the group of instances with similar characteristics, as the branching factor or the number of stages increases, the
estimate of the average budget available per product at every stage. This trend suggests that the restricted NA dual bound performs better with tighter budget, when there is a higher dependency between stages, in which case the importance of intermediate decisions becomes more pronounced.

**5.2.3. Optimality Gap** For instances with larger number of realizations, we compare the optimality gap from the benchmark methods (OPT<sup>T</sup>) with the gap obtained by applying the newly proposed methods, namely two-stage LDRs and LDDRs (OPT<sup>N</sup>):

$$ \text{OPT}^{\text{T}} = 100 \left( \frac{\nu^{\text{PI}} - \nu^{\text{LDR}}}{\nu^{\text{LDR}}} \right), \quad \text{OPT}^{\text{N}} = 100 \left( \frac{\bar{\nu}_{R}^{\text{NA-DO}} - \nu^{\text{2S-LDR}}}{\nu^{\text{2S-LDR}}} \right). $$

Table 1 presents the optimality gaps from the benchmark and proposed models (whose running times are provided in Appendix D.1). Achieving an average gap reduction of 81.6%, there is a

<table><caption>**Table 1** Optimality gaps for larger instances of the newsvendor problem</caption><thead><tr><th rowspan="2">Instance</th><th rowspan="2">T</th><th rowspan="2">BR</th><th rowspan="2">|Ξ<sup>T</sup>|</th><th rowspan="2">I</th><th rowspan="2">B</th><th colspan="2">Primal Bounds</th><th colspan="2">Dual Bounds</th><th colspan="2">Optimality Gap</th><th rowspan="2">Gap reduction</th></tr><tr><th>ν<sup>LDR</sup></th><th>ν<sup>2S-LDR</sup></th><th>ν̄<sup>NA-DO</sup></th><th>ν<sup>PI</sup></th><th>OPT<sup>N</sup></th><th>OPT<sup>T</sup></th></tr></thead><tbody><tr><td>1</td><td></td><td></td><td></td><td>3</td><td>200</td><td>5648.8</td><td>8142.9</td><td>9002.8</td><td>9353.0</td><td>10.6%</td><td>65.6%</td><td>83.9%</td></tr><tr><td>2</td><td></td><td></td><td></td><td>3</td><td>300</td><td>9143.5</td><td>13687.6</td><td>17853.0</td><td>17853.0</td><td>30.4%</td><td>95.3%</td><td>68.1%</td></tr><tr><td>3</td><td>4</td><td>10</td><td>1000</td><td>4</td><td>200</td><td>-104.4</td><td>440.0</td><td>642.4</td><td>919.0</td><td>46.0%</td><td>980.4%</td><td>95.3%</td></tr><tr><td>4</td><td></td><td></td><td></td><td>4</td><td>300</td><td>11029.5</td><td>15854.7</td><td>18432.0</td><td>18432.0</td><td>16.3%</td><td>67.1%</td><td>75.8%</td></tr><tr><td>5</td><td></td><td></td><td></td><td>5</td><td>300</td><td>724.9</td><td>6125.5</td><td>6614.0</td><td>7368.0</td><td>8.0%</td><td>916.5%</td><td>99.1%</td></tr><tr><td>6</td><td></td><td>4</td><td>15</td><td>3375</td><td>3</td><td>200</td><td>5111.3</td><td>7072.0</td><td>8222.0</td><td>8222.0</td><td>16.3%</td><td>60.9%</td><td>73.3%</td></tr><tr><td>7</td><td></td><td>4</td><td>20</td><td>8000</td><td>3</td><td>300</td><td>9030.9</td><td>13040.7</td><td>17615.0</td><td>17615.0</td><td>35.1%</td><td>95.1%</td><td>63.1%</td></tr><tr><td>8</td><td></td><td>5</td><td>5</td><td>625</td><td>3</td><td>300</td><td>11134.5</td><td>15053.4</td><td>16216.4</td><td>16680.0</td><td>7.7%</td><td>49.8%</td><td>84.5%</td></tr><tr><td>9</td><td></td><td></td><td></td><td></td><td>4</td><td>300</td><td>3651.5</td><td>9494.8</td><td>10403.1</td><td>10628.0</td><td>9.6%</td><td>191.1%</td><td>95.0%</td></tr><tr><td>10</td><td></td><td>5</td><td>6</td><td>1296</td><td>3</td><td>400</td><td>15013.4</td><td>20271.7</td><td>25157.0</td><td>25157.0</td><td>24.1%</td><td>67.6%</td><td>64.3%</td></tr><tr><td>11</td><td></td><td></td><td></td><td>3</td><td>400</td><td>15492.3</td><td>21457.9</td><td>28163.0</td><td>28163.0</td><td>31.2%</td><td>81.8%</td><td>61.8%</td></tr><tr><td>12</td><td>6</td><td>4</td><td>1024</td><td>4</td><td>400</td><td>7124.3</td><td>14805.2</td><td>15473.0</td><td>15473.0</td><td>4.5%</td><td>117.2%</td><td>96.2%</td></tr><tr><td>13</td><td></td><td></td><td></td><td>4</td><td>500</td><td>14445.4</td><td>24887.7</td><td>32973.0</td><td>32973.0</td><td>32.5%</td><td>128.3%</td><td>74.7%</td></tr><tr><td>14</td><td></td><td></td><td></td><td>3</td><td>300</td><td>24.3</td><td>1495.9</td><td>2150.0</td><td>2383.0</td><td>43.7%</td><td>9692.2%</td><td>99.5%</td></tr><tr><td>15</td><td>7</td><td>3</td><td>729</td><td>3</td><td>400</td><td>12994.7</td><td>17774.1</td><td>19554.5</td><td>19983.0</td><td>10.0%</td><td>53.8%</td><td>81.4%</td></tr><tr><td>16</td><td></td><td></td><td></td><td>4</td><td>400</td><td>3267.5</td><td>4303.1</td><td>5118.1</td><td>5427.0</td><td>18.9%</td><td>66.1%</td><td>71.3%</td></tr><tr><td>17</td><td></td><td>8</td><td>3</td><td>2187</td><td>3</td><td>500</td><td>16149.1</td><td>25321.4</td><td>30259.0</td><td>30259.0</td><td>19.5%</td><td>87.4%</td><td>77.7%</td></tr><tr><td>18</td><td></td><td></td><td></td><td></td><td>4</td><td>600</td><td>16892.3</td><td>27608.1</td><td>28747.0</td><td>28747.0</td><td>4.1%</td><td>70.2%</td><td>94.1%</td></tr></tbody></table>

considerable value in using the two-stage LDRs and LDDRs in devising policies for the multistage newsvendor problem. The optimality gaps OPT<sup>N</sup> range from 4% to 46%. There are some interesting cases such as instance 3 where LDR policies can even lead to a profit loss in the worst case. In the majority of the instances, both two-stage LDRs and LDDRs contribute to the improvement of the gap, although the primal side clearly has the larger impact. For instance, in five instances, the dual bound $\bar{\nu}_R^{\text{NA-DO}}$ and the benchmark $\nu^\text{PI}$ are the same. In some instances, such as instances 8, 12 and 18, the PI bound might already be strong enough so that using LDDRs does not make a tangible difference in strengthening it.
<table><caption>Table 2: Quality of the primal and dual bounds for the scenario-tree instances of the location-transportation problem</caption><thead><tr><th>T</th><th>I</th><th>J</th><th>BR</th><th>|Ξ|</th><th>ν<sup>LDR</sup></th><th>ν<sup>2S-LDR</sup></th><th>ν̄<sup>NA</sup><sub>R</sub></th><th>ν<sup>PI</sup></th></tr></thead><tbody><tr><td rowspan="10">3</td><td rowspan="10">5</td><td rowspan="10">10</td><td>3</td><td>9</td><td>2.9%</td><td>0.0%</td><td>1.4%</td><td>4.7%</td></tr><tr><td>4</td><td>16</td><td>0.8%</td><td>0.0%</td><td>3.4%</td><td>3.4%</td></tr><tr><td>5</td><td>25</td><td>4.6%</td><td>0.0%</td><td>1.5%</td><td>13.2%</td></tr><tr><td>6</td><td>36</td><td>5.6%</td><td>0.0%</td><td>17.4%</td><td>27.5%</td></tr><tr><td>7</td><td>49</td><td>8.2%</td><td>0.0%</td><td>9.4%</td><td>9.4%</td></tr><tr><td>8</td><td>64</td><td>8.8%</td><td>0.0%</td><td>5.2%</td><td>14.6%</td></tr><tr><td>9</td><td>81</td><td>6.5%</td><td>0.0%</td><td>6.1%</td><td>6.1%</td></tr><tr><td>10</td><td>100</td><td>6.6%</td><td>0.0%</td><td>8.0%</td><td>10.2%</td></tr><tr><td rowspan="10">3</td><td rowspan="10">10</td><td rowspan="10">10</td><td>3</td><td>9</td><td>2.3%</td><td>0.0%</td><td>2.8%</td><td>13.2%</td></tr><tr><td>4</td><td>16</td><td>2.7%</td><td>0.1%</td><td>1.2%</td><td>8.4%</td></tr><tr><td>5</td><td>25</td><td>6.1%</td><td>0.0%</td><td>12.0%</td><td>19.3%</td></tr><tr><td>6</td><td>36</td><td>9.8%</td><td>0.0%</td><td>10.2%</td><td>22.6%</td></tr><tr><td>7</td><td>49</td><td>13.7%</td><td>0.1%</td><td>6.7%</td><td>14.1%</td></tr><tr><td>8</td><td>64</td><td>13.9%</td><td>0.3%</td><td>5.5%</td><td>5.5%</td></tr><tr><td>9</td><td>81</td><td>12.2%</td><td>0.2%</td><td>8.1%</td><td>15.5%</td></tr><tr><td>10</td><td>100</td><td>12.3%</td><td>0.2%</td><td>10.1%</td><td>19.0%</td></tr><tr><td rowspan="10">3</td><td rowspan="10">20</td><td rowspan="10">20</td><td>3</td><td>9</td><td>2.5%</td><td>0.7%</td><td>9.5%</td><td>11.3%</td></tr><tr><td>4</td><td>16</td><td>3.9%</td><td>0.5%</td><td>7.0%</td><td>7.0%</td></tr><tr><td>5</td><td>25</td><td>5.8%</td><td>0.5%</td><td>28.4%</td><td>28.4%</td></tr><tr><td>6</td><td>36</td><td>6.6%</td><td>0.5%</td><td>13.0%</td><td>19.2%</td></tr><tr><td>7</td><td>49</td><td>9.1%</td><td>0.5%</td><td>11.7%</td><td>31.3%</td></tr><tr><td rowspan="3">4</td><td rowspan="3">5</td><td rowspan="3">10</td><td>3</td><td>27</td><td>8.0%</td><td>0.0%</td><td>2.3%</td><td>8.0%</td></tr><tr><td>4</td><td>64</td><td>6.8%</td><td>0.0%</td><td>2.0%</td><td>2.0%</td></tr><tr><td>5</td><td>125</td><td>13.5%</td><td>0.1%</td><td>5.0%</td><td>5.0%</td></tr><tr><td rowspan="3">4</td><td rowspan="3">10</td><td rowspan="3">10</td><td>3</td><td>27</td><td>11.9%</td><td>0.1%</td><td>17.3%</td><td>17.3%</td></tr><tr><td>4</td><td>64</td><td>13.3%</td><td>0.0%</td><td>10.4%</td><td>21.0%</td></tr><tr><td>5</td><td>125</td><td>13.1%</td><td>0.4%</td><td>9.4%</td><td>10.0%</td></tr><tr><td rowspan="3">4</td><td rowspan="3">20</td><td rowspan="3">20</td><td>3</td><td>27</td><td>9.0%</td><td>0.9%</td><td>8.5%</td><td>18.1%</td></tr><tr><td>4</td><td>64</td><td>11.9%</td><td>0.5%</td><td>9.7%</td><td>21.7%</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr><td>5</td><td>5</td><td>10</td><td>3</td><td>81</td><td>14.7%</td><td>1.0%</td><td>4.4%</td><td>6.3%</td></tr><tr><td rowspan="2">5</td><td rowspan="2">10</td><td rowspan="2">10</td><td>4</td><td>256</td><td>13.9%</td><td>0.5%</td><td>11.2%</td><td>19.1%</td></tr><tr><td></td><td></td><td></td><td></td><td></td><td></td></tr></tbody></table>

$\nu^{LDR}$ bound gets noticeably worse. In contrast, $\nu^{2S-LDR}$ stays very close to the optimal value, with an average relative distance of 0.2% among all the instances (compared to 8.7% for $\nu^{LDR}$). We do not observe the same trend among the dual bounds, and their relative distance to the optimal value fluctuates even between two instances that only differ in the number of scenarios. Nevertheless, $\bar{\nu}_R^{NA}$ considerably outperforms the gap from PI, with an average improvement of 39.4%.

### 4.3.3. Budgeted-Uncertainty Instances
The $\nu^{2S-LDR}$ bound for the location-transportation instances of this section is obtained using the C&CG method of Section 2.2 and detailed models are given in Appendix A.1. Similar to the scenario-tree instances, we use the
Lastly, for the primal bounds, we observe that the solution times for obtaining the traditional and two-stage LDR bounds are comparable. For the dual bounds, the times for obtaining the PI bound are quite small whereas for our dual problem the solution time significantly increases with problem size with several instances not being solved to optimality within the 1-hour time limit.

## 5.3. Robust Location-Transportation Problem

The two-stage robust location-transportation problem studied by Zeng and Zhao (2013) is as follows. Given a set of $I$ potential facilities with building cost $f_i$ and unit capacity cost $a_i, i \in [I]$, we have to meet the uncertain demand of a set of customers $J$ with unit transportation cost $c_{ij}, i \in [I], j \in [J]$. The goal is to decide which facilities to open and their initial capacities, such that the worst-case total cost of facility deployment and future transportation is minimized. Letting $d_{jt}(\xi^t)$ be the demand of customer $j$ at stage $t$, we define the MSARO extension of the problem:

$$ \min z \qquad (33a) $$

$$ \text{s.t.} \quad z \ge \sum_{i \in [I]} (f_i y_i + a_i s_{i1}) + \sum_{t \in [2,T]} \sum_{i \in [I]} \left( a_i s_{it}(\xi^t) + \sum_{j \in [J]} c_{ij} x_{ijt}(\xi^t) \right) \quad \xi^T \in \Xi^T \qquad (33b) $$

$$ s_{i1} \le K_i y_i \qquad i \in [I] \qquad (33c) $$

$$ s_{it}(\xi^t) = s_{i1} - \sum_{j \in [J]} x_{ijt}(\xi^t) \qquad i \in [I], t = 2, \xi^t \in \Xi^t \qquad (33d) $$

$$ s_{it}(\xi^t) = s_{i,t-1}(\xi^{t-1}) - \sum_{j \in [J]} x_{ijt}(\xi^t) \qquad i \in [I], t \in [3, T], \xi^t \in \Xi^t \qquad (33e) $$

$$ \sum_{i \in [I]} x_{ijt}(\xi^t) \ge d_{jt}(\xi^t) \qquad j \in [J], t \in [2, T], \xi^t \in \Xi^t \qquad (33f) $$

$$ y \in \{0, 1\}^I, \quad s_1 \in \mathbb{R}_+^I \qquad (33g) $$

$$ s_t(\xi^t) \in \mathbb{R}_+^I, \quad x_t(\xi^t) \in \mathbb{R}_+^{I \times J} \qquad t \in [2, T], \xi^t \in \Xi^t, \qquad (33h) $$

where $y_i$ is a binary variable equal to 1 if facility $i$ is built, $s_{i1}$ determines the initial capacity of facility $i$, $s_{it}(\xi^t)$ is the state variable calculating the remaining capacity of facility $i$ at stage $t$, and $x_{ijt}(\xi^t)$ is the amount of goods transported from facility $i$ to customer $j$ at stage $t$. Objective function (33a) together with constraint (33b) measures the worst-case cost. Constraints (33c) bound the initial capacity of the facilities, whereas constraints (33d)-(33e) are the state equations calculating the remaining capacities. Constraints (33f) ensure that the customer demands are met.

### 5.3.1. Problem Instances
Our instances are generated using the instance parameters described in (Zeng and Zhao 2013) for the two-stage problem. For number of stages $T \in [3, 5]$, we have five combinations for the number of facilities and customers, $(I, J) \in \{(5, 5), (5, 7), (5, 10), (10, 10), (20, 20)\}$. Fixed installation, unit capacity, and unit transportation