On Distributionally Robust Multistage Convex Optimization: Data-driven Models and Performance

In this section, we present a data-driven model for DR-MCO and some of its properties. Let $\mathcal{T}\mathrel{\mathop{\ordinarycolon}}=\{1,\dots,T\}$ denote the set of stage indices. In each stage $t\in\mathcal{T}$, we use $\mathcal{X}_{t}\subset\mathbb{R}^{d_{t}}$ to denote the convex state space and $x_{t}$ its elements, which is known as the state vector. We denote the set of uncertainties before stage $t$ as $\Xi_{t}\subseteq\mathbb{R}^{\delta_{t}}$ and its elements as $\xi_{t}$. For simplicity, we use the notation $\mathcal{X}_{0}=\{x_{0}\}$ and $\Xi_{1}=\{\xi_{1}\}$ to denote parameter sets of the given initial state. After each stage $t\in\mathcal{T}\setminus\{T\}$, the uncertainty $\xi_{t}$ is assumed to be distributed according to an unknown probability measure $p_{t}$ taken from a subset of the probability measures $\mathcal{P}_{t}\subset\mathcal{M}^{\mathrm{Prob}}(\Xi_{t})$ supported on $\Xi_{t}$. The cost in each stage $t\in\mathcal{T}$ is described through a nonnegative, lower semicontinuous local cost function $f_{t}(x_{t-1},x_{t};\xi_{t})$, which is assumed to be convex in $x_{t-1}$ and $x_{t}$ for every $\xi_{t}\in\Xi_{t}$. We allow $f_{t}$ to take $+\infty$ to model constraints relating the states $x_{t-1},x_{t}$ and the uncertainty $\xi_{t}$ (so that $\mathcal{X}_{t}$ can be independent of $x_{t-1}$ and $\xi_{t}$). The DR-MCO can be written in a nested formulation as follows.

Here, $\mathbb{E}_{\xi_{t}\sim p_{t}}$ is the expectation with respect to variable $\xi_{t}$ distributed according to the probability measure $p_{t}$. We remark that both the uncertainty sets $\Xi_{t}$ and the ambiguity sets $\mathcal{P}_{t}$ are independent between stages, which are usually referred to as stagewise independence.

Based on stagewise independence in the nested formulation, we can write the following recursion that is equivalent to (1) using the (worst-case expected) cost-to-go functions,

for each $t\in\mathcal{T}$ and we set by convention $\mathcal{Q}_{T}(x_{T})\mathrel{\mathop{\ordinarycolon}}=0$ for any $x_{T}\in\mathcal{X}_{T}$. To simplify the notation, we also define the following value function for each stage $t\in\mathcal{T}$:

Using these value functions, we may write the optimal value of the DR-MCO (1) as $Q_{1}(x_{0};\xi_{1})$ and further simplify the recursion (2) as

While there are many different choices of the ambiguity set $\mathcal{P}_{t}$ for each stage $t\in\mathcal{T}$ (see e.g., [47]), we focus on the data-driven ambiguity sets constructed as follows. Suppose we have the knowledge of $n_{t}$ samples $\hat{\xi}_{t,1},\dots,\hat{\xi}_{t,n_{t}}$ of the uncertainty $\xi_{t}$. The empirical probability measure is given by $\hat{\nu}_{t}\mathrel{\mathop{\ordinarycolon}}=\frac{1}{n_{t}}\sum_{k=1}^{n_{t}}\delta_{\hat{\xi}_{t,k}}$, where for each $k=1,\dots,n_{t}$, $\delta_{\hat{\xi}_{t,k}}$ is the Dirac probability measure supported at the point $\hat{\xi}_{t,k}\in\Xi_{t}$, i.e., $\int_{\Xi_{t}}f\operatorname{d\!}\delta_{\hat{\xi}_{t,k}}=f(\hat{\xi}_{t,k})$ for any compactly supported function $f$ on $\Xi_{t}$. Such an empirical probability measure $\hat{\nu}_{t}$ captures the information from the sample data and is often used to build the sample average approximation for multistage stochastic optimization [43].

Fix any distance function $d_{t}(\cdot,\cdot)$ on $\Xi_{t}$, the Wasserstein distance of order 1 (a.k.a, Kantorovich-Rubinstein distance) is defined as

for any two probability measures $\mu,\nu\in\mathcal{M}^{\mathrm{Prob}}(\Xi_{t})$, where $P^{i}_{*}(\pi)$ is the pushforward measure induced by the projection maps $P^{i}\mathrel{\mathop{\ordinarycolon}}\Xi_{t}\times\Xi_{t}\to\Xi_{t}$ by sending $P^{i}(\xi^{1},\xi^{2})=\xi^{i}$, for $i=1$ or $2$. That is, the joint probability measure $\pi$ in (5) has marginal probability measures equal to $\mu$ and $\nu$.

It can be shown that $W_{t}$ is indeed a distance on the space of probability measures $\mathcal{M}^{\mathrm{Prob}}(\Xi_{t})$ [46, Definition 6.1] except that it may take the value of $+\infty$. Thus it is natural to restrict our attention to the convex subset of probability measures with finite distance to a Dirac measure on $\Xi_{t}$

Note that any continuous function $g(\xi)$ that satisfies $\lvert g(\xi)\rvert\leq C(1+d_{t}(\bar{\xi},\xi))$ for some $C>0$ and $\bar{\xi}\in\Xi_{t}$ would be integrable for any probability measure in $\mathcal{W}_{t}$. Now given any such continuous functions $g_{t,1},\dots,g_{t,m_{t}}$ on $\Xi_{t}$ and a real vector $\rho_{t}\mathrel{\mathop{\ordinarycolon}}=(\rho_{t,j})_{t=0}^{m_{t}}\in\mathbb{R}^{m_{t}+1}$, we define the Wasserstein ambiguity set $\mathcal{P}_{t}$ as

The first inequality constraint in the definition (7) bounds the Wasserstein distance of the probability measure $p\in\mathcal{W}_{t}$ from the empirical measure $\hat{\nu}_{t}$, while the second set of constraints on $p$ are defined as $\langle g_{t,j},p\rangle\mathrel{\mathop{\ordinarycolon}}=\int_{\Xi_{t}}g_{t,j}(\xi)\operatorname{d\!}p(\xi)$, which could be bounds on the moments.

It is well studied that the recursion for data-driven Wasserstein DR-MCO problems (2) can be reformulated as a finite-dimensional minimum-supremum optimization problem, under the assumption of strict feasibility.

We provide the derivation and proof details in Section A. As a corollary, we can prove a special version of the Kantorovich-Rubinstein duality formula [46, Remark 6.5].

A major motivation for using Wasserstein DR-MCO models is the out-of-sample performance guarantee, which ensures that the decisions evaluated on the true probability distribution would perform no worse than the in-sample training with high probability. To begin with, we say that a probability measure $\mu\in\mathcal{W}_{t}$ is sub-Gaussian if $\int_{\Xi_{t}}\exp(C\xi^{2})\operatorname{d\!}\mu(\xi)<+\infty$ for some constant $C>0$, or it has finite third moments if $\int_{\Xi_{t}}\lVert\xi\rVert^{3}\operatorname{d\!}\mu(\xi)<+\infty$. Our discussion is based on the following specialized version of the measure concentration inequality [16, Theorem 2].

The measure concentration bound in (9) becomes slightly more intricate when the dimension of the uncertainty $\delta_{t}=2$ (see the details in [16]), so we focus our discussion below on the other cases.

The out-of-sample performance refers to the evaluation of the solutions and policies obtained from solving DR-MCO (1) on the true probability measures $\nu_{t}$ for each $t\in\mathcal{T}$. To be precise, fix any optimal policy given by the DR-MCO (1), i.e., an optimal initial stage state $x_{1}^{*}\in\operatorname*{arg\,min}_{x_{1}\in\mathcal{X}_{1}}\{f_{1}(x_{0},x_{1};\xi_{1})+\mathcal{Q}_{1}(x_{1})\}$, and a collection of mappings $x_{t}^{*}\mathrel{\mathop{\ordinarycolon}}\mathcal{X}_{t-1}\times\Xi_{t}\to\mathcal{X}_{t}$ for $t=2,\dots,T$, such that

1. 1.

   $x_{t}^{*}(x_{t-1},\cdot)$ is Borel measurable for any $x_{t-1}\in\mathcal{X}_{t-1}$;

2. 2.

   $x_{t}^{*}(x_{t-1},\xi_{t})\in\operatorname*{arg\,min}_{x_{t}\in\mathcal{X}_{t}}\{f_{t}(x_{t-1},x_{t};\xi_{t})+\mathcal{Q}_{t}(x_{t})\}$ for any $x_{t-1}\in\mathcal{X}_{t-1}$ and $\xi_{t}\in\Xi_{t}$.

We define recursively the costs in the out-of-sample evaluation cost-to-go functions as

with $\mathcal{Q}_{T}^{\mathrm{Eval}}(x_{T})=0$ for any $x_{T}\in\mathcal{X}_{T}$. Then the out-of-sample evaluation mean cost associated with this policy is defined as $v^{\mathrm{Eval}}\mathrel{\mathop{\ordinarycolon}}=f_{1}(x_{0},x_{1}^{*};\xi_{1})+\mathcal{Q}_{1}^{\mathrm{Eval}}(x_{1}^{*})$. The next theorem provides a lower bound on the probability that the event $v^{\mathrm{Eval}}\leq Q_{1}(x_{0};\xi_{1})$ happens, which is often known as the out-of-sample performance guarantee.

While Theorem 4 is a direct consequence of Theorem 3, it shows some interesting aspects of DR-MCO using Wasserstein ambiguity sets. First, to get a certain probabilistic bound for the out-of-sample performance guarantee, we may need to increase the number of samples $n_{t}$ or the Wasserstein distance bound $\rho_{t,0}$ for a larger number of stages $T$. Second, for probability measures that are not sub-Gaussian (or more generally those with heavy tails, see [15]), for which we would like to apply the second condition, such increase require the product $n_{t}\cdot\rho_{t,0}^{2}$ to grow approximately on the order of $\mathcal{O}(T)$ when $\alpha$ is close to 1. Therefore, it is sometimes useful to take larger values of $\rho_{t,0}$ for out-of-sample performance guarantees, especially when the number of sample $n_{t}$ is limited, or when the true probability measures $\nu_{t}$ are not sub-Gaussian.

Another well studied approach to guarantee out-of-sample performance is the multistage robust convex optimization (MRCO) model [18, 48]. However, MRCO considers only the worst-case outcomes of the uncertainties and thus can be overly conservative. To be precise, we define the nominal MSCO from the empirical measures $\hat{\nu}_{t}$ without risk aversion by the following recursion

where

and $\mathcal{Q}_{T}^{\mathrm{Nomin}}(x_{T})=0$ for any $x_{T}\in\mathcal{X}_{T}$. Any MRCO that is built directly from data could have a much larger optimal cost than the nominal MSCO with a probability growing with the numbers of samples $n_{t}$ and stages $T$, as illustrated by the following example.

In contrast, the difference of the optimal values between the DR-MCO and the nominal MSCO $Q_{1}(x_{0};\xi_{1})-Q_{1}^{\mathrm{Nomin}}(x_{0};\xi_{1})$ can be bounded by the Wasserstein distances $\rho_{t,0}$, as shown below.

Theorem 5 shows that the conservatism of the DR-MCO can be adjusted linearly with the Wasserstein distance bound $\rho_{t,0}$, assuming the Lipschitz continuity of the value functions in the uncertainties. However, as the out-of-sample performance guarantee in Theorem 4 depends on some unknown constants in Theorem 3, it is not easy to numerically determine the Wasserstein distance bounds $\rho_{t,0}$. We discuss some practical choices of the bounds $\rho_{t,0}$ in Section 4.

Recall that for any convex function $\mathcal{Q}\mathrel{\mathop{\ordinarycolon}}\mathcal{X}\to\mathbb{R}\cup\{+\infty\}$, an affine function $\mathcal{V}\mathrel{\mathop{\ordinarycolon}}\mathcal{X}\to\mathbb{R}$ is called a (valid) linear cut if $\mathcal{Q}(x)\geq\mathcal{V}(x)$ for all $x\in\mathcal{X}$. A collection of such valid linear cuts $\{\mathcal{V}^{j}\}_{1\leq j\leq i}$ defines a valid under-approximation $\underline{\mathcal{Q}}^{i}(x)\mathrel{\mathop{\ordinarycolon}}=\max_{1\leq j\leq i}\mathcal{V}^{j}(x)$ of $\mathcal{Q}(x)$. Similarly by convexity, given a collection of overestimate values $v^{j}\geq\mathcal{Q}(x^{j})$ for $j=1,\dots,i$, we can define a valid over-approximation by the convex envelope $\overline{\mathcal{Q}}^{i}(x)\mathrel{\mathop{\ordinarycolon}}=\mathrm{conv}_{1\leq j\leq i}(v^{j}+\iota_{x^{j}}(x))$, where $\iota_{x^{j}}(x)=0$ when $x=x^{j}$ and $+\infty$ otherwise, is the convex indicator function centered at $x^{j}$. The validness of these approximations $\underline{\mathcal{Q}}(x)\leq\mathcal{Q}(x)\leq\overline{\mathcal{Q}}(x)$ for all $x\in\mathcal{X}$ suggests that we may use them in the place of $\mathcal{Q}$ for recursive updates during a stagewise decomposition algorithm.

To see how the recursive updates may work, let us assume temporarily in this subsection that the worse-case probability measure $p_{t}^{*}$ in (4) exists and can be found. Given any under-approximation $\underline{\mathcal{Q}}_{t}$ for $\mathcal{Q}_{t}$, for any $\hat{\xi}_{t}\in\Xi_{t}$, we can generate a linear cut $V_{t}(x_{t-1};\xi_{t})\mathrel{\mathop{\ordinarycolon}}=\underline{\hbox{$v$}}(\hat{\xi}_{t})+u(\hat{\xi}_{t})^{\mathsf{T}}(x_{t-1}-\tilde{x}_{t-1})$ for the value function $Q_{t}(x_{t-1};\hat{\xi}_{t})$ with some $\tilde{x}_{t-1}\in\mathcal{X}_{t-1}$, where $\underline{\hbox{$v$}}(\hat{\xi}_{t})$ and $u(\hat{\xi}_{t})$ are the optimal value and an optimal dual solution to the following Lagrangian dual of the under-approximation problem that is parametrized by $\hat{\xi}_{t}$:

assuming that $u(\hat{\xi}_{t})$ exists. Then we can aggregate these linear cuts $V_{t}(x_{t-1};\xi_{t})$ into a linear cut $\mathcal{V}_{t-1}(\cdot)\mathrel{\mathop{\ordinarycolon}}=\mathbb{E}_{\xi_{t}\sim p_{t}^{*}}V_{t}(\cdot;\xi_{t})$ for $\mathcal{Q}_{t-1}$, where the expectation is taken componentwise with respect to the probability measure $p_{t}^{*}$. In this way, we can use the under-approximation $\underline{\mathcal{Q}}_{t}$ for $\mathcal{Q}_{t}$ to update the under-approximation $\underline{\mathcal{Q}}_{t-1}$ for $\mathcal{Q}_{t-1}$ by the aggregated linear cut $\mathcal{V}_{t-1}$.

Likewise, given any over-approximation $\overline{\mathcal{Q}}_{t}$ for $\mathcal{Q}_{t}$, for any $\hat{\xi}_{t}\in\Xi_{t}$ and $\tilde{x}_{t-1}\in\mathcal{X}_{t-1}$, we can solve the following over-estimation problem

which gives an overestimate value $\bar{v}(\hat{\xi}_{t})\geq Q_{t}(\tilde{x}_{t-1};\hat{\xi}_{t})$. Again we can aggregate the overestimate value by setting $v_{t-1}\mathrel{\mathop{\ordinarycolon}}=\mathbb{E}_{\xi_{t}\sim p_{t}^{*}}\bar{v}(\xi_{t})$, which by definition satisfies $v_{t-1}\geq\mathcal{Q}_{t-1}(\tilde{x}_{t-1})$ and thus can be used to update the over-approximation $\overline{\mathcal{Q}}_{t-1}$ for $\mathcal{Q}_{t-1}$.

There are, however, some potential issues with this recursive approximation method. First, the supremum in the Lagrangian dual problem (14) may not be attained, which could happen if $\underline{\hbox{$v$}}(\hat{\xi}_{t})=+\infty$ (i.e., $\tilde{x}_{t-1}$ is an infeasible state for $\hat{\xi}_{t}$) or if any neighborhood of $\tilde{x}_{t-1}$ contains such an infeasible state. In this case, we may fail to generate a linear cut $V_{t}(\cdot;\hat{\xi}_{t})$. Second, Lipschitz constants of the linear cuts $\mathcal{V}_{t-1}$ may be affected by the under-approximation $\underline{\mathcal{Q}}_{t}$, causing a worse approximation quality. This may happens when $\tilde{x}_{t-1}$ is an extreme point of $\mathcal{X}_{t-1}$. In fact, it is shown in [48] that the Lipschitz constants of $\underline{\mathcal{Q}}_{t}$ could exceed those of $\mathcal{Q}_{t}$ and grow with the total number of stages $T$. Third, the over-approximation function evaluates to $+\infty$ at any point that is not in the convex hull of previously visited points, which makes the gap $\overline{\mathcal{Q}}_{t}(x)-\underline{\mathcal{Q}}_{t}(x)$ less useful as an estimate of the quality of the current solutions or policies.

To remedy these issues, we consider a technique called Lipschitzian regularization. Given regularization factors $M_{t}>0$, we define the regularized local cost function as

and the regularized value function

recursively for $t=T,T-1,\dots,2$, where $\mathcal{Q}_{t}^{\mathrm{R}}$ is the regularized expected cost-to-go function defined as

for $t\leq T-1$, and $\mathcal{Q}_{T}^{\mathrm{R}}(x_{T})\equiv 0$ for any $x_{T}\in\mathcal{X}_{T}$. It is then straightforward to check that $Q_{t}^{\mathrm{R}}(x_{t-1};\xi_{t})$ is uniformly $M_{t}$-Lipschitz continuous in $x_{t-1}$ for each $\xi_{t}\in\Xi_{t}$, and consequently $\mathcal{Q}_{t-1}^{\mathrm{R}}$ is also $M_{t}$-Lipschitz continuous. In the definitions (14) and (15), we can replace accordingly the original cost functions $f_{t}$ with the regularized cost function $f_{t}^{\mathrm{R}}$, which would guarantee the $M_{t}$-Lipschitz continuity of the generated cuts and allow us to enhance the over-approximation to be $M_{t}$-Lipschitz continuous.

Lipschitzian regularization in general only gives under-approximations of the true value and expected cost-to-go functions. We need the following assumption to preserve the optimality and feasibility of the solutions.

We remark by the following proposition that Assumption 2 can be satisfied in any problem that already have uniformly Lipschitz continuous value function $Q_{t}(\cdot;\xi_{t})$ for all $\xi_{t}\in\Xi_{t}$.