Stochastic Decomposition Method for Two-Stage Distributionally Robust Optimization

Harsha Gangammanavar and Manish Bansal Department of Engineering Management, Information, and Systems, Southern Methodist University, Dallas, TX 75275 ([email protected]).Department of Industrial and Systems Engineering, Virginia Tech, Blacksburg, VA 24061 ([email protected]).

(First Submitted: November 4, 2020)

Abstract. In this paper, we present a sequential sampling-based algorithm for the two-stage distributionally robust linear programming (2-DRLP) models. The 2-DRLP models are defined over a general class of ambiguity sets with discrete or continuous probability distributions. The algorithm is a distributionally robust version of the well-known stochastic decomposition algorithm of Higle and Sen (Math. of OR 16(3), 650-669, 1991) for a two-stage stochastic linear program. We refer to the algorithm as the distributionally robust stochastic decomposition (DRSD) method. The key features of the algorithm include (1) it works with data-driven approximations of ambiguity sets that are constructed using samples of increasing size and (2) efficient construction of approximations of the worst-case expectation function that solves only two second-stage subproblems in every iteration. We identify conditions under which the ambiguity set approximations converge to the true ambiguity sets and show that the DRSD method asymptotically identifies an optimal solution, with probability one. We also computationally evaluate the performance of the DRSD method for solving distributionally robust versions of instances considered in stochastic programming literature. The numerical results corroborate the analytical behavior of the DRSD method and illustrate the computational advantage over an external sampling-based decomposition approach (distributionally robust L-shaped method).

1 Introduction

Many applications of practical interest have been formulated as stochastic programming (SP) models. The models with recourse, particularly in a two-stage setting, have gained wide acceptance across application domains. These two-stage stochastic linear programs (2-SLPs) can be stated as follows:

\displaystyle\min~{}\{c^{\top}x+\mathbb{E}_{P^{\star}}[Q(x,\tilde{\omega})]~{}|~{}x\in\mathcal{X}\}.

(1)

One needs to have complete knowledge of the probability distribution $P^{\star}$ to formulate the above problem. Alternatively, one must have an appropriate means to simulate observations of the random variable so that a sample average approximation (SAA) problem with a finite number of scenarios can be formulated and solved. In many practical applications, distribution associated with random parameters in the optimization model is not precisely known. It either has to be estimated from data or constructed by expert judgments, which tend to be subjective. In any case, identifying a distribution using available information may be cumbersome at best. Stochastic min-max programming that has gained significant attention in recent years under the name of distributionally robust optimization (DRO) is intended to alleviate the ambiguity is distributional information.

In this paper, we study a particular manifestation of the DRO problem in the two-stage setting, viz., the two-stage distributionally robust linear programming (2-DRLP) problem. This problem is stated as:

\displaystyle\min~{}\{f(x)=c^{\top}x+\mathbb{Q}(x)~{}|~{}x\in\mathcal{X}\}.

(2)

Here, $c$ is the coefficient vector of a linear cost function and $\mathcal{X}\subseteq\mathbb{R}^{d_{x}}$ is the feasible set of the first-stage decision vector. The feasible region $\mathcal{X}$ takes the form of a compact polyhedron $\mathcal{X}=\{x~{}|~{}Ax\geq b,x\geq 0\}$ , where $A\in\mathbb{R}^{m_{1}\times d_{x}}$ and $b\in\mathbb{R}^{m_{1}}$ . The function $\mathbb{Q}(x)$ is the worst-case expected recourse cost, which is formally defined as follows:

\displaystyle\mathbb{Q}(x)=\max_{P\in\mathfrak{P}}~{}\bigg{\{}\mathcal{Q}(x;P):=\mathbb{E}_{P}[Q(x,\tilde{\omega})]\bigg{\}}.

(3)

The random vector $\tilde{\omega}\in\mathbb{R}^{d_{\omega}}$ is defined on a measurable space $(\Omega,\cal F)$ , where $\Omega$ is the sample space equipped with the sigma-algebra $\mathcal{F}$ and $\mathfrak{P}$ is a set of continuous or discrete probability distributions defined on the measurable space $(\Omega,\mathcal{F})$ . The set of probability distributions is often referred to as the ambiguity set. The expectation operation $\mathbb{E}_{P}[\cdot]$ is taken with respect to the probability distribution $P\in\mathfrak{P}$ . For a given $x\in\mathcal{X}$ , we refer to the optimization problem in (3) as the distribution separation problem. For a given realization $\omega$ of the random vector $\tilde{\omega}$ , the recourse cost in (3) is the optimal value of the following second-stage linear program:

	$\displaystyle Q(x,\omega):=\min\quad$	$\displaystyle g(\omega)^{\top}y$		(4)
	s.t.	$\displaystyle y\in\mathcal{Y}(x,\omega):=\big{\{}W(\omega)y=r(\omega)-T(\omega)x,~{}y\geq 0\big{\}}\subset\mathbb{R}^{d_{y}}.$

The second-stage parameters $g\in\mathbb{R}^{d_{y}}$ , $W\in\mathbb{R}^{m\times d_{y}}$ , $r\in\mathbb{R}^{m}$ , and $T\in\mathbb{R}^{m\times d_{x}}$ can depend on uncertainty.

Most data-driven approaches for 2-SLP, such as SAA, tackle the problem in two steps. In the first step, an uncertainty representation is generated using a finite set of observations that serves as an approximation of $\Omega$ and the corresponding empirical distribution serves as an approximation of $P^{\star}$ . For a given uncertainty representation, one obtains a deterministic approximation of (1). In the second step, the approximate problem is solved using deterministic optimization methods. Such a two-step approach may lead to poor out-of-sample performance, forcing the entire process to be repeated from scratch with an improved uncertainty representation. Since sampling is performed prior to the optimization step, the two-step approach is also referred to as the external sampling procedure.

The data-driven approaches for DRO problems avoid working with a fixed approximation of $P^{\star}$ in the first step. However, the ambiguity set is still defined either over the original sample space $\Omega$ or a finite approximation of it. Therefore, the resulting ambiguity set is a deterministic representation of the true ambiguity set in (2). Once again, deterministic optimization tools are employed to solve the DRO problem. In many data-driven settings, prior knowledge of the sample space may not be available, and using a finite sample to approximate the original sample space may result in similar out-of-sample performance as in the case of the external sampling approach for 2-SLP.

1.1 Contributions

In light of the above observations regarding the two-step procedure for data-driven optimization, the main contributions of this manuscript are highlighted in the following.

1.

A Sequential Sampling Algorithm: We present a sequential sampling approach for solving 2-DRLP. We refer to this algorithm as the distributionally robust stochastic decomposition (DRSD) algorithm following its risk-neutral predecessor, the two-stage stochastic decomposition (SD) method [20]. The algorithm uses a sequence of ambiguity sets that evolve over the course of the algorithm due to the sequential inclusion of new observations. While the simulation step improves the representation of the ambiguity set, the optimization step improves the solution in an online manner. Therefore, the DRSD method concurrently performs simulation and optimization steps. Moreover, the algorithm design does not depend on any specific ambiguity set description, and hence, is suitable for a general family of ambiguity sets.
2.

Convergence Analysis: The DRSD method is an inexact bundle method that creates outer linearization for the dynamically evolving approximation of the first-stage problem. We provide the asymptotic analysis of DRSD and identify conditions on ambiguity sets under which the sequential sampling approach identifies an optimal solution to the 2-DRLP problem in (2) with probability one.
3.

Computational Evidence of Performance: We provide the first evidence that illustrates the advantages of a sequential sampling approach for DRO through computational experiments conducted on well-known instances in SP literature.

1.2 Related work

The DRSD method has its roots in two-stage SP, in particular two-stage SD. In this subsection, we review the related two-stage SP literature along with decomposition and reformulation-based approaches for 2-DRLP.

For 2-SLP problems with finite support, including the SAA problem, the L-shaped method due to Van Slyke and Wets [40] has proven to be very effective. Other algorithms for 2-SLP’s such as the Dantzig-Wolfe decomposition [12] and the progressive hedging (PH) algorithm [29] also operate on problems with finite support. The well-established theory of SAA (see Chapter 5 in [37]) supports the external sampling procedure for 2-SLP. The quality of the solution obtained by solving an SAA problem is assessed using the procedures developed, e.g., in [4]. When the quality of the SAA solution is not acceptable, a new SAA is constructed with a larger number of observations. Prior work, such as [5] and [31], provide rules on how to choose the sequence of sample sizes in a sequential SAA procedure.

In contrast to the above, SD-based methods incorporate one new observation in every iteration to create approximations of the dynamically updated SAA of (1). First proposed in [20], this method has seen significant development in the past three decades with the introduction of quadratic regularization term in [21], statistical optimality rules [23], and extensions to multistage stochastic linear programs [18, 34]. The DRSD method presented in this manuscript extends the sequential sampling approach (i.e., SD) for 2-SLPs to DRO problems. Since the simulation of new observations and optimization steps are carried out in every iteration of the SD-based methods, they can also be viewed as internal sampling methods.

The concept of DRO dates back to the work of Scarf [33], and has gained significant attention in recent years. We refer the reader to [27] for a comprehensive treatment on various aspects of the DRO. The algorithmic works on DRO are either decomposition-based or reformulation-based approaches. The decomposition-based methods for 2-DRLP mimic the two-stage SP approach of using a deterministic representation of the sample space using a finite number of observations. As a consequence, the SP solution methods with suitable adaptation can be applied to solve the 2-DRLP problems. For instance, Breton and El Hachem [10, 11] apply the PH algorithm for a 2-DRLP model with a moment-based ambiguity set. Riis and Anderson [28] extend the L-shaped method for 2-DRLP with continuous recourse and moment-based ambiguity set. Bansal et.al. [1] extend the algorithm in [28], which they refer to as the distributionally robust (DR) L-shaped method, to solve 2-DRLPs, with ambiguity set defined over a polytope, in finite iterations. Further extensions of this decomposition approach are presented in [1] and [2] for DRO with mixed-binary recourse and disjunctive programs, respectively.

Another predominant approach to solve 2-DRLP problems is to reformulate the distribution separation problem in (3) as a minimization problem, pose the problem in (2) as a single deterministic optimization problem, and use off-the-shelf deterministic optimization tools to solve the reformulation. For example, Shapiro and Kleywegt [38] and Shapiro and Ahmed [35] provided approaches for 2-DRLP problem with moment matching set to derive an equivalent stochastic program with a certain reference distribution. Bertsimas et al. [7] provided tight semidefinite programming reformulations for 2-DRLP where the ambiguity set is defined using multivariate distributions with known first and second moments. Likewise, Hanasusanto and Kuhn [19] provided a conic programming reformulation for 2-DRLP problem where the ambiguity set comprises of a 2-Wasserstein ball centered at a discrete distribution. Xie [41] provided similar reformulations to tractable convex programs for 2-DRLP problems with ambiguity set defined using $\infty-$ Wasserstein metric. Jiang and Guan [24] reduced the worst-case expectation in 2-DRLP, where the ambiguity set is defined using $l_{1}$ -norm on the space of all (continuous and discrete) probability distributions, to a convex combination of CVaR and an essential supremum. By taking the dual of inner maximization, Love and Bayraksan [3] demonstrated that 2-DRLP where the ambiguity set is defined using $\phi$ -divergence and finite sample space is equivalent to 2-SLP with a coherent risk measure. When reformulations result in equivalent stochastic programs (in [24, 3, 36], for instance), a SAA of the reformulation is used to obtain an approximate solution.

Data-driven approaches for DRO have been presented for specific ambiguity sets. In [13], problems with ellipsoidal moment-based ambiguity set whose parameters are estimated using sampled data are addressed. Esfahani et. al. tackled data-driven problems with Wasserstein metric-based ambiguity sets with convex reformulations in [26]. In both these works, the authors provide finite-sample performance guarantees that probabilistically bound the gap between approximate and true DRO problems. Sun and Xu presented asymptotic convergence analysis of DRO problems with ambiguity sets that are based on moments and mixture distributions constructed using a finite set of observations in [39]. A practical approach to incorporate the results of these works to identify a high-quality DRO solution will be similar to the sequential SAA procedure for SP in [5]. Such an approach will involve the following steps performed in a series – a deterministic representation of ambiguity set using sampled observations, applying appropriate reformulation, and solving the resulting deterministic optimization problem. If the quality of the solution is deemed insufficient, then the entire series of steps is repeated with an improved representation of the ambiguity set (possibly with a larger number of observations).

Organization

The remainder of the paper is organized as follows. In §2, we present the two key ideas of the DRSD, viz., the sequential approximation of the ambiguity set and the recourse function. We provide a detailed description of the DRSD method in §3. We show the convergence of the value functions and solutions generated by the DRSD method in §4. We present results of our computational experiments in §5, and finally we conclude and discuss about potential extensions of this paper in §6.

Notations and Definitions

We define the ambiguity sets over $\mathcal{M}$ , the set of all finite signed measures on the measurable space $(\Omega,\mathcal{F})$ . A nonnegative measure (written as $P\succeq 0$ ) that satisfies $P(\Omega)=1$ is a probability measure. For probability distributions $P,P^{\prime}\in\mathfrak{P}$ , we define

\displaystyle\text{dist}(P,P^{\prime}):=\sup_{F\in\mathcal{F}}\Big{|}\mathbb{E}_{P}[F(\tilde{\omega})]-\mathbb{E}_{P^{\prime}}[F(\tilde{\omega})]\Big{|}

(5)

as the uniform distance of expectation, where $\mathcal{F}$ is a class of measurable functions. The above is the distance with $\zeta$ -structure that is used for the stability analysis in SP [30]. The distance between a single probability distribution $P$ to a set of distributions $\mathfrak{P}$ is given as $\text{dist}(P,\mathfrak{P})=\inf_{P^{\prime}\in\mathfrak{P}}d(P,P^{\prime})$ . The distance between two sets of probability distributions $\mathfrak{P}$ and $\widehat{\mathfrak{P}}$ is given as

\displaystyle\mathbb{D}(\mathfrak{P},\widehat{\mathfrak{P}}):=\sup_{P\in\widehat{\mathfrak{P}}}\text{dist}(P,\mathfrak{P}).

Finally, the Hausdorff distance between $\mathfrak{P}$ and $\widehat{\mathfrak{P}}$ is defined as

\displaystyle\mathbb{H}(\mathfrak{P},\widehat{\mathfrak{P}}):=\max\{\mathbb{D}(\mathfrak{P},\widehat{\mathfrak{P}}),~{}\mathbb{D}(\widehat{\mathfrak{P}},\mathfrak{P})\}.

With suitable definitions for the set $\mathcal{F}$ , the distance in (5) accepts the bounded Lipschitz, the Kantorovich and the $p$ -th order Fourier-Mourier metrics (see [30]).

2 Approximating Ambigiuty Set and Recourse Function

In this section, we present the building blocks for the DRSD method. Specifically, we present the procedures to obtain approximations of the ambiguity set $\mathfrak{P}$ and the recourse function $Q(x,\omega)$ . These procedures will be embedded within a sequential sampling-based approach. Going forward we make the following assumptions on the 2-DRLP models:

(A1)

The first-stage feasible region $\mathcal{X}$ is a non-empty and compact set.
(A2)

The recourse function satisfies relatively complete recourse. The dual feasible region of the recourse problem is nonempty compact polyhedral set. The transfer (or technology) matrix satisfies $\sup_{P\in\mathfrak{P}}\mathbb{E}_{P}[T(\tilde{\omega})]<\infty$ .
(A3)

The randomness affects the right-hand sides of constraints in (4).
(A4)

Sample space $\Omega$ is a compact metric space and the ambiguity set $\mathfrak{P}\neq\emptyset$ .

As a consequence of (A2), the recourse function satisfies $Q(x,\tilde{\omega})<\infty$ with probability one for all $x\in\mathcal{X}$ . It also implies that the second-stage feasible region, i.e., $\{y~{}:~{}Wy=r(\omega)-T(\omega)x,~{}y\geq 0\}$ , is non-empty for all $x\in\mathcal{X}$ and every $\omega\in\Omega$ . The non-empty dual feasible region $\Pi$ implies that there exists a $L>-\infty$ such that $Q(x,\tilde{\omega})>L$ . Without loss of generality, we assume that $Q(x,\tilde{\omega})\geq 0$ . As a consequence of (A3), the cost coefficient vector $g$ and the recourse matrix $W$ are not affected by uncertainty. Problems that satisfy (A4) are said to have a fixed recourse. Finally, the compactness of the support $\Omega$ guarantees that every probability measure $P\in\mathfrak{P}$ is tight.

2.1 Approximating the Ambiguity Set

The DRO approach assumes only partial knowledge about the underlying uncertainty that is captured by a suitable description of the ambiguity set. An ambiguity set must capture the true distribution with absolute or high degree of certainty, and must be computationally manageable. In this section we present a family of of ambiguity sets that are of interest to us in this paper.

The computational aspects of solving a DRO problem relies heavily on the structure of the ambiguity set. The description of these structures involve parameters which are determined based on practitioner’s risk-preferences. The ambiguity set descriptions that are prevalent in the literature include moment-based ambiguity sets with linear constraints (e.g., [14]) or conic constraints (e.g., [13]); Kantorovich distance or Wasserstein metric-based ambiguity sets [25]; $\zeta$ -structure metrics [42], $\phi$ -divergences such as $\chi^{2}$ distance and Kullback-Leibler divergence [6]; Prokhorov metrics [15], among others. Although the design of the DRSD method can work with any ambiguity set description defined over a compact sample space, we use 2-DRLPs with moment-based ambiguity sets and Wasserstein distance-based ambiguity sets to illustrate the algorithm in details.

In a data-driven setting, the parameters used in the description of ambiguity sets are estimated using a finite set of independent observations which can either be past realizations of the random variable $\tilde{\omega}$ or simulated by an oracle. We will denote such a sample by $\Omega^{k}\subseteq\Omega$ with $\Omega^{k}=\{\omega^{j}\}_{j=1}^{k}$ . Naturally, we can view $\Omega^{k}$ as a random sample and define the empirical frequency

\displaystyle\hat{p}^{k}(\omega^{j})=\frac{\kappa(\omega^{j})}{k}\qquad\text{for all }\omega^{j}\in\Omega^{k},

(6)

where $\kappa(\omega^{j})$ denotes the number of times $\omega^{j}$ is observed in the sample. Since in sequential sampling setting, the sample set is updated within the optimization algorithm, it is worthwhile to note that the empirical frequency can be updated using the following recursive equations:

\displaystyle\hat{p}^{k}(\omega)=\left\{\begin{array}[]{ll}\theta^{k}\hat{p}^{k-1}(\omega)&\text{if}~{}\omega\in\Omega^{k-1},\omega\neq\omega^{k}\\ \theta^{k}\hat{p}^{k-1}(\omega)+(1-\theta^{k})&\text{if}~{}\omega\in\Omega^{k-1},\omega=\omega^{k}\\ (1-\theta^{k})&\text{if}~{}\omega\notin\Omega^{k-1},\omega=\omega^{k}.\end{array}\right.

(10)

where $\theta^{k}\in[0,1]$ . We will succinctly denote the the above operation using the mapping $\Theta^{k}:\mathbb{R}^{|\Omega^{k-1}|}\rightarrow\mathbb{R}^{|\Omega^{k}|}$ .

In remainder of this section, we will present alternate descriptions of ambiguity sets and show the construction of what we will refer to as approximate ambiguity sets, denoted by $\mathfrak{P}^{k}$ . Let $\mathcal{F}^{k}=\sigma(\omega^{j}~{}|~{}j\leq k)$ be the $\sigma$ -algebra generated by the observations in the sample $\Omega^{k}$ . Notice that $\mathcal{F}^{k-1}\subseteq\mathcal{F}^{k}$ , and hence, $\{\mathcal{F}^{k}\}_{k\geq 1}$ is a filtration. We will define the approximate ambiguity sets over the measurable space $(\Omega^{k},\mathcal{F}^{k})$ . These sets should be interpreted to include all distributions that could have generated using the sample $\Omega^{k}$ , which share a certain relationship with sample statistics. We will use $\mathcal{M}^{k}$ to denote the finite signed measures on $(\Omega^{k},\mathcal{F}^{k})$ .

2.1.1 Moment-based Ambiguity Sets

Given the first $q$ moments associated with the random variable $\tilde{\omega}$ , the moment-based ambiguity set can be defined as

\displaystyle\mathfrak{P}_{\text{mom}}=\left\{P\in\mathcal{M}\left|\begin{array}[]{l}\int_{\Omega}dP(\tilde{\omega})=1,\\ \int_{\Omega}\psi_{i}(\tilde{\omega})dP(\tilde{\omega})=b_{i}\qquad i=1,\ldots,q\end{array}\right.\right\}.

(13)

While the first constraint ensures the definition of a probability measure, the moment requirements are guaranteed by the second constraints. Here, $\psi_{i}(\tilde{\omega})$ denote real valued measurable function on $(\Omega,\mathcal{F})$ and $b_{i}\in\mathbb{R}$ be a scalar for $i=1,\ldots,q$ . Existence of moments ensures that $b_{i}<\infty$ for all $i=1,\ldots,q$ . Notice that the description of the ambiguity set requires explicit knowledge of the following statistics: support $\Omega$ and the moments $b_{i}$ for $i=1,\ldots,q$ . In the data-driven setting, the support is approximated by $\Omega^{k}$ and the sample moments $\hat{b}_{i}^{k}=(1/k)\sum_{j=1}^{k}\psi_{i}(\omega^{j})$ are used to define the following approximate ambiguity set

\displaystyle\widehat{\mathfrak{P}}_{\textup{mom}}^{k}=\left\{P\in\mathcal{M}^{k}\left|\begin{array}[]{l}\sum_{\omega\in\Omega^{k}}p(\omega)=1,\\ \sum_{\omega\in\Omega^{k}}p(\omega)\psi_{i}(\omega)=\hat{b}_{i}^{k}\qquad i=1,\ldots,q\end{array}\right.\right\}.

(16)

The following result characterizes the relationship between distributions drawn from the above approximate ambiguity set, as well as asymptotic behavior of the sequence $\{\widehat{\mathfrak{P}}_{\textup{mom}}^{k}\}_{k\geq 1}$ .

Proposition 2.1.

For any $P\in\widehat{\mathfrak{P}}_{\textup{mom}}^{k-1}$ , we have $\Theta^{k}(P)\in\widehat{\mathfrak{P}}_{\textup{mom}}^{k}$ where $\theta^{k}=\frac{k-1}{k}$ . Further, suppose $\widehat{\mathfrak{P}}_{\textup{mom}}^{k}\neq\emptyset$ for all $k\geq 0$ , $\mathbb{H}(\widehat{\mathfrak{P}}_{\textup{mom}}^{k},\mathfrak{P}_{\textup{mom}})\rightarrow 0$ as $k\rightarrow\infty$ , almost surely.

Proof.

See Appendix §A. ∎

In the context of DRO, similar ambiguity sets have been studied in [8, 14] where only the first moment (i.e., $q=1$ ) was considered. The form of ambiguity set above also relates to those used in [13, 28, 33, 39] where constraints were imposed only on the mean and covariance. In the data driven setting of [13] and [39], the statistical estimates are used in constructing the approximate ambiguity set as in the case of (16). However, the ambiguity sets in these previous works are defined over the original sample space $\Omega$ , as opposed to $\Omega^{k}$ that is used in (16). This marks a critical deviation in the way the approximate ambiguity set are constructed.

Remark 2.1.

When the moment information is available about the underlying distribution $P^{\star}$ , an approximate moment-based ambiguity set with constant parameters in (16) (i.e., with $\hat{b}_{i}^{k}=b_{i}$ for all $k$ ) can be constructed. Such an approximate ambiguity sets defined over $\Omega^{k}$ is studied in [28]. Notice that these approximate ambiguity sets satisfy $\cup_{k\geq 1}\widehat{\mathfrak{P}}^{k}\subseteq\mathfrak{P}$ and $\widehat{\mathfrak{P}}^{k}\subseteq\widehat{\mathfrak{P}}^{k+1}$ , for all $k\geq 1$ .

2.1.2 Wasserstein distance-based Ambiguity Sets

We next present approximations of another class of ambiguity sets that has gained significant attention in the DRO literature, viz., the Wasserstein distance-based ambiguity sets. Consider probability distributions $\mu_{1},\mu_{2}\in\mathcal{M}$ , and a function $\nu:\Omega\times\Omega\rightarrow\mathbb{R}_{+}\cup\{\infty\}$ such that $\nu$ is symmetric, $\nu^{\frac{1}{r}}(\cdot)$ satisfies triangle inequality for $1\leq r<\infty$ , and $\nu(\omega_{1},\omega_{2})=0$ whenever $\omega_{1}=\omega_{2}$ . If $\Pi(\mu_{1},\mu_{2})$ denotes the joint distribution of random vectors $\omega_{1}$ and $\omega_{2}$ with marginals $\mu_{1}$ and $\mu_{2}$ , respectively, then the Wasserstein metric of order $r$ is given by

\displaystyle d_{\text{w}}(\mu_{1},\mu_{2})=\inf_{\eta\in\Pi(\mu_{1},\mu_{2})}\bigg{\{}\int_{\Omega\times\Omega}\nu(\omega_{1},\omega_{2})\eta(d\omega_{1},d\omega_{2})\bigg{\}}.

(17)

In the above definition, the decision variable $\eta\in\Pi$ can be viewed as a plan to transport goods/mass from an entity whose spatial distribution is given by the measure $\mu_{1}$ to another entity with spatial distribution $\mu_{2}$ . Therefore, the $d_{\text{w}}(\mu_{1},\mu_{2})$ measures the optimal transport cost between the measures. Notice that an arbitrary norm $\|\bullet\|^{r}$ on $\mathbb{R}^{d_{\omega}}$ satisfies the requirement of the function $c(\cdot)$ . In this paper, we will use $\|\bullet\|$ as the choice of our metric, in which case we obtain the Wasserstein metric of order 1. Using this metric, we define an ambiguity set as follows:

\displaystyle\mathfrak{P}_{\textup{w}}=\{P\in\mathcal{M}~{}|~{}d_{W}(P,P^{*})\leq\epsilon\}

(18)

for a given $\epsilon>0$ and a reference distribution $P^{*}$ . As was done in §2.1.1, we define approximate ambiguity sets defined over the measurable space $(\Omega^{k},\mathcal{F}^{k})$ as follows:

\displaystyle\widehat{\mathfrak{P}}_{\textup{w}}^{k}=\{P\in\mathcal{M}^{k}~{}|~{}d_{W}(P,\widehat{P}^{k})\leq\epsilon\}.

(19)

For the above approximate ambiguity set, the distribution separation problem in (3) takes the following form:


$\displaystyle\max~{}$	$\displaystyle\sum_{\omega\in\Omega^{k}}p(\omega)Q(x,\omega)$	(20a)
subject to	$\displaystyle P\in\widehat{\mathfrak{P}}^{k}_{\textup{w}}=\left\{P\in\mathcal{M}^{k}\left\|\begin{array}[]{l}\sum_{\omega\in\Omega^{k}}p(\omega)=1\\ \sum_{\omega^{\prime}\in\Omega^{k}}\eta(\omega,\omega^{\prime})=p(\omega)\qquad\forall\omega\in\Omega^{k},\\ \sum_{\omega\in\Omega^{k}}\eta(\omega,\omega^{\prime})=\hat{p}^{k}(\omega^{\prime})\qquad\forall\omega^{\prime}\in\Omega^{k},\\ \sum_{(\omega,\omega^{\prime})\in\Omega^{k}\times\Omega^{k}}\\|\omega-\omega^{\prime}\\|\eta(\omega,\omega^{\prime})\leq\epsilon\\ \eta(\omega,\omega^{\prime})\geq 0\quad\forall\omega,\omega^{\prime}\in\Omega^{k}\end{array}\right.\right\}.$	(20g)

The following result characterizes the distributions drawn from the approximate ambiguity sets of the form in (19), or equivalently (20g).

Proposition 2.2.

The sequence of Wasserstein distance-based approximate ambiguity set satisfies the following properties (1) for any $P\in\widehat{\mathfrak{P}}_{\textup{w}}^{k-1}$ , we have $\Theta^{k}(P)\in\widehat{\mathfrak{P}}_{\textup{w}}^{k}$ where $\theta^{k}=\frac{k-1}{k}$ , and (2) $\mathbb{H}(\widehat{\mathfrak{P}}_{\textup{w}}^{k},\mathfrak{P}_{\textup{w}})\rightarrow 0$ as $k\rightarrow\infty$ , almost surely.

Proof.

See appendix §A. ∎

The approximate ambiguity set in [26] is a ball constructed in the space of probability distributions that are defined over the sample space $\Omega$ and whose radius reduces with increase in the number of observations. Using Wasserstein balls of shrinking radii, the authors of [26] show that the optimal value of the sequence of DRO problems converges to the optimal value of the expectation-valued SP problem in (1) associated with the true distribution $P^{\star}$ . The approximate ambiguity set in (19), on the other hand, uses a constant radius for all $k\geq 1$ . In this regard, we consider settings where the ambiguity is not necessarily resolved with increasing number of observations. This is reflected in the approximate ambiguity sets (16) and (19) that converge to their respective true ambiguity sets (13) and (18), respectively.

2.2 Approximating the Recourse Problem

Cutting plane methods for the 2-SLPs use an outer linearization-based approximation of the first-stage objective function in (1). In such algorithms, the challenging aspect of computing the expectation is addressed by taking advantage of the structure of the recourse problem (4). Specifically, for a given $\omega$ , the recourse value $Q(\cdot,\omega)$ is known to be convex in the right-hand side parameters that includes the first-stage decision vector $x$ . Additionally, if the support of $\tilde{\omega}$ is finite and (A2) holds, then the function $Q(\cdot,\omega)$ is polyhedral. Under assumptions (A2) and (A4), these structural property of convexity extends to the expected recourse value $\mathcal{Q}(x)$ .

Due to strong duality of linear programs, the recourse value is also equal to the optimal value of the dual of (4), i.e.,

	$\displaystyle Q(x,\omega)=\max~{}$	$\displaystyle\pi^{\top}[r(\omega)-T(\omega)x]$		(21)
	subject to	$\displaystyle\pi\in\Pi:=\{\pi~{}\|~{}W^{\top}\pi\leq g\}.$

Due to (A2) and (A4), the dual feasible region $\Pi$ is a polytope that is not affected by the uncertainty. If $\Pi\subset\Pi$ denotes the set of all extreme points of the polytope $\Pi$ , then the recourse value can also be expressed as the pointwise maximum of affine functions computed using elements of the set $\Pi$ as given below.

\displaystyle Q(x,\omega)=\max_{\pi\in\Pi}\pi^{\top}[r(\omega)-T(\omega)x].

(22)

The outer linearization approaches tend to approximate the above form of recourse function by identifying the extreme points (optimal solutions to (21)) at a sequence of candidate (or trial) solutions $\{x^{k}\}$ , and generating the corresponding affine functions. If $\pi(x^{k},\omega)$ is the optimal dual obtained by solving (21) with $x^{k}$ as input, then the affine function $\alpha^{k}(\omega)+(\beta^{k}(\omega))^{\top}x$ is obtained by computing the coefficients $\alpha^{k}(\omega)=(\pi^{k}(\omega))^{\top}r(\omega)$ and $\beta^{k}(\omega)=C(\omega)^{\top}\pi^{k}(\omega)$ . Following linear programming duality, notice that this affine function is a supporting hyperplane to $Q(x,\omega)$ at $x^{k}$ , and lower bounds the function at every other $x\in\mathcal{X}$ .

If the support $\Omega$ is finite, then one can solve a dual subproblem for all $\omega\in\Omega$ with the candidate solution as input, generate the affine functions, and collate them together to obtain an approximation of the first-stage objective function. This is the essence of the L-shaped method applied to 2-SLP in (1). In each iteration of the L-shaped method, the affine functions generated using a candidate solution $x^{k}$ and information gathered from individual observations are weighed by the probability density of the observation to update the approximate first-stage objective function. Notice that even when SAA of the first-stage objective function of the 2-SLP, using a sample $\Omega_{N}\subset\Omega$ of size $N$ , the L-shaped method can be applied. A similar approximation strategy is used in the DR L-shaped method for 2-DRLP problems.

Alternatively, we can consider the following approximation of the recourse function expressed in the form given in (22):

\displaystyle Q^{k}(x,\omega)=\max_{\pi\in\Pi^{k}}\pi^{\top}[r(\omega)-C(\omega)x].

(23)

Notice that the above approximation is built using only a subset $\Pi^{k}\subset\Pi$ of extreme points, and therefore, satisfies $Q^{k}(x,\omega)\leq Q(x,\omega)$ . Since $Q(x,\omega)\geq 0$ , we begin with $\Pi^{0}=\{0\}$ . Subsequently, we construct a sequence of sets $\{\Pi^{k}\}$ such that $\Pi^{k}\subseteq\Pi^{k+1}\subseteq\ldots\subset\Pi$ that ensures $Q^{k}(x,\omega)\geq 0$ for all $k$ . The following result from [20] captures the behavior of the sequence of approximation $\{Q^{k}\}$ .

Proposition 2.3.

The sequence $\{Q^{k}(x,\omega)\}_{k\geq 1}$ converges uniformly to a continuous function on $\mathcal{X}$ for any $\omega\in\Omega$ .

Proof.

See Appendix A. ∎

The approximations of the form in (23) is one of the principal features of the SD algorithm (see [20, 21]). While the L-shaped and DR L-shaped methods require a finite support for $\tilde{\omega}$ , SD is applicable even for problems with continuous support. The algorithm uses an “incremental” SAA for the first-stage objective function by adding one new observation in each iteration. Therefore, the first-stage objective function approximation used in SD is built using the recourse problem approximation in (23) and the incremental SAA. This approximation is given by:

\displaystyle\mathcal{Q}^{k}(x)=c^{\top}x+\frac{1}{k}\sum_{j=1}^{k}Q^{k}(x,\omega^{j}).

(24)

The affine functions generated in SD provide an outer linearization for the approximation in (24). The sequence of sets that grow monotonically in size, viz. $\{\Pi^{k}\}$ , is generated by adding one new vertex to the previous set $\Pi^{k-1}$ to obtain the updated set $\Pi^{k}$ . The newly added vertex is the optimal dual solution obtained by solving (21) with the most recent observation $\omega^{k}$ and candidate solution $x^{k}$ as input.

We refer the reader to [9], [1, 28], and [20, 22] for the a detailed exposition of the L-shaped, the DR L-Shaped, and the SD methods, respectively, and note only the key differences between these methods. Firstly, the sample used to in the (DR) L-shaped method is fixed prior to optimization. In SD, this sample is updated dynamically throughout the course of the algorithm. Secondly, exact subproblem optimization for all observations in the sample is used in every iteration of the (DR) L-shaped method. On the contrary, exact optimization is used only for the subproblems corresponding to the latest observation, and an “argmax” procedure (to be described in the next section) is used for observations encountered in earlier iterations.

3 Distributionally Robust Stochastic Decomposition

In this paper, we focus on a setting where the ambiguity set $\mathfrak{P}$ is approximated by a sequence of ambiguity sets $\{\widehat{\mathfrak{P}}^{k}\}_{k>0}$ such that the following properties are satisfied: ( $i$ ) for any $P\in\widehat{\mathfrak{P}}^{k-1}$ , there exists $\theta^{k}\in[0,1]$ such that $\Theta^{k}(P)\in\widehat{\mathfrak{P}}^{k}$ and ( $ii$ ) $\mathbb{H}(\widehat{\mathfrak{P}}^{k},\mathfrak{P})\rightarrow 0$ as $k\rightarrow\infty$ , almost surely. The moment-based ambiguity set $\mathfrak{P}_{\text{mom}}$ and Wasserstein distance-based ambiguity set $\mathfrak{P}_{w}$ are two sets that satisfy these properties (Propositions 2.1 and 2.2, respectively). Recall that the approximate ambiguity set in an iteration (say $k$ ) is constructed using a finite set of observations $\Omega^{k}$ that progressively grow in size. Note that the sequence of approximate ambiguity sets may not necessarily converge to a single distribution. In other words, we do not assume that increasing sample size will overcome ambiguity asymptotically, as is the case in [26, 42].

Algorithm 1 Distributionally Robust Stochastic Decomposition

1:Input: Incumbent solution

\hat{x}^{1}\in\mathcal{X}

; initial sample

\Omega^{0}\subseteq\Omega

; stopping tolerance

\tau>0

;

\gamma\in(0,1]

, and minimum iterations

k^{\min}

2:Initialization: Set iteration counter

k\leftarrow 1

;

\Pi^{0}=\emptyset

;

\mathcal{L}^{0}=\emptyset

, and

f^{0}(x)=0

3:while (

k\geq k^{\min}

and

f^{k-1}(\hat{x}^{k-1})-f^{k-1}(x^{k})\geq\tau f^{k-1}(\hat{x}^{k-1})

) do

4: Solve the master problem

\mathcal{M}^{k}

(25) to obtain a candidate solution

x^{k}

5: Generate a scenario

\omega^{k}\in\Omega

to get sample

\Omega^{k}\leftarrow\Omega^{k-1}\cup\{\omega^{k}\}

6: for

\omega\in\Omega^{k}

7: if

\omega=\omega^{k}

then

8: Solve the second-stage linear program (4) with (

x^{k},\omega

) as input;

9: Obtain the optimal value

Q(x^{k},\omega)

and optimal dual solution

\pi(x^{k},\omega)

;

10: Update dual vertex set

\Pi^{k}\leftarrow\Pi^{k-1}\cup\{\pi(x^{k},\omega)\}

11: else

12: Use the argmax procedure (26) to identify dual vertex

\pi(x^{k},\omega^{k})

;

13: Store

Q^{k}(x^{k},\omega)=(\pi(x^{k},\omega))^{\top}[r(\omega)-T(\omega)x^{k}]

14: end if

15: end for

16: Solve the distribution separation problem using the ambiguity set

\widehat{\mathfrak{P}}^{k}

and

\{Q^{k}(x^{k},\omega)\}_{\omega\in\Omega^{k}}

to get an extremal distribution

P^{k}:=(p^{k}(\omega))_{\omega\in\Omega_{k}}

17: Derive affine function

{\ell}_{k}^{k}(x)={\alpha}_{k}^{k}+({\beta}_{k}^{k})^{\top}x

using

\{\pi(x^{k},\omega)\}_{\omega\in\Omega^{k}}

and

P^{k}

to get lower bound approximation of

\mathbb{Q}^{k}(x)

as in (29);

18: Perform Steps 6-17 with

\hat{x}^{k-1}

(incumbent solution) to obtain

\hat{\ell}_{k}^{k}(\cdot)

19: for

\ell_{j}^{k-1}\in\mathcal{L}^{k-1}

20: Update previously generated affine functions

\ell^{k-1}_{j}(x)={\alpha}_{k-1}^{j}+({\beta}_{k-1}^{j})^{\top}x

\alpha^{k}_{j}=\theta^{k}\alpha^{k-1}_{j}\text{ and }\beta^{k}_{j}=\theta^{k}\alpha^{k-1}_{j};

21: Set

{\ell}^{k}_{j}(x)={\alpha}^{k}_{j}+({\beta}^{k}_{j})^{\top}x

that provides lower bound approx. of

\mathbb{Q}^{k}(x)

;

22: end for

23: Build a collection of these affine functions, denoted by

\mathcal{L}^{k}

;

24: Update approximation of the first-stage objective function:

\displaystyle c^{\top}x+\mathbb{Q}^{k}(x)\geq f^{k}(x)=c^{\top}x+\max_{j\in\mathcal{L}^{k}}~{}\{\alpha^{k}_{j}+(\beta^{k}_{i})^{\top}x\};

25: if Incumbent update rule (33) is satisfied then

26:

\hat{x}^{k}\leftarrow\hat{x}^{k-1}

and

i_{k}\leftarrow i_{k-1}

;

27: else

28:

\hat{x}^{k}\leftarrow x^{k}

and

i_{k}\leftarrow k

;

29: end if

30: Update the master problem by replacing

f^{k-1}(x)

with

f^{k}(x)

to obtain

\mathcal{M}^{k+1}

;

31:

k\leftarrow k+1

;

32:end while

33:return

x^{k}

The pseudocode of the DRSD method is given in Algorithm 1. We present the main steps of the DRSD method in iteration $k$ (Steps 4-31 of Algorithm 1). At the beginning of iteration $k$ , we have a certain approximation of the first-stage objective function that we denote as $f^{k-1}(x)$ , a finite set of observations $\Omega^{k-1}$ and an incumbent solution $\hat{x}^{k-1}$ . We will use the term incumbent solution to refer to the best solution discovered by the algorithm until iteration $k$ . The solution identified in the current iteration will be referred to as the candidate solution and denote it as $x^{k}$ (without $\hat{\bullet}$ ).

Iteration $k$ begins by first identifying the candidate solution by solving the following the master problem (Step 4):

\displaystyle x^{k}\in\arg\min~{}\{f^{k-1}(x)~{}|~{}x\in\mathcal{X}\},

(25)

denoted by $\mathcal{M}^{k}$ . Following this, a new observation $\omega^{k}\in\Omega$ is realized, and added to the current sample of observations $\Omega^{k-1}$ to get $\Omega^{k}=\Omega^{k-1}\cup\{\omega^{k}\}$ (Step 5).

In order to build the first-stage objective function approximation, we rely upon the recourse function approximation presented in Section 2.2. For the most recent observation $\omega^{k}$ and the candidate solution $x^{k}$ , we evaluate the recourse function value $Q(x^{k},\omega^{k})$ by solving (4), and obtain the dual optimum solution $\pi(x^{k},\omega^{k})$ . Likewise, we obtain dual optimum solution $\pi(\hat{x}^{k-1},\omega^{k})$ by solving (4) for incumbent solution $\hat{x}^{k-1}$ (Steps 7–10). These dual vectors are added to a set $\Pi^{k-1}$ of previously discovered optimal dual vectors. In other words, we recursively update $\Pi^{k}\leftarrow\Pi^{k-1}\cup\{\pi(x^{k},\omega^{k}),\pi(\hat{x}^{k-1},\omega^{k})\}$ . For all other observations ( $\omega\in\Omega^{k},\ \omega\neq\omega^{k}$ ), we identify a dual vector in $\Pi^{k}$ that provides the best lower bounding approximation at $\{Q(x^{k},\omega)\}$ using the following operation (Steps 11–13):

\pi(x^{k},\omega)\in\arg\max~{}\{\pi^{\top}[r(\omega)-T(\omega)x^{k}]~{}|~{}\pi\in\Pi^{k}\}.

(26)

Note that the calculations in (26) are carried out only for previous observations as $\pi(x^{k},\omega^{k})$ provides the best lower bound at $Q(x^{k},\omega^{k})$ . Further, notice that

\pi(x^{k},\omega)^{\top}[r(\omega)-T(\omega)x^{k}]=Q^{k}(x^{k},\omega),

the approximate recourse function value at $x^{k}$ defined in (23), for all $\omega\in\Omega^{k}$ , and $Q^{k}(x^{k},\omega^{k})=Q(x^{k},\omega^{k})$ .

Using $\{Q^{k}(x^{k},\omega^{j})\}_{j=1}^{k}$ , we solve a distribution separation problem (in Step 16):

\displaystyle\mathbb{Q}^{k}(x^{k})=\max~{}\bigg{\{}\sum_{\omega\in\Omega^{k}}p(\omega)Q^{k}(x^{k},\omega)~{}|~{}p(\omega)\in\widehat{\mathfrak{P}}^{k}\bigg{\}}.

(27)

Let $P^{k}=(p^{k}(\omega))_{\omega\in\Omega^{k}}$ denote the optimal solution of the above problem which we identify as the maximal/extremal probability distribution. Since the problem is solved over measures $\mathcal{M}^{k}$ that are defined only over the observed set $\Omega^{k}$ , the maximal probability distribution has weights $p^{k}(\omega^{j})$ for $\omega^{j}\in\Omega^{k}$ , and $p^{k}(\omega)=0$ for $\omega\in\Omega\setminus\Omega^{k}$ . Notice that the problem in (3) differs from the distribution separation problem (27) as the latter uses the recourse function approximation $Q^{k}(\cdot)$ and approximate ambiguity set $\widehat{\mathfrak{P}}^{k}$ as opposed to the true recourse function $Q(\cdot)$ and ambiguity set $\mathfrak{P}$ , respectively. For the ambiguity sets in (2.1) the distribution separation problem is a deterministic linear program. In general, the distribution separation problems associated with well-known ambiguity sets remain deterministic convex optimization problems [27], and off-the-shelf solvers can used to obtain the extremal distribution.

In Step 17 of Algorithm 1, we use the dual vectors $\{\pi(x^{k},\omega^{j})\}_{j\leq k}$ and the maximal probability distribution $P^{k}$ to generate a lower bounding affine function:

\displaystyle\mathbb{Q}^{k}(x)=\max_{P\in\widehat{\mathfrak{P}}^{k}}\mathbb{E}_{P}[Q^{k}(x,\tilde{\omega})]\geq\sum_{\omega^{j}\in\Omega^{k}}p^{k}(\omega^{j})\cdot(\pi(x^{k},\omega^{j}))^{\top}[r(\omega^{j})-T(\omega^{j})x],

(28)

for the worst case expected recourse function measured with respect to the maximal probability distribution $P^{k}\in\widehat{\mathfrak{P}}^{k}$ which is obtained by solving the distribution separation problem (27). We denote the coefficients of the affine function on the right-hand side of (28) by

\displaystyle\alpha_{k}^{k}=\sum_{\omega^{j}\in\Omega^{k}}p^{k}(\omega^{j})\pi(x^{k},\omega^{j})^{\top}r(\omega^{j})\text{ and }\beta_{k}^{k}=-\sum_{\omega^{j}\in\Omega^{k}}p^{k}(\omega^{j})T(\omega^{j})^{\top}\pi(x^{k},\omega^{j}),

(29)

and succinctly write the affine function as $\ell_{k}^{k}(x)=\alpha_{k}^{k}+(\beta_{k}^{k})^{\top}x$ . Similar calculations are carried out using the incumbent solution $\hat{x}$ to identify a maximal probability distribution and a lower bounding affine function resulting in the affine function $\hat{\ell}_{k}^{k}(x)=\hat{\alpha}_{k}^{k}+(\hat{\beta}_{k}^{k})^{\top}x$ .

While the latest affine functions provide lower bound on $\mathbb{Q}^{k}$ , the affine functions generated in previous iteration are not guaranteed to lower bound $\mathbb{Q}^{k}$ . To see this, let us consider the moment-based approximate ambiguity sets $\{\widehat{\mathfrak{P}}^{k}_{\text{mom}}\}_{k\geq 1}$ . Let $P^{j}_{\text{mom}}\in\widehat{\mathfrak{P}}^{j}_{\text{mom}}$ be the maximal distribution identified in an iteration $j<k$ which was used to compute the affine function $\ell_{j}^{j}(x)$ . By assigning $p^{j}(\omega)=0$ for all new observations encountered after iteration $j$ , i.e., $\omega\in\Omega^{k}\setminus\Omega^{j}$ , we can construct a probability distribution $\bar{P}=((p^{j}(\omega))_{\omega\in\Omega^{j}},(0)_{\omega\in\Omega^{k}\setminus\Omega^{j}})\in\mathbb{R}^{|\Omega^{k}|}_{+}$ . This reconstructed distribution satisfies $\sum_{\omega\in\Omega^{k}}\bar{p}(\omega)=1$ . However, it is easy to see that it does not satisfy $\sum_{\omega\in\Omega^{k}}\psi_{i}(\omega)\bar{p}(\omega)=\hat{b}^{j}_{i}=\hat{b}^{k}_{i}$ for all $i=1,\ldots,q$ . Therefore, $\bar{P}\notin\mathfrak{P}^{k}$ . In other words, while the coefficients $(\alpha_{j}^{j},\beta_{j}^{j})$ are $\mathcal{F}^{j}$ -measurable, the corresponding measure is not feasible to the approximate ambiguity set $\mathfrak{P}^{k}$ . Therefore, $\ell_{j}^{j}(x)$ is not a valid lower bound to $\mathbb{Q}^{k}$ . The arguments for the Wasserstein-based approximate ambiguity set are more involved, but persistence of a similar issue can be demonstrated.

In order to address this, we update the previously generated affine functions $\ell^{k-1}_{j}(x)={\alpha}_{k-1}^{j}+({\beta}_{k-1}^{j})^{\top}x$ for $j<k$ , as follows (Steps 19 - 22):

\displaystyle\alpha_{j}^{k}=\theta^{k}\alpha_{j}^{k-1},\ \ \beta_{j}^{k}=\theta^{k}\beta_{j}^{k-1},\text{ and }{\ell}^{k}_{j}(x)={\alpha}_{k}^{j}+({\beta}_{k}^{j})^{\top}x\qquad\text{ for all }j<k,

(30)

such that $\ell^{k}_{j}(x)$ provides lower bound approximation of $\mathbb{Q}^{k}(x)$ for all $j\in\{1,\ldots,k-1\}$ . Similarly, we update the affine functions $\hat{\ell}^{k}_{j}(x)$ , $j<k$ , associated with incumbent solution (Step 18). Recall that for $\widehat{\mathfrak{P}}_{\text{mom}}$ and $\widehat{\mathfrak{P}}_{w}$ , the parameter $\theta^{k}=\frac{k-1}{k}$ (Propositions 2.1 and 2.2). The candidate and the incumbent affine functions ( $\ell_{k}^{k}(x)$ and $\hat{\ell}_{k}^{k}(x)$ , respectively), as well as the updated collection of previously generated affine functions are used to build the set of affine functions which we will denote by $\mathcal{L}^{k}$ (Step 23). The lower bounding property of this first-stage objective function approximation is captured in the following result.

Theorem 3.1.

Under assumption (A2), the first-stage objective function approximation in (32) satisfies

\displaystyle f^{k}(x)\leq c^{\top}x+\mathbb{Q}^{k}(x)\text{ for all }x\in\mathcal{X}\text{ and }k\geq 1.

Proof.

For non-empty approximate ambiguity set $\widehat{\mathfrak{P}}^{1}$ of ambiguity set $\mathfrak{P}$ , the construction of the affine function ensures that $\ell_{1}^{1}(x)\leq\mathbb{Q}^{1}(x)$ . Assume that $\ell(x)\leq\mathbb{Q}^{k-1}(x)$ for all $\ell\in\mathcal{L}^{k-1}$ and $k>1$ . The maximal nature of the probability distribution $P^{k}$ satisfies:

\displaystyle\sum_{\omega\in\Omega^{k}}p^{k}(\omega)Q^{k}(x,\omega)

\displaystyle\geq\sum_{\omega\in\Omega^{k}}p(\omega)Q^{k}(x,\omega)\qquad\forall P\in\widehat{\mathfrak{P}}^{k}.

Using above and the monotone property of the approximate recourse function, we have

	$\displaystyle\sum_{\omega\in\Omega^{k}}p^{k}(\omega)Q^{k}(x,\omega)$	$\displaystyle\geq\sum_{\omega\in\Omega^{k}}p(\omega)Q^{k-1}(x,\omega)$
		$\displaystyle=\sum_{\omega\in\Omega^{k}\setminus\omega^{k}}p(\omega)Q^{k-1}(x,\omega)+p(\omega^{k})Q^{k-1}(x,\omega^{k}),$		(31)

for all $\{p(\omega)\}_{\omega\in\Omega^{k}}\in\widehat{\mathfrak{P}}^{k}$ . Based on the properties of $\mathfrak{P}$ and $\{\widehat{\mathfrak{P}}^{k}\}_{k\geq 1}$ (similar to Propositions 2.1 and 2.2), we know that for every $P\in\widehat{\mathfrak{P}}^{k-1}$ we can construct a probability distribution in $\mathfrak{P}^{k}$ using the mapping $\Theta^{k}$ defined by (10). Considering a probability distribution $P^{\prime}=\{p^{\prime}(\omega)\}_{\omega\in\Omega^{k-1}}\in\widehat{\mathfrak{P}}^{k-1}$ such that $\Theta^{k}(P^{\prime})\in\widehat{\mathfrak{P}}^{k}$ , the inequality (31) reduces to

	$\displaystyle\sum_{\omega\in\Omega^{k}}p^{k}(\omega)Q^{k}(x,\omega)$	$\displaystyle\geq\sum_{\omega\in\Omega^{k}\setminus\omega^{k}}[\theta^{k}p^{\prime}(\omega)Q^{k-1}(x,\omega)]+[\theta^{k}p^{\prime}(\omega^{k})+(1-\theta^{k})]Q^{k-1}(x,\omega^{k})$
		$\displaystyle=\theta^{k}\bigg{[}\sum_{\omega\in\Omega^{k-1}}p^{\prime}(\omega)Q^{k-1}(x,\omega)\bigg{]}+(1-\theta^{k})Q^{k-1}(x,\omega^{k})$
		$\displaystyle\geq\theta^{k}\bigg{[}\sum_{\omega\in\Omega^{k-1}}p^{\prime}(\omega)Q^{k-1}(x,\omega)\bigg{]}.$

The last inequality is due to assumption (A2), i.e., $Q(x,\omega^{k})\geq 0$ and the construction of recourse function approximation $Q^{k}$ described in §2.2. Since $\ell(x)$ lower bounds the term in bracket, we have

\displaystyle\sum_{\omega\in\Omega^{k}}p^{k}(\omega)Q^{k}(x,\omega)

\displaystyle\geq\theta^{k}\ell(x).

Using the same arguments for all $\ell\in\mathcal{L}^{k-1}$ , and the fact that the $\ell_{k}^{k}(x)$ and $\hat{\ell}(x)$ are constructed as lower bounds to the $\mathbb{Q}^{k}$ , we have $f^{k}(x)\leq c^{\top}x+\mathbb{Q}^{k}(x)$ . This completes the proof by induction. ∎

Using the collection of affine functions $\mathcal{L}^{k}$ , we update approximation of the first-stage objective function in Step 24, as follows:

\displaystyle f^{k}(x)=c^{\top}x+\max_{i\in\mathcal{L}^{k}}~{}\{\alpha^{i}+(\beta^{i})^{\top}x\}.

(32)

Once the approximation is updated, the performance of the candidate solution is compared relative to the incumbent solution (Steps 25–32). This comparison is performed by verifying if the following inequality

\displaystyle f^{k}(x^{k})-f^{k}(\hat{x}^{k-1})<\gamma[f^{k-1}(x^{k})-f^{k-1}(\hat{x}^{k-1}],

(33)

where parameter $\gamma\in(0,1]$ , is satisfied. If so, the candidate solution is designated to be the next incumbent solution, i.e., $\hat{x}^{k}=x^{k}$ . If the inequality is not satisfied, the precious incumbent solution is retained as $\hat{x}^{k}=\hat{x}^{k-1}$ . This completes an iteration of the DRSD method.

Remark 3.1.

The algorithm design can be extended to incorporate 2-DRLP where the relatively complete recourse assumption of (A2) and/or assumption (A3) is not satisfied. For problems where relatively complete recourse condition is not met, a candidate solution may lead to one or more subproblems to be infeasible. In this case, the dual extreme rays can be used to compute a feasibility cut that is included in the first-stage approximation. The argmax procedure in (26) is only valid when assumption (A3) is satisfied. In problems where the uncertainty also affects the cost coefficients, the argmax procedure presented in [17] can be utilized. These algorithmic enhancements can be incorporated without affecting the convergence properties of DRSD that we present in the next section.

4 Convergence Analysis

In this section we provide the convergence result of the sequential sampling-based approach to solve DRO problems. In order to facilitate the exposition of our theoretical results, we will define certain quantities for notational convenience that are not necessarily computed during the course of the algorithm in the form presented in the previous section. Our convergence results are built upon stability analyses presented in [39] and convergence analysis of the SD algorithm in [20].

We define a function which is defined over the approximate ambiguity set using the recourse function $Q(\cdot,\cdot)$ , that is

\displaystyle g^{k}(x):=c^{\top}x+\max_{P\in\widehat{\mathfrak{P}}^{k}}\mathbb{E}_{P}[Q(x,\tilde{\omega})]

(34)

for a fixed $x\in\mathcal{X}$ . We begin by analyzing the behavior of the sequence $\{g^{k}\}_{k\geq 1}$ as $k\rightarrow\infty$ . In particular, we will assess the sequence of function evaluations at a converging subsequence of first-stage solutions. The result is captured in the following proposition.

Proposition 4.1.

Suppose $\{\hat{x}^{k_{n}}\}$ denotes a subsequence of $\{\hat{x}^{k}\}$ such that $\hat{x}^{k_{n}}\rightarrow\bar{x}$ , then $\lim_{n\rightarrow\infty}|g^{k_{n}}(\hat{x}^{k_{n}})-f(\bar{x})|=0$ , with probability one.

Proof.

Consider the ambiguity set $\widehat{\mathfrak{P}}^{k}$ . For $i=1,2$ and $x_{i}\in\mathcal{X}$ , let $P(x_{i})\in\mathrm{argmax}_{P\in\widehat{\mathfrak{P}}^{k}}\{\mathbb{E}_{P}[Q(x_{i},\tilde{\omega})]$ . Then,

	$\displaystyle g^{k}(x_{1})=~{}$	$\displaystyle c^{\top}x_{1}+\mathbb{E}_{P(x_{1})}[Q(x_{1},\tilde{\omega})]$
	$\displaystyle\geq~{}$	$\displaystyle c^{\top}x_{1}+\mathbb{E}_{P(x_{2})}[Q(x_{1},\tilde{\omega})]$
	$\displaystyle=~{}$	$\displaystyle c^{\top}x_{2}+\mathbb{E}_{P(x_{2})}[Q(x_{2},\tilde{\omega}]+c^{\top}(x_{1}-x_{2})+$
		$\displaystyle\hskip 85.35826pt\mathbb{E}_{P(x_{2})}[Q(x_{1},\tilde{\omega})]-\mathbb{E}_{P(x_{2})}[Q(x_{2},\tilde{\omega})]$
	$\displaystyle=~{}$	$\displaystyle g^{k}(x_{2})+c^{\top}(x_{1}-x_{2})+\mathbb{E}_{P(x_{2})}[Q(x_{1},\tilde{\omega})]-\mathbb{E}_{P(x_{2})}[Q(x_{2},\tilde{\omega})].$

The inequality in the above follows from optimality of $P(x_{1})$ . The above implies that

$\displaystyle g^{k}(x_{2})-g^{k}(x_{1})\leq~{}$	$\displaystyle c^{\top}(x_{2}-x_{1})+\mathbb{E}_{P(x_{2})}[Q(x_{2},\tilde{\omega})]-\mathbb{E}_{P(x_{2})}[Q(x_{1},\tilde{\omega})]$
$\displaystyle\leq~{}$	$\displaystyle\|c^{\top}(x_{2}-x_{1})\|+\bigg{\|}\mathbb{E}_{P(x_{2})}[Q(x_{2},\tilde{\omega})]-\mathbb{E}_{P(x_{2})}[Q(x_{1},\tilde{\omega})]\bigg{\|}$
$\displaystyle\leq~{}$	$\displaystyle(\\|c\\|+C)\\|x_{2}-x_{1}\\|.$	(35)

The second relationship is due to the triangular inequality. The third inequality follows from uniform Lipschitz continuity of recourse function $Q(x,\tilde{\omega})$ , with probability one, under assumption (A2) which implies that there exists a constant $C$ such that $|\mathbb{E}_{P}[Q(x_{1},\tilde{\omega})]-\mathbb{E}_{P}[Q(x_{2},\tilde{\omega})]|\leq C\|x_{1}-x_{2}\|$ for any $P$ . Therefore, the function $g^{k}(x)$ is equi-continuous on $x\in\mathcal{X}$ . Starting with $x_{2}$ and using the same arguments, we have

\displaystyle g^{k}(x_{1})-g^{k}(x_{2})\leq~{}

\displaystyle(\|c\|+C)\|x_{2}-x_{1}\|.

(36)

Now consider ambiguity sets $\mathfrak{P}$ and $\widehat{\mathfrak{P}}^{k}$ . Note that

	$\displaystyle\|f(x)-g^{k}(x)\|=$	$\displaystyle~{}\bigg{\|}\max_{P\in\mathfrak{P}}\mathbb{E}_{P}[Q(x,\tilde{\omega})]-\max_{P^{\prime}\in\widehat{\mathfrak{P}}^{k}}\mathbb{E}_{P^{\prime}}[Q(x,\tilde{\omega})]\bigg{\|}\qquad\forall x\in\mathcal{X}$
	$\displaystyle\leq$	$\displaystyle~{}\max_{P\in\mathfrak{P}}\min_{P^{\prime}\in\widehat{\mathfrak{P}}^{k}}\big{\|}\mathbb{E}_{P}[Q(x,\tilde{\omega})]-\mathbb{E}_{P^{\prime}}[Q(x,\tilde{\omega}]\big{\|}\qquad\forall x\in\mathcal{X}$
	$\displaystyle\leq$	$\displaystyle~{}\max_{P\in\mathfrak{P}}\min_{P^{\prime}\in\widehat{\mathfrak{P}}^{k}}\sup_{x\in\mathcal{X}}\big{\|}\mathbb{E}_{P}[Q(x,\tilde{\omega})]-\mathbb{E}_{P^{\prime}}[Q(x,\tilde{\omega})]\big{\|}.$

Using the definition of deviation between ambiguity sets $\mathfrak{P}$ and $\widehat{\mathfrak{P}}^{k}$ as well as its relationship with Hausdorff distance, we have

\displaystyle|f(x)-g^{k}(x)|=\mathbb{D}(\mathfrak{P},\widehat{\mathfrak{P}}^{k})\leq\mathbb{H}(\mathfrak{P},\widehat{\mathfrak{P}}^{k}).

(37)

For $\hat{x}^{k_{n}}$ and $\bar{x}$ , combining (35), (36), and (37), we have

	$\displaystyle\|f(\bar{x})-g^{k_{n}}(\hat{x}^{k_{n}})\|\leq$	$\displaystyle~{}\|f(\bar{x})-g^{k_{n}}(\bar{x})\|+\|g^{k_{n}}(\bar{x})-g^{k_{n}}(\hat{x}^{k_{n}})\|$
	$\displaystyle\leq$	$\displaystyle~{}\mathbb{H}(\mathfrak{P},\widehat{\mathfrak{P}}^{k})+(\\|c\\|+C)\\|\bar{x}-\hat{x}^{k_{n}}\\|.$

As $n\rightarrow\infty$ , the family of ambiguity sets considered satisfy $\mathbb{H}(\mathfrak{P},\widehat{\mathfrak{P}}^{k_{n}})\rightarrow 0$ and the hypothesis informs us that $\hat{x}^{k_{n}}\rightarrow\bar{x}$ . Therefore, we conclude that $g^{k_{n}}(\hat{x}^{k_{n}})\rightarrow f(\bar{x})$ as $n\rightarrow\infty$ . ∎

Notice that, the behavior of the approximate ambiguity sets defined in §2.1, in particular, the condition $\mathbb{H}(\mathfrak{P},\widehat{\mathfrak{P}}^{k})\rightarrow 0$ as $k\rightarrow\infty$ plays a central role in the above proof. Recall that for the moment and Wasserstein distance-based ambiguity sets, the condition is established in propositions 2.1 and 2.2, respectively. It is also worthwhile to note that under the foregoing conditions, (37) also implies uniform convergence of the sequence $\{g^{k}\}$ to $f(x)$ , with probability one.

The above result applies to any algorithm that generates a converging sequence of iterates $\{x^{k}\}$ and a corresponding sequence of extremal distributions. Such an algorithm is guaranteed to exhibit convergence to the optimal distributionally robust objective function value. Therefore, this result is applicable to the sequence of instances constructed using external sampling and solved, for example, using reformulation-based methods. Such an approach was adopted in [28] and [39]. The analysis in [28] relies upon two rather restrictive assumptions. The first assumption is that for all $P\in\mathfrak{P}$ there exists a sequence of measures $\{P^{k}\}$ such that $P^{k}\in\widehat{\mathfrak{P}}^{k}$ and converges weakly to $P$ . The second assumption requires the approximate ambiguity sets to be strict subsets of the true ambiguity set, i.e., $\widehat{\mathfrak{P}}^{k}\subset\mathfrak{P}$ . Both of these assumptions are very difficult to satisfy in a data-driven setting (also see Remark 2.1).

The analysis in [39], on the other hand, does not make the above assumptions. Therefore, their analysis is more broadly applicable in settings where external sampling is used to generate $\Omega^{k}$ . DRO instances are constructed based on statistics estimated using $\Omega^{k}$ and solved to optimality for each $k\geq 1$ . They show the convergence of optimal objective function values and optimal solution sets of approximate problems to the optimal objective function value and solutions of the true DRO problem, respectively. In this regard, the result in Proposition 4.1 can alternatively be derived using Theorem 1(i) in [39]. While the above function is not computed during the course of the sequential sampling algorithm, it provides the necessary benchmark for our convergence analysis.

One of the main point of deviation in our analysis stems from the fact that we use the objective function approximations that are built based on the approximate recourse function in (23). In order to study the piecewise affine approximation of the first-stage objective function, we introduce another benchmark function

\displaystyle\phi^{k}(x):=c^{\top}x+\max_{P\in\widehat{\mathfrak{P}}^{k}}\mathbb{E}_{P}[Q^{k}(x,\tilde{\omega})].

(38)

Notice that the above function uses the approximations for the ambiguity set (as in the case of (34)) as well as the approximation for recourse function. This construction ensures for all $x\in\mathcal{X}$ and $k\geq 1$ that $\phi^{k}(x)\leq g^{k}(x)$ that follows from the fact that $Q^{k}(x,\tilde{\omega})\leq Q(x,\tilde{\omega})$ , almost surely. Further, following the result in Theorem 3.1 ensures that $f^{k}(x)\leq\phi^{k}(x)$ . Putting these together, we obtain the following relationship

\displaystyle f^{k}(x)\leq\phi^{k}(x)\leq g^{k}(x)\qquad\forall x\in\mathcal{X},k\geq 1.

(39)

While the previous proposition was focused on the upper limit in the above relationship, we present the asymptotic behavior of the $\{f^{k}\}$ sequence in the following results.

Lemma 4.2.

Suppose $\{\hat{x}^{k_{n}}\}$ denotes a subsequence of $\{\hat{x}^{k}\}$ such that $\hat{x}^{k_{n}}\rightarrow\bar{x}$ , $\lim_{n\rightarrow\infty}f^{k_{n}}(\hat{x}^{k_{n}})-f(\bar{x})=0$ , with probability one.

Proof.

From Proposition 4.1, we have $\lim_{n\rightarrow\infty}|f(\bar{x})-g^{k_{n}}(\hat{x}^{k_{n}})|=0$ . Therefore, there exists $N_{1}<\infty$ and $\epsilon_{1}>0$ such that

\displaystyle\bigg{|}\max_{P\in\mathfrak{P}}\mathbb{E}_{P}[Q(\bar{x},\tilde{\omega})]-\max_{P\in\widehat{\mathfrak{P}}^{k_{n}}}\mathbb{E}_{P}[Q(\hat{x}^{k_{n}},\tilde{\omega})]\bigg{|}<\epsilon_{1}/2\quad\forall n>N_{1}.

(40)

Now consider,

	$\displaystyle\max_{P\in\widehat{\mathfrak{P}}^{k_{n}}}\mathbb{E}_{P}[Q(\hat{x}^{k_{n}},\tilde{\omega})]$	$\displaystyle-\max_{P\in\widehat{\mathfrak{P}}^{k_{n}}}\mathbb{E}_{P}[Q^{k_{n}}(\hat{x}^{k_{n}},\tilde{\omega})]$
	$\displaystyle=$	$\displaystyle\max_{P\in\widehat{\mathfrak{P}}^{k_{n}}}\big{(}\mathbb{E}_{P}[Q(\hat{x}^{k_{n}},\tilde{\omega})]-\mathbb{E}_{P}[Q^{k_{n}}(\hat{x}^{k_{n}},\tilde{\omega})]\big{)}$
	$\displaystyle\leq$	$\displaystyle~{}\max_{P\in\widehat{\mathfrak{P}}^{k_{n}}}\big{\|}\mathbb{E}_{P}[Q(\hat{x}^{k_{n}},\tilde{\omega})]-\mathbb{E}_{P}[Q^{k_{n}}(\hat{x}^{k_{n}},\tilde{\omega})]\big{\|}$
	$\displaystyle=$	$\displaystyle~{}\max_{P\in\widehat{\mathfrak{P}}^{k_{n}}}\mathbb{E}_{P}[\|Q(\hat{x}^{k_{n}},\tilde{\omega})-Q^{k_{n}}(\hat{x}^{k_{n}},\tilde{\omega})\|].$

The last equality follows from the fact that $Q(x,\tilde{\omega})\geq Q^{k}(x,\tilde{\omega})$ for all $x\in\mathcal{X}$ and $k\geq 1$ , almost surely. Moreover, because of the uniform convergence of $\{Q^{k}\}$ (Proposition 2.3), the sequence of approximate functions $\{\phi^{k}\}$ also convergences uniformly. This implies that, there exists $N_{2}<\infty$ such that

\displaystyle\bigg{|}\max_{P\in\widehat{\mathfrak{P}}^{k_{n}}}\mathbb{E}_{P}[Q(\hat{x}^{k_{n}},\tilde{\omega})]-\max_{P\in\widehat{\mathfrak{P}}^{k_{n}}}\mathbb{E}_{P}[Q^{k_{n}}(\hat{x}^{k_{n}},\tilde{\omega})]\bigg{|}<\epsilon_{1}/2.

(41)

Let $N=\max\{N_{1},N_{2}\}$ . Using (40) and (41), we have for all $n>N$

\displaystyle\bigg{|}\max_{P\in\mathfrak{P}}\mathbb{E}_{P}[Q(\bar{x},\tilde{\omega})]-\max_{P\in\widehat{\mathfrak{P}}^{k_{n}}}\mathbb{E}_{P}[Q^{k_{n}}(x^{k_{n}},\tilde{\omega})]\bigg{|}<\epsilon_{1}.

This implies that $|f(\bar{x})-\phi^{k_{n}}(\hat{x}^{k_{n}})|\rightarrow 0$ as $n\rightarrow\infty$ . Based on (26), we have $Q^{k_{n}}(\hat{x}^{k_{n}},\omega)=(\pi(\hat{x}^{k_{n}},\omega))^{\top}[r(\omega)-T(\omega)\hat{x}^{k_{n}}]\geq(\pi(\hat{x}^{k_{n}},\omega))^{\top}[r(\omega)-T(\omega)x]$ for all $x\in\mathcal{X}$ and $\omega\in\Omega^{k_{n}}$ . Let

\displaystyle\alpha^{k_{n}}_{k_{n}}=\sum_{\omega\in\Omega^{k_{n}}}p^{k_{n}}(\omega)(\pi(\hat{x}^{k_{n}},\omega))^{\top}r(\omega)\text{ and }\beta^{k_{n}}_{k_{n}}=-\sum_{\omega\in\Omega^{k_{n}}}p^{k_{n}}(\omega)T(\omega)^{\top}\pi(\hat{x}^{k_{n}},\omega),

where $\{p^{k_{n}}(\omega)\}_{\omega\in\Omega^{k_{n}}}$ is an optimal solution of the distributional separation problem (27) where index $k$ is replaced by $k_{n}$ . Then, the affine function $\alpha^{k_{n}}_{k_{n}}+(c+\beta^{k_{n}}_{k_{n}})^{\top}x$ provides a lower bound approximation for function $\phi^{k_{n}}(x)$ , i.e.,

\displaystyle\phi^{k_{n}}(x)\geq\alpha^{k_{n}}_{k_{n}}+(c+\beta^{k_{n}}_{k_{n}})^{\top}x\qquad\text{ for all }x\in\mathcal{X},

with strict equality holding only at $\hat{x}^{k_{n}}$ . Therefore, using the definition of $f^{k}(x)$ we have $\lim_{n\rightarrow\infty}\alpha^{k_{n}}_{k_{n}}+(c+\beta^{k_{n}}_{k_{n}})^{\top}\hat{x}^{k_{n}}=\lim_{n\rightarrow\infty}f^{k_{n}}(\hat{x}^{k_{n}})=\lim_{n\rightarrow\infty}\phi^{k_{n}}(\hat{x}^{k_{n}})=f(\bar{x})$ , almost surely. This completes the proof. ∎

The above result characterizes the behavior of the sequence of affine functions generated during the course of the algorithm. In particular, the sequence $\{f^{k}(\hat{x}^{k})\}_{k\geq 1}$ accumulates at the objective value of the original DRO problem (2). Recall that the candidate solution $x^{k}$ is a minimizer of $f^{k-1}(x)$ and an affine function is generated at this point such that $f^{k}(x^{k})=\phi^{k}(x^{k})$ in all iterations $k\geq 1$ . However, due to the update procedure in (30) the quality of the estimates at $x^{k}$ gradually diminishes leading to a large value for $(\phi^{k}(x^{k})-f^{k}(x^{k}))$ as $k$ increases. This emphasizes the role of the incumbent solution and computing the incumbent affine function $\hat{\ell}(x)$ during the course of the algorithm. By updating the incumbent solution and frequently reevaluating the affine functions at the incumbent solution, we can ensure that the approximation is “sufficiently good” in the neighborhood of the incumbent solution. In order to assess the improvement of approximation quality, we define

\displaystyle\delta^{k}:=f^{k-1}(x^{k})-f^{k-1}(\hat{x}^{k-1})\leq 0\qquad\forall k\geq 1.

(42)

The inequality follows from the optimality of $x^{k}$ with respect to the objective function $f^{k-1}$ . The quantity $\delta^{k}$ measures the error in objective function estimate at the candidate solution with respect to the estimate at the current incumbent solution. The following result captures the asymptotic behavior of this error term.

Lemma 4.3.

Let $\mathcal{K}$ denotes a sequence of iterations where the incumbent solution changes. There exists a subsequence of iterations, denoted as $\mathcal{K}^{*}\subseteq\mathcal{K}$ , such that $\lim_{k\in\mathcal{K}^{*}}\delta^{k}=0$ .

Proof.

We will consider two cases depending on whether the set $\mathcal{K}$ is finite or not. First, suppose that $|\mathcal{K}|$ is not finite. By the incumbent update rule and (42),

\displaystyle f^{k_{n}}(x^{k_{n}})-f^{k_{n}}(\hat{x}^{k_{n}-1})<\gamma[f^{k_{n}-1}(x^{k_{n}})-f^{{k_{n}}-1}(\hat{x}^{k_{n}-1})]=\gamma\delta^{k_{n}}\leq 0\qquad\forall k_{n}\in\mathcal{K}.

Subsequently, we have $\limsup_{n\rightarrow\infty}\delta^{k_{n}}\leq 0$ . Since $x^{k_{n}}=\hat{x}^{k_{n}}$ and $\hat{x}^{k_{n}-1}=\hat{x}^{k_{n-1}}$ , we have

\displaystyle f^{k_{n}}(\hat{x}^{k_{n}})-f^{k_{n}}(\hat{x}^{k_{n-1}})\leq\gamma\delta^{k_{n}}\leq 0.

The left-hand side of the above inequality captures the improvement in the objective function value at the current incumbent solution over the previous incumbent solution. Using the above, we can write the average improvement attained over $n$ incumbent changes as

\displaystyle\frac{1}{n}\sum_{j=1}^{n}\bigg{[}f^{k_{j}}(\hat{x}^{k_{j}})-f^{k_{j}}(\hat{x}^{k_{j-1}})\bigg{]}\leq\frac{1}{n}\sum_{j=1}^{n}\gamma\delta^{k_{j}}\leq 0\qquad\text{ for all }n.

This implies that

\displaystyle\frac{1}{n}\underbrace{\bigg{(}f^{k_{n}}(\hat{x}^{k_{n}})-f^{k_{1}}(\hat{x}^{k_{0}})\bigg{)}}_{(a)}+\frac{1}{n}\bigg{[}\sum_{j=1}^{n-1}\underbrace{\bigg{(}f^{k_{j}}(\hat{x}^{k_{j}})-f^{k_{j+1}}(\hat{x}^{k_{j}})\bigg{)}}_{(b)}\bigg{]}\leq\frac{1}{n}\sum_{j=1}^{n}\gamma\delta^{k_{j}}\leq 0,\ \ \forall n.

Under the assumption that the dual feasible region is non-empty and bounded (this is ensured by relatively complete recourse, (A2)), $\{f^{k}\}$ is a sequence of Lipschitz continuous functions. This along with compactness of $\mathcal{X}$ (A1), implies that $f^{k_{n}}(\hat{x}^{k_{n}})-f^{k_{1}}(\hat{x}^{k_{0}})$ is bounded from above. Hence, the term (a) reduces to zero as $n\rightarrow\infty$ . The term (b) converges to zero, with probability one, due to uniform convergence of $\{f^{k}\}$ . Since $\gamma\in(0,1]$ , we have

\displaystyle\lim_{n\rightarrow\infty}\frac{1}{n}\sum_{j=1}^{n}\delta^{k_{j}}=0

with probability one. Further,

\displaystyle\lim_{n\rightarrow\infty}\frac{1}{n}\sum_{j=1}^{n}\delta^{k_{j}}\leq\limsup_{n\rightarrow\infty}\delta^{k_{n}}\leq 0.

Thus, there exists a subsequence indexed by the set $\mathcal{K}^{*}$ such that $\lim_{k\in\mathcal{K}^{*}}\delta^{k}=0$ , with probability one.

Now if $|\mathcal{K}|$ is finite, then there exists $\hat{x}$ and $K<\infty$ such that for all $k\geq K$ , we have $\hat{x}^{k}=\hat{x}$ . Notice that, if $\lim_{k\in\mathcal{K}^{*}}x^{k}=\bar{x}$ , uniform convergence of the sequence $\{f^{k}\}$ and Lemma 4.2 ensure that


	$\displaystyle\lim_{k\in\mathcal{K}^{}}f^{k}(x^{k})=\lim_{k\in\mathcal{K}^{}}f^{k-1}(x^{k})=f(\bar{x})$		(43a)
	$\displaystyle\lim_{k\in\mathcal{K}^{}}f^{k}(\hat{x})=\lim_{k\in\mathcal{K}^{}}f^{k-1}(\hat{x})=f(\hat{x}).$		(43b)

Further, since the incumbent is not updated in iterations $k\geq K$ , we must have from the update rule in (33) that

\displaystyle f^{k}(x^{k})-f^{k}(\hat{x})\geq\gamma[f^{k-1}(x^{k})-f^{k-1}(\hat{x})]=\gamma\delta^{k}\quad\text{ for all }k\geq K.

Using (43), we have

\displaystyle\lim_{k\in\mathcal{K}^{*}}\big{(}f^{k}(x^{k})-f^{k}(\hat{x})\big{)}

\displaystyle\geq\gamma\lim_{k\in\mathcal{K}^{*}}\big{(}f^{k-1}(x^{k})-f^{k-1}(\hat{x})\big{)},

which implies

\displaystyle f(\bar{x})-f(\hat{x})

\displaystyle\geq\gamma(f(\bar{x})-f(\hat{x})).

Since $\gamma\in(0,1]$ , we must have $f(\bar{x})-f(\hat{x})=0$ . Hence, $\lim_{k\in\mathcal{K}^{*}}\delta^{k}=f(\bar{x})-f(\hat{x})=0$ , with probability one. ∎

Equipped with the results in lemmas 4.2 and 4.3, we state the main theorem which establishes the existence of a subsequence of incumbent solution sequence for which every accumulation point is an optimal solution to (2).

Theorem 4.4.

Let $\{x^{k}\}_{k=1}^{\infty}$ and $\{\hat{x}^{k}\}_{k=1}^{\infty}$ be the sequence candidate and incumbent solutions generated by the algorithm. There exists a subsequence $\{\hat{x}^{k}\}_{k\in\mathcal{K}}$ for which every accumulation point is an optimal solution of 2-DRLP (2), with probability one.

Proof.

Let $x^{*}\in\mathcal{X}$ be an optimal solution of (2). Consider a subsequence indexed by $\mathcal{K}$ such that $\lim_{k\in\mathcal{K}}\hat{x}^{k}=\bar{x}$ . Compactness of $\mathcal{X}$ ensures the existence of accumulation point $\bar{x}\in\mathcal{X}$ and therefore,

\displaystyle f(x^{*})\leq f(\bar{x}).

(44)

From Theorem 3.1, we have

	$\displaystyle f^{k}(x)\leq~{}$	$\displaystyle c^{\top}x+\mathbb{Q}^{k}(x)$
	$\displaystyle\leq~{}$	$\displaystyle c^{\top}x+\max_{P\in\widehat{\mathfrak{P}}^{k}}\mathbb{E}_{P}[Q(x,\tilde{\omega})]=g^{k}(x)\qquad\forall k,x\in\mathcal{X}.$

Thus, using the uniform convergence of $\{g^{k}\}$ (Proposition 4.1) we have

\displaystyle\limsup_{k\in\mathcal{K}^{\prime}}f^{k}(x^{*})\leq\lim_{k\in\mathcal{K}^{\prime}}g^{k}(x^{*})=f(x^{*})

(45)

for all subsequences indexed by $\mathcal{K}^{\prime}\subseteq\{1,2,\ldots\}$ , with probability one. Recall that,

\displaystyle\delta^{k}=f^{{k}-1}(x^{k})-f^{{k}-1}(\hat{x}^{{k}-1})

\displaystyle\leq f^{{k}-1}(x^{*})-f^{{k}-1}(\hat{x}^{{k}-1})\qquad\text{for all }k\geq 1.

The inequality in the above follows from the optimality of $x^{k}$ with respect to $f^{k-1}(x)$ . Taking limit over $\mathcal{K}$ , we have

	$\displaystyle\lim_{k\in\mathcal{K}}\delta^{k}\leq~{}$	$\displaystyle\lim_{k\in\mathcal{K}}\big{(}f^{k-1}(x^{*})-f^{{k}-1}(\hat{x}^{{k}-1})\big{)}$
	$\displaystyle\leq~{}$	$\displaystyle\limsup_{k\in\mathcal{K}}f^{k-1}(x^{*})-\liminf_{k\in\mathcal{K}}f^{{k}-1}(\hat{x}^{{k}-1})$
	$\displaystyle\leq~{}$	$\displaystyle f(x^{*})-f(\bar{x}).$

The last inequality follows from (45) and $\lim_{k\in\mathcal{K}}f^{k-1}(\hat{x}^{k-1})=f(\bar{x})$ (Lemma 4.2). From Lemma 4.3, there exists a subsequence indexed by $\mathcal{K}^{*}\subseteq\mathcal{K}$ such that $\lim_{k\in\mathcal{K}^{*}}\delta^{k}=0$ . Therefore, if $\{\hat{x}^{k}\}_{k\in\mathcal{K}^{*}}\rightarrow\bar{x}$ , we have

\displaystyle f(x^{*})-f(\bar{x})\geq 0.

Using (44) and the above inequality, we conclude that $\bar{x}$ is an optimal solution with probability one. ∎

5 Computational Experiment

In this section, we evaluate the effectiveness and efficiency of the DRSD method in solving 2-DRLP problems. For our preliminary experiments, we consider 2-DRLP problems with moment-based ambiguity set $\mathfrak{P}_{\text{mom}}$ for the first two moments ( $q=2$ ). We report results from the computational experiments conducted on four well-known SP test problems: CEP, PGP, BAA, and STORM. The test problems have supports of size $216$ , $576$ , $10^{18}$ and $10^{81}$ , respectively.

We use an external sampling-based approach as a benchmark for comparison. The external sampling-based instances involve constructing approximate problems of the form (34) with a pre-determined number of observations $N\in\{100,250,500,1000,2000\}$ (that might not be unique). The resulting instances are solved using the DR L-Shaped method. Specifically, we compare the solution quality provided by these methods along with the solution time.

We conduct $30$ independent replications for each problem instance with different seeds for the random number generator. The algorithms are implemented in the C programming language, and the experiments are conducted on a $64$ -bit Intel core i7 - 4770 CPU at $3.4$ GHz $\times~{}8$ machine with $32$ GB memory. All linear programs, i.e., master problem, subproblems, and distribution separation problem, are solved using CPLEX 12.10 callable subroutines. The results from the experiments are presented in Tables 1 and 2. The results for the instances with a finite sample size obtained from the DR L-shaped method are labeled as DRLS- $N$ , where $N$ denotes the number of observations used to approximate the ambiguity set. The DRSD method is run for a minimum of $256$ iteration, while no such limit is imposed on the DR L-shaped method.

Method	# Iterations	Objective Estimate	# Unique Obs.
PGP
DRLS-100	$17.867$ ( $\pm 0.92$ )	$457.610$ ( $\pm 3.28$ )	$38.233$ ( $\pm 1.10$ )
DRLS-250	$20.467$ ( $\pm 0.73$ )	$462.922$ ( $\pm 2.28$ )	$53.267$ ( $\pm 1.74$ )
DRLS-500	$20.467$ ( $\pm 0.63$ )	$464.704$ ( $\pm 1.95$ )	$68.667$ ( $\pm 1.62$ )
DRLS-1000	$20.500$ ( $\pm 0.84$ )	$466.104$ ( $\pm 1.78$ )	$85.833$ ( $\pm 1.78$ )
DRLS-2000	$20.900$ ( $\pm 0.57$ )	$469.579$ ( $\pm 2.53$ )	$102.867$ ( $\pm 1.87$ )
DRSD	$503.667$ ( $\pm 687.01$ )	$463.185$ ( $\pm 16.28$ )	$65.667$ ( $\pm 10.04$ )
CEP
DRLS-100	$2.600$ ( $\pm 0.19$ )	$658817.129$ ( $\pm 14457.30$ )	$80.700$ ( $\pm 1.21$ )
DRLS-250	$2.267$ ( $\pm 0.17$ )	$680735.606$ ( $\pm 10511.47$ )	$147.900$ ( $\pm 2.05$ )
DRLS-500	$2.067$ ( $\pm 0.09$ )	$683252.019$ ( $\pm 5948.62$ )	$195.167$ ( $\pm 1.44$ )
DRLS-1000	$2.000$ ( $\pm 0.00$ )	$679665.728$ ( $\pm 4926.88$ )	$214.100$ ( $\pm 0.47$ )
DRLS-2000	$2.000$ ( $\pm 0.00$ )	$680744.118$ ( $\pm 3872.39$ )	$215.967$ ( $\pm 0.07$ )
DRSD	$257.000$ ( $\pm 0.00$ )	$681772.359$ ( $\pm 10348.32$ )	$150.033$ ( $\pm 2.12$ )
BAA
DRLS-100	$247.000$ ( $\pm 5.17$ )	$248677.974$ ( $\pm 1238.19$ )	$100.000$ ( $\pm 0.00$ )
DRLS-250	$240.233$ ( $\pm 5.74$ )	$249173.795$ ( $\pm 770.38$ )	$250.000$ ( $\pm 0.00$ )
DRLS-500	$236.067$ ( $\pm 6.28$ )	$249827.472$ ( $\pm 499.59$ )	$500.000$ ( $\pm 0.00$ )
DRLS-1000	$229.067$ ( $\pm 6.03$ )	$250640.029$ ( $\pm 355.53$ )	$1000.000$ ( $\pm 0.00$ )
DRLS-2000	$219.900$ ( $\pm 5.36$ )	$251405.806$ ( $\pm 243.54$ )	$2000.000$ ( $\pm 0.00$ )
DRSD	$316.367$ ( $\pm 31.53$ )	$250235.142$ ( $\pm 737.21$ )	$315.367$ ( $\pm 31.53$ )
STORM
DRLS-100	$11.667$ ( $\pm 0.51$ )	$15742456.082$ ( $\pm 12191.85$ )	$100.000$ ( $\pm 0.00$ )
DRLS-250	$11.167$ ( $\pm 0.52$ )	$15781724.510$ ( $\pm 8753.59$ )	$250.000$ ( $\pm 0.00$ )
DRLS-500	$11.733$ ( $\pm 0.59$ )	$15797019.902$ ( $\pm 5345.83$ )	$500.000$ ( $\pm 0.00$ )
DRLS-1000	$11.767$ ( $\pm 0.52$ )	$15806575.387$ ( $\pm 3771.91$ )	$1000.000$ ( $\pm 0.00$ )
DRLS-2000	$11.900$ ( $\pm 0.46$ )	$15817039.917$ ( $\pm 2502.88$ )	$2000.000$ ( $\pm 0.00$ )
DRSD	$516.200$ ( $\pm 108.53$ )	$15786864.662$ ( $\pm 9155.49$ )	$515.200$ ( $\pm 108.53$ )

Table 1: Comparison of results obtained from DR L-Shaped method and DRSD method

Table 1 shows the average number of iterations, the average objective function value, and the average number of unique observations in Columns 2, 3, and 4, respectively. The values in the parenthesis are the half-widths of the corresponding confidence intervals. Notice that for DRSD, the number of iterations is also equal to the number of observations used to approximate the ambiguity set. The results show the ability of DRSD to dynamically determine the number of observations by assessing the progress made during the algorithm. The objective function estimate obtained using the DRSD is comparable to the objective function estimate obtained using the DR L-shaped method of comparable size. For instance, the DRSD objective function estimate for STORM that is based upon a sample of size $516.2$ (on average) is within $0.1\%$ of the objective function value estimate of DRLS-500. These results show that the optimal objective function estimate obtained from DRSD are comparable to those obtained using an external sampling-based approach.

PGP
Method	Total	Master	Subproblem	Optimality	Argmax	Separation
DRLS-100	$0.0517$	$0.0019$	$0.0422$	$0.0000$	$-$	$0.0023$
DRLS-250	$0.0770$	$0.0021$	$0.0651$	$0.0000$	$-$	$0.0025$
DRLS-500	$0.0959$	$0.0020$	$0.0822$	$0.0000$	$-$	$0.0026$
DRLS-1000	$0.1214$	$0.0021$	$0.1050$	$0.0000$	$-$	$0.0030$
DRLS-2000	$0.1463$	$0.0022$	$0.1270$	$0.0000$	$-$	$0.0033$
DRSD	$0.3498$	$0.1225$	$0.0549$	$0.0003$	$0.0018$	$0.0697$
CEP
DRLS-100	$0.0154$	$0.0002$	$0.0086$	$0.0000$	$-$	$0.0004$
DRLS-250	$0.0237$	$0.0001$	$0.0124$	$0.0000$	$-$	$0.0004$
DRLS-500	$0.0283$	$0.0001$	$0.0138$	$0.0000$	$-$	$0.0003$
DRLS-1000	$0.0285$	$0.0001$	$0.0134$	$0.0000$	$-$	$0.0003$
DRLS-2000	$0.0286$	$0.0001$	$0.0137$	$0.0000$	$-$	$0.0003$
DRSD	$0.1154$	$0.0299$	$0.0262$	$0.0001$	$0.0009$	$0.0292$
BAA
DRLS-100	$2.6497$	$0.0778$	$2.1698$	$0.0067$	$-$	$0.1837$
DRLS-250	$6.1376$	$0.0756$	$5.2157$	$0.0078$	$-$	$0.3252$
DRLS-500	$11.1848$	$0.0655$	$9.7257$	$0.0056$	$-$	$0.4624$
DRLS-1000	$21.1755$	$0.0616$	$18.5514$	$0.0043$	$-$	$0.8286$
DRLS-2000	$44.2579$	$0.0685$	$39.2027$	$0.0056$	$-$	$1.2710$
DRSD	$2.0943$	$0.1570$	$0.0577$	$0.0002$	$0.0550$	$0.9003$
STORM
DRLS-100	$0.4336$	$0.0022$	$0.3197$	$0.0001$	$-$	$0.0391$
DRLS-250	$1.0080$	$0.0021$	$0.7429$	$0.0001$	$-$	$0.0885$
DRLS-500	$2.1167$	$0.0025$	$1.5629$	$0.0001$	$-$	$0.1957$
DRLS-1000	$4.3179$	$0.0027$	$3.1438$	$0.0001$	$-$	$0.4571$
DRLS-2000	$9.0015$	$0.0027$	$6.2653$	$0.0001$	$-$	$1.3203$
DRSD	$30.4394$	$0.7815$	$0.3084$	$0.0003$	$0.4781$	$23.8544$

Table 2: Computational time comparison between DR L-shaped and DRSD

Table 2 shows the average total computational time (Column 2) for each instance. The table also includes the average time spent to solve the Master problem (first-stage approximate problem) and the subproblems, to verify the optimality conditions, to complete the argmax procedure (only for DRSD), and to solve the distribution separation problem (Columns 3–7, respectively). The results for small scale instances (PGP and CEP) show that both DRSD and the DR L-shaped method take a fraction of a second, but the computational time for DRSD is higher than the DR L-shaped method for all $N$ . This behavior can be attributed to the fact that (i) the subproblems are relatively easy to solve and the computational effort to solve all the subproblems does not increase significantly with $N$ , and (ii) the DRSD is run for a minimum number of iterations (256) and thereby, contributing to the total time taken to solve master problems and distribution separation problems. This observation is in-line with our computational experience with the SD method for 2-SLPs. It is important to note that, while the computational time for the DR L-shaped method on an individual instance may be lower, but the iterative procedure necessary to identify a sufficient sample size may require solving several instances with an increasing sample size. This may result in a significantly higher cumulative computational time. The DRSD method, and the sequential sampling idea in general, mitigates the need for this iterative process.

On the other hand, for large-scale problems (BAA and STORM), we observe a noticeable increase in the computational time for the DR L-shaped method with an increase in $N$ . A significant portion of this time is spent on solving the subproblems. Since the DRSD solves only two subproblems in each iteration, the time taken to solve the subproblems is significantly less in comparison to the DR L-shaped method for all $N$ where all subproblems corresponding to unique observations are solved in each iteration. Notice that for STORM, the average number of iterations taken by DRSD is at least $43$ times the average number of iterations taken by DR L-shaped for each $N$ . This is reflected in the significant increase in the computational time for solving the master problems and the distributional separation problems. In contrast, for BAA, the average number of DRSD iterations is only $28\%$ higher than the DR L-shaped iterations. As a result, the increase in the computational times due to an increase in the master and distributional separation problem solution times is overshadowed by the computational gains attained by solving only two subproblems in each iteration. Moreover, the computational time associated with solving the distribution separation problem can be reduced by using column-generation procedures that take advantage of the problem structure. Such an implementation was not undertaken for our current experiments and is a fruitful future research avenue.

6 Conclusions

We presented a new decomposition approach for solving two-stage distributionally robust linear programs (2-DRLPs) with a general ambiguity set that is defined using continuous and/or discrete probability distributions with a very large sample space. Since this approach extended the stochastic decomposition approach of Higle and Sen [20] for 2-DRLPs with a singleton ambiguity set, we referred to it as Distributionally Robust Stochastic Decomposition (DRSD) method. The DRSD is a sequential sampling-based approach that allowed sampling within the optimization step where only two second-stage subproblems are solved in each iteration. In contrast, an external sampling procedure utilizes the distributionally robust L-shaped method [1] for solving 2-DRLP with a finite number of scenarios, where in each iteration all subproblems are solved. While the design of DRSD accommodates general ambiguity sets, we provided its asymptotic convergence analysis for a family of ambiguity sets that includes the well-known moment-based and Wasserstein metric-based ambiguity sets. Furthermore, we performed computational experiments to evaluate the efficiency and effectiveness of solving distributionally robust variants of four well-known stochastic programming test problems that have supports of size ranging from $216$ to $10^{81}$ . Based on our results, we observed that the objective function estimates obtained using the DRSD and the DR L-shaped method are statistically comparable. These DRSD estimates are obtained while providing computational improvements that are critical for large-scale problem instances.

Appendix A Proofs

In this appendix, we provide the proofs for the propositions related to the asymptotic behavior of the ambiguity sets approximations defined in §2.1 and the recourse function approximation presented in §2.2.

Proof.

(Proposition 2.1) For $P\in\widehat{\mathfrak{P}}_{\textup{mom}}^{k-1}$ , it is easy to verify that $p^{\prime}=\Theta P$ satisfies the support constraint, viz., $\sum_{\omega\in\Omega^{k}}p^{\prime}(\omega)=1$ . Now consider for $i=1,\ldots,q$ , we have

	$\displaystyle\sum_{\omega\in\Omega^{k}}p^{\prime}(\omega)\psi_{i}(\omega)=~{}$	$\displaystyle\sum_{\omega\in\Omega^{k-1},\omega\neq\omega^{k}}p^{\prime}(\omega)\psi_{i}(\omega)+p^{\prime}(\omega^{k})\psi_{i}(\omega^{k})$
	$\displaystyle=~{}$	$\displaystyle\theta^{k}\sum_{\omega\in\Omega^{k-1},\omega\neq\omega^{k}}p(\omega)\psi_{i}(\omega)+\theta^{k}p(\omega^{k})\psi(\omega^{k})+(1-\theta^{k})\psi_{i}(\omega^{k})$
	$\displaystyle=~{}$	$\displaystyle\theta^{k}\sum_{\omega\in\Omega^{k-1}}p(\omega)\psi_{i}(\omega)+(1-\theta^{k})\psi_{i}(\omega^{k})$
	$\displaystyle=~{}$	$\displaystyle\hat{b}_{i}^{k-1}+(1-\theta^{k})\psi_{i}(\omega^{k})$
	$\displaystyle=~{}$	$\displaystyle\sum_{\omega\in\Omega^{k-1}}\hat{p}^{k-1}(\omega)\psi_{i}(\omega)+(1-\theta^{k})\psi_{i}(\omega^{k})=\hat{b}^{k}_{i}.$

This implies that $\Theta^{k}(P)\in\widehat{\mathfrak{P}}_{\textup{mom}}^{k}$ .

Using Proposition 4 in [39], there exists a positive constant $\chi$ such that

\displaystyle 0\leq\mathbb{H}(\widehat{\mathfrak{P}}_{\textup{mom}}^{k},\mathfrak{P}_{\textup{mom}})\leq\chi\|\hat{\mathbf{b}}^{k}-\mathbf{b}\|.

Here, $\mathbf{b}=(b_{i})_{i=1}^{q}$ and $\hat{\mathbf{b}}^{k}=(\hat{b}_{i}^{k})_{i=1}^{q}$ , and $\|\cdot\|$ denotes the Euclidean norm. Since the approximate ambiguity sets are constructed using independent and identically distributed samples of $\tilde{\omega}$ , using law of large numbers, we have $\hat{b}_{i}^{k}\rightarrow b_{i}$ for all $i=1,\ldots,q$ . This completes the proof. ∎

Proof.

(Proposition 2.2) Consider approximate ambiguity sets $\widehat{\mathfrak{P}}^{k-1}_{\text{w}}$ and $\widehat{\mathfrak{P}}^{k}_{\text{w}}$ of the form given in (20g). Let $P^{\prime}=(p^{\prime}(\omega))_{\omega\in\Omega^{k-1}}\in\widehat{\mathfrak{P}}^{k-1}_{\text{w}}$ , and let the reconstructed probability distribution be denoted by $P$ . We can easily check that $P=\Theta(P^{\prime})$ is indeed a probability distribution. With $P$ fixed, it suffices now to show that the polyhedron

\displaystyle\mathcal{E}(P,\widehat{P}^{k})=\left\{\eta\in\mathbb{R}^{\Omega^{k}\times\Omega^{k}}\left|\begin{array}[]{l}\sum_{\omega^{\prime}\in\Omega^{k}}\eta(\omega,\omega^{\prime})=p(\omega)\qquad\forall\omega\in\Omega^{k},\\ \sum_{\omega\in\Omega^{k}}\eta(\omega,\omega^{\prime})=\hat{p}^{k}(\omega^{\prime})\qquad\forall\omega^{\prime}\in\Omega^{k},\\ \sum_{(\omega,\omega^{\prime})\in\Omega^{k}\times\Omega^{k}}\|\omega-\omega^{\prime}\|\eta(\omega,\omega^{\prime})\leq\epsilon\end{array}\right.\right\}.

(49)

is non-empty. Since $P^{\prime}\in\widehat{\mathfrak{P}}^{k-1}_{\text{w}}$ , there exist $\eta^{\prime}(\omega,\omega^{\prime})$ for all $(\omega,\omega^{\prime})\in\Omega^{k-1}\times\Omega^{k-1}$ such that the constraints in the description of the approximate ambiguity set in (20g) are satisfied. We will do this by analyzing two possibilities,

We encounter a previously seen observation, i.e., $\omega^{k}\in\Omega^{k-1}$ and $\Omega^{k}=\Omega^{k-1}$ . Let $\eta(\omega,\omega^{\prime})=\theta^{k}\eta^{\prime}(\omega,\omega^{\prime})$ for $\omega,\omega^{\prime}\in\Omega^{k-1}$ and $\omega\neq\omega^{\prime}\neq\omega^{k}$ ; and $\eta(\omega^{k},\omega^{k})=\theta^{k}\eta^{\prime}(\omega^{k},\omega^{k})+(1-\theta^{k})$ . We will verify the feasibility of this choice by verifying the three sets of constraints in (49).

	$\displaystyle\sum_{\omega^{\prime}\in\Omega^{k}}\eta(\omega,\omega^{\prime})=~{}$	$\displaystyle\sum_{\omega^{\prime}\in\Omega^{k}\setminus\{\omega^{k}\}}\eta(\omega,\omega^{\prime})+\eta(\omega,\omega^{k})$
	$\displaystyle=~{}$	$\displaystyle\sum_{\omega^{\prime}\in\Omega^{k-1}\setminus\{\omega^{k}\}}\theta^{k}\eta^{\prime}(\omega,\omega^{\prime})+\theta^{k}\eta^{\prime}(\omega,\omega^{k})+\mathbf{1}_{\omega=\omega^{k}}(1-\theta^{k})$
	$\displaystyle=~{}$	$\displaystyle\theta^{k}\bigg{(}\sum_{\omega^{\prime}\in\Omega^{k-1}}\eta^{\prime}(\omega,\omega^{\prime})\bigg{)}+\mathbf{1}_{\omega=\omega^{k}}(1-\theta^{k})$
	$\displaystyle=~{}$	$\displaystyle\theta^{k}p^{\prime}(\omega)+\mathbf{1}_{\omega=\omega^{k}}(1-\theta^{k})=p(\omega)\qquad\forall\omega\in\Omega^{k}.$

	$\displaystyle\sum_{\omega\in\Omega^{k}}\eta(\omega,\omega^{\prime})=~{}$	$\displaystyle\sum_{\omega\in\Omega^{k}\setminus\{\omega^{k}\}}\eta(\omega,\omega^{\prime})+\eta(\omega^{k},\omega^{\prime})$
	$\displaystyle=~{}$	$\displaystyle\sum_{\omega\in\Omega^{k-1}\setminus\{\omega^{k}\}}\theta^{k}\eta^{\prime}(\omega,\omega^{\prime})+\theta^{k}\eta(\omega^{k},\omega)+\mathbf{1}_{\omega^{\prime}=\omega^{k}}(1-\theta^{k})$
	$\displaystyle=~{}$	$\displaystyle\theta^{k}\bigg{(}\sum_{\omega\in\Omega^{k-1}}\eta(\omega,\omega^{\prime})\bigg{)}+\mathbf{1}_{\omega^{\prime}=\omega^{k}}(1-\theta^{k})$
	$\displaystyle=~{}$	$\displaystyle\theta^{k}\hat{p}^{k-1}(\omega^{\prime})+\mathbf{1}_{\omega^{\prime}=\omega^{k}}(1-\theta^{k})=\hat{p}^{k}(\omega^{\prime})\qquad\forall\omega^{\prime}\in\Omega^{k}.$

Finally,

	$\displaystyle\sum_{(\omega,\omega^{\prime})\in\Omega^{k}\times\Omega^{k}}\\|\omega-\omega^{\prime}\\|$	$\displaystyle\eta(\omega,\omega^{\prime})$
	$\displaystyle=~{}$	$\displaystyle\sum_{\begin{subarray}{c}(\omega,\omega^{\prime})\in\Omega^{k-1}\times\Omega^{k-1}\\ \omega\neq\omega^{\prime}\neq\omega^{k}\end{subarray}}\theta^{k}\\|\omega-\omega^{\prime}\\|\eta(\omega,\omega^{\prime})+\\|\omega^{k}-\omega^{k}\\|\eta(\omega^{k},\omega^{k})$
	$\displaystyle=~{}$	$\displaystyle\theta^{k}\bigg{(}\sum_{(\omega,\omega^{\prime})\in\Omega^{k-1}\times\Omega^{k-1}}\\|\omega-\omega^{\prime}\\|\eta(\omega,\omega^{\prime})\bigg{)}\leq\theta^{k}\epsilon\leq\epsilon.$

Since all the three constraints are satisfied, the chosen values for $\eta$ is an element of the polyhedron $\mathcal{E}$ , and therefore, $\mathcal{E}\neq\emptyset$ .

We encounter a new observation, i.e., $\omega^{k}\notin\Omega^{k-1}$ . Let $\eta(\omega,\omega^{\prime})=\theta^{k}\eta^{\prime}(\omega,\omega^{\prime})$ for $\omega,\omega^{\prime}\in\Omega^{k-1}$ , $\eta(\omega^{k},\omega^{\prime})=0$ for $\omega^{\prime}\in\Omega^{k-1}$ , $\eta(\omega,\omega^{k})=0$ for $\omega\in\Omega^{k-1}$ , and $\eta(\omega^{k},\omega^{k})=(1-\theta^{k})$ . Let us verify the three conditions defining (49) with this choice for $\eta$ variables.

	$\displaystyle\sum_{\omega^{\prime}\in\Omega^{k}}\eta(\omega,\omega^{\prime})=~{}$	$\displaystyle\sum_{\omega^{\prime}\in\Omega^{k}\setminus\{\omega^{k}\}}\eta(\omega,\omega^{\prime})+\eta(\omega,\omega^{k})$
	$\displaystyle=~{}$	$\displaystyle\sum_{\omega^{\prime}\in\Omega^{k-1}}\theta^{k}\eta^{\prime}(\omega,\omega^{\prime})+\mathbf{1}_{\omega=\omega^{k}}(1-\theta^{k})$
	$\displaystyle=~{}$	$\displaystyle\theta^{k}p^{\prime}(\omega)+\mathbf{1}_{\omega=\omega^{k}}(1-\theta^{k})=p(\omega).$
	$\displaystyle\sum_{\omega\in\Omega^{k}}\eta(\omega,\omega^{\prime})=~{}$	$\displaystyle\sum_{\omega\in\Omega^{k}\setminus\{\omega^{k}\}}\eta(\omega,\omega^{\prime})+\eta(\omega^{k},\omega^{\prime})$
	$\displaystyle=~{}$	$\displaystyle\sum_{\omega^{\prime}\in\Omega^{k-1}}\theta^{k}\eta^{\prime}(\omega,\omega^{\prime})+\mathbf{1}_{\omega^{\prime}=\omega^{k}}(1-\theta^{k})$
	$\displaystyle=~{}$	$\displaystyle\theta^{k}\hat{p}^{k-1}++\mathbf{1}_{\omega^{\prime}=\omega^{k}}(1-\theta^{k})=\hat{p}^{k}(\omega^{\prime})$

Consider,

	$\displaystyle\sum_{(\omega,\omega^{\prime})\in\Omega^{k}\times\Omega^{k}}\\|\omega-\omega^{\prime}\\|$	$\displaystyle\eta(\omega,\omega^{\prime})$
	$\displaystyle=~{}$	$\displaystyle\sum_{(\omega,\omega^{\prime})\in\Omega^{k-1}\times\Omega^{k-1}}\theta^{k}\\|\omega-\omega^{\prime}\\|\eta^{\prime}(\omega,\omega^{\prime})+\\|\omega^{k}-\omega^{k}\\|\eta(\omega^{k},\omega^{k})$
		$\displaystyle\hskip 28.45274pt+\sum_{\omega\in\Omega^{k}}\\|\omega-\omega^{k}\\|\eta(\omega,\omega^{k})+\sum_{\omega^{\prime}\in\Omega^{k}}\\|\omega^{k}-\omega^{\prime}\\|\eta(\omega^{k},\omega^{\prime})$
	$\displaystyle\leq~{}$	$\displaystyle\theta^{k}\epsilon\leq\epsilon.$

Therefore, the value of the $\eta$ variables satisfy the constraints and $\mathcal{E}\neq\emptyset$ . This implies that $\Theta^{k}(P)\in\widehat{\mathfrak{P}}_{\textup{w}}^{k}$ .

Next, let us consider a distribution $Q\in\widehat{\mathfrak{P}}_{\textup{w}}^{k}$ . Then,

\displaystyle d_{\textup{w}}(Q,P^{*})\leq~{}

\displaystyle d_{\textup{w}}(Q,\widehat{P}^{k})+d_{\textup{w}}(\widehat{P}^{k},P^{*})\leq\epsilon+d_{\textup{w}}(\widehat{P}^{k},P^{*}).

The above inequality is a consequence of the triangle inequality of Wasserstein distance. Since $P\in\widehat{\mathfrak{P}}^{k}_{\textup{w}}$ , we have $d_{\textup{w}}(P^{*},P)\leq\epsilon$ . Under compactness assumption for $\Omega$ , $d>2$ , and $\mathbb{E}_{P^{*}}[\text{exp}(\|\tilde{\omega}\|^{a})]<\infty$ , Theorem 2 in [16] guarantees

\displaystyle\text{Prob}\big{[}d_{\textup{w}}(\widehat{P}^{k},P^{*})\leq\delta\big{]}\leq\left\{\begin{array}[]{ll}C~{}\text{exp}(-ck\delta^{d})&\text{if}~{}\delta>1\\ C~{}\text{exp}(-ck\delta^{a})&\text{if}~{}\delta\leq 1\end{array}\right.

for all $k\geq 1$ . This implies that the $\lim_{k\rightarrow\infty}d_{\textup{w}}(\widehat{P}^{k},P^{*})=0$ , almost surely. Consequently, we obtain that $d_{\textup{w}}(Q,P^{*})\leq\epsilon$ (or equivalently $Q\in\mathfrak{P}_{\textup{w}}$ ) as $k\rightarrow\infty$ , almost surely. This completes the proof. ∎

Proof.

(Proposition 2.3) Recall that $\mathcal{X}\times\Omega$ is a compact set because of Assumptions (A1)) and (A4), and $\{Q^{k}\}$ is a sequence of continuous (piecewise linear and convex) functions. Further, the construction of the set of dual vertices satisfies $\Pi^{0}=\{{\bf 0}\}\subseteq\ldots\subseteq\Pi^{k}\subseteq\Pi^{k+1}\subseteq\ldots\subseteq\Pi$ which ensures that $0\leq Q^{k}(x,\omega)\leq Q^{k+1}(x,\omega)\leq Q(x,\omega)$ for all $(x,\omega)\in(\mathcal{X},\Omega)$ . Since $\{Q^{k}\}$ increases monotonically and is bounded by a finite function $Q$ (due to (A2)), this sequence pointwise converges to some function $\xi(x,\omega)\leq Q(x,\omega)$ . Once again due to (A2), we know that the set of dual vertices $\Pi$ is finite and since $\Pi^{k}\subseteq\Pi^{k+1}\subseteq\Pi$ , the set $\lim_{k\rightarrow\infty}\Pi^{k}:=\overline{\Pi}~{}(\subseteq\Pi)$ is also a finite set. Clearly,

\displaystyle\xi(x,\omega)=\lim_{k\rightarrow\infty}Q^{k}(x,\omega)=\max~{}\{\pi^{\top}[r(\omega)-T(\omega)x]~{}|~{}\pi\in\overline{\Pi}\}

is the optimal value of a LP, and hence, is a continuous function. The compactness of $\mathcal{X}\times\Omega$ , and continuity, monotonicity and pointwise convergence of $\{Q^{k}\}$ to $\xi$ guarantees that the sequence uniformly converges to $\xi$ (Theorem 7.13 in [32]). ∎

References

[1] M. Bansal, K. Huang, and S. Mehrotra. Decomposition Algorithms for Two-Stage Distributionally Robust Mixed Binary Programs. SIAM Journal on Optimization, pages 2360–2383, January 2018.
[2] M. Bansal and S. Mehrotra. On solving two-stage distributionally robust disjunctive programs with a general ambiguity set. European Journal of Operational Research, 279(2):296–307, December 2019.
[3] G. Bayraksan and D. K. Love. Data-Driven Stochastic Programming Using Phi-Divergences. In The Operations Research Revolution, INFORMS TutORials in Operations Research, pages 1–19. INFORMS, September 2015.
[4] G. Bayraksan and D.P. Morton. Assessing solution quality in stochastic programs. Mathematical Programming, 108(2):495–514, Sep 2006.
[5] G. Bayraksan and D.P. Morton. A sequential sampling procedure for stochastic programming. Operations Research, 59(4):898–913, 2011.
[6] A. Ben-Tal, D. den Hertog, A. De Waegenaere, B. Melenberg, and G. Rennen. Robust Solutions of Optimization Problems Affected by Uncertain Probabilities. Management Science, 59(2):341–357, November 2012.
[7] D. Bertsimas, X. V. Doan, K. Natarajan, and C. Teo. Models for minimax stochastic linear optimization problems with risk aversion. Mathematics of Operations Research, 35(3):580–602, 2010.
[8] D. Bertsimas and I. Popescu. Optimal inequalities in probability theory: A convex optimization approach. SIAM Journal on Optimization, 15(3):780–804, 2005.
[9] J. R. Birge and F. Louveaux. Introduction to Stochastic Programming. Springer Series in Operations Research and Financial Engineering. Springer, 2011.
[10] M. Breton and S. El Hachem. Algorithms for the solution of stochastic dynamic minimax problems. Computational Optimization and Applications, 4(4):317–345, October 1995.
[11] M. Breton and S. El Hachem. A scenario aggregation algorithm for the solution of stochastic dynamic minimax problems. Stochastics and Stochastic Reports, 53(3-4):305–322, May 1995.
[12] G. B. Dantzig and P. Wolfe. Decomposition principle for linear programs. Operations Research, 8(1):101–111, 1960.
[13] E. Delage and Y. Ye. Distributionally robust optimization under moment uncertainty with application to data-driven problems. Operations Research, 58:595–612, 2010.
[14] J. Dupacová. The minimax approach to stochastic programming and an illustrative application. Stochastics, 20:73–88, 1987.
[15] E. Erdogan and G. Iyengar. Ambiguous chance constrained problems and robust optimization. Mathematical Programming, 107(1-2):37–61, 2006.
[16] N. Fournier and A. Guillin. On the rate of convergence in wasserstein distance of the empirical measure. Probability Theory and Related Fields, 162(3-4):707–738, 2015.
[17] H. Gangammanavar, Y. Liu, and S. Sen. Stochastic decomposition for two-stage stochastic linear programs with random cost coefficients. INFORMS Journal on Computing, 2020.
[18] H. Gangammanavar and S. Sen. Stochastic dynamic linear programming: A sequential sampling algorithm for multistage stochastic linear programming. Technical report, Southern Methodist University, 2019. Available on Optimization Online: http://www.optimization-online.org/DB_HTML/2019/09/7395.html.
[19] G. A. Hanasusanto and D. Kuhn. Conic programming reformulations of two-stage distributionally robust linear programs over wasserstein balls. Operations Research, 66(3):849–869, 2018.
[20] J. L. Higle and S Sen. Stochastic decomposition: An algorithm for two-stage linear programs with recourse. Mathematics of Operations Research, 16(3):650–669, 1991.
[21] J. L. Higle and S Sen. Finite master programs in regularized stochastic decomposition. Mathematical Programming, 67(1-3):143–168, 1994.
[22] J. L. Higle and S Sen. Stochastic Decomposition: A Statistical Method for Large Scale Stochastic Linear Programming. Kluwer Academic Publishers, Boston, MA., 1996.
[23] J. L. Higle and S Sen. Statistical approximations for stochastic linear programming problems. Annals of Operations Research, 85(0):173–193, 1999.
[24] R. Jiang and Y. Guan. Risk-averse two-stage stochastic program with distributional ambiguity. Operations Research, 66(5):1390–1405, 2018.
[25] Sanjay Mehrotra and He Zhang. Models and algorithms for distributionally robust least squares problems. Mathematical Programming, 148(1–2):123–141, 2014.
[26] P. Mohajerin Esfahani and D. Kuhn. Data-driven distributionally robust optimization using the wasserstein metric: performance guarantees and tractable reformulations. Mathematical Programming, 171(1):115–166, Sep 2018.
[27] H. Rahimian and S. Mehrotra. Distributionally robust optimization: A review. arXiv preprint arXiv:1908.05659, 2019.
[28] M. Riis and K. A. Andersen. Applying the minimax criterion in stochastic recourse programs. European Journal of Operational Research, 165(3):569–584, September 2005.
[29] R. T. Rockafellar and Roger J.-B. Wets. Scenarios and Policy Aggregation in Optimization Under Uncertainty. Mathematics of Operations Research, 16(1):119–147, February 1991.
[30] W. Römisch. Stability of stochastic programming problems. Handbooks in operations research and management science, 10:483–554, 2003.
[31] J. O. Royset and R. Szechtman. Optimal budget allocation for sample average approximation. Operations Research, 61(3):762–776, 2013.
[32] W. Rudin. Principles of mathematical analysis. McGraw-Hill Book Co., New York, third edition, 1976. International Series in Pure and Applied Mathematics.
[33] H. Scarf. A min-max solution of an inventory problem. In Studies in the Mathematical Theory of Inventory and Production, chapter 12, pages 201–209. RAND Corporation, Santa Monica CA, 1958.
[34] S. Sen and Z. Zhou. Multistage Stochastic Decomposition: A bridge between Stochastic Programming and Approximate Dynamic Programming. SIAM Journal on Optimization, 24(1):127–153, 2014.
[35] A. Shapiro and S. Ahmed. On a class of minimax stochastic programs. SIAM Journal on Optimization, 14(4):1237–1249, 2004.
[36] A. Shapiro and S. Ahmed. On a class of minimax stochastic programs. SIAM Journal on Optimization, 14(4):1237–1249, 2004.
[37] A. Shapiro, D. Dentcheva, and A. Ruszczynski. Lectures on Stochastic Programming: Modeling and Theory, Second Edition. Society for Industrial and Applied Mathematics, Philadelphia, PA, USA, 2014.
[38] A. Shapiro and A. Kleywegt. Minimax analysis of stochastic problems. Optimization Methods and Software, 17(3):523–542, 2002.
[39] H. Sun and H. Xu. Convergence analysis for distributionally robust optimization and equilibrium problems. Mathematics of Operations Research, 41(2):377–401, May 2016.
[40] R. M. Van Slyke and R. J. B. Wets. L-shaped linear programs with applications to optimal control and stochastic programming. SIAM Journal on Applied Mathematics, 17(4):638–663, 1969.
[41] W. Xie. Tractable reformulations of distributionally robust two-stage stochastic programs with $\infty-$ wasserstein distance. arXiv preprint arXiv:1908.08454, 2019.
[42] C. Zhao and Y. Guan. Data-driven risk-averse two-stage stochastic program with $\zeta$ -structure probability metrics. Available at http://www.optimization-online.org/DB_FILE/2015/07/5014, 2015.

	$\displaystyle\|f(x)-g^{k}(x)\|=$	$\displaystyle~{}\bigg{\|}\max_{P\in\mathfrak{P}}\mathbb{E}_{P}[Q(x,\tilde{\omega})]-\max_{P^{\prime}\in\widehat{\mathfrak{P}}^{k}}\mathbb{E}_{P^{\prime}}[Q(x,\tilde{\omega})]\bigg{\|}\qquad\forall x\in\mathcal{X}$
	$\displaystyle\leq$	$\displaystyle~{}\max_{P\in\mathfrak{P}}\min_{P^{\prime}\in\widehat{\mathfrak{P}}^{k}}\big{\|}\mathbb{E}_{P}[Q(x,\tilde{\omega})]-\mathbb{E}_{P^{\prime}}[Q(x,\tilde{\omega}]\big{\|}\qquad\forall x\in\mathcal{X}$
	$\displaystyle\leq$	$\displaystyle~{}\max_{P\in\mathfrak{P}}\min_{P^{\prime}\in\widehat{\mathfrak{P}}^{k}}\sup_{x\in\mathcal{X}}\big{\|}\mathbb{E}_{P}[Q(x,\tilde{\omega})]-\mathbb{E}_{P^{\prime}}[Q(x,\tilde{\omega})]\big{\|}.$

	$\displaystyle\|f(\bar{x})-g^{k_{n}}(\hat{x}^{k_{n}})\|\leq$	$\displaystyle~{}\|f(\bar{x})-g^{k_{n}}(\bar{x})\|+\|g^{k_{n}}(\bar{x})-g^{k_{n}}(\hat{x}^{k_{n}})\|$
	$\displaystyle\leq$	$\displaystyle~{}\mathbb{H}(\mathfrak{P},\widehat{\mathfrak{P}}^{k})+(\\|c\\|+C)\\|\bar{x}-\hat{x}^{k_{n}}\\|.$