On Distributionally Robust Chance Constrained Programs with Wasserstein Distance

Introduction

We study distributional robust chance constrained programs (DRCCPs) of the form:

We denote the feasible region induced by DRCC (1c) as

2 Assumptions

In this paper, we consider Wasserstein ambiguity set ${\mathcal{P}}$, i.e., we make the following assumption on the ambiguity set ${\mathcal{P}}$.

The Wasserstein ambiguity set ${\mathcal{P}}$ is defined as

This assumption has been studied in recent DRCCP literature ;

By making this assumption, it might cause the DRCC (1c) to be more conservative than the general setting studied in ;

In practice, one needs to choose a proper Wasserstein radius $\delta$ through cross validation to alleviate the over-conservatism caused by Assumption (A2), which will be illustrated in Section 5.

Finally, we suppose that Assumptions (A1) and (A2) hold throughout the paper.

3 Related Literature

When DRCC set $Z$ is not convex, many inner convex approximations have been proposed. In , the authors proposed to aggregate the multiple uncertain constraints with positive scalars in to a single constraint, and then use conditional value-at-risk (${\bf{CVaR}}$) approximation scheme to develop an inner approximation of $Z$. This approximation is shown to be exact for single DRCCP when ${\mathcal{P}}$ is specified by first and second moments in or, more generally, by convex moment constraints in . In , the authors provided several sufficient conditions under which the well-known Bonferroni approximation of joint DRCCP is exact and yields a convex reformulation.

4 Contributions

In this paper, we study approximations and exact reformulations of DRCCP under Wasserstein ambiguity set. In particular, our main contributions are summarized as below.

We derive a deterministic equivalent reformulation for set $Z$ and show that this reformulation admits a conditional value-at-risk (${\bf{CVaR}}$) interpretation, i.e.,

where $f(\cdot,\cdot)$ is defined in Theorem 1.

We show that set $Z$, once bounded, is mixed integer representable with big-M coefficients and $N$ additional binary variables.

We derive inner and outer approximations based upon ${\bf{CVaR}}$ interpretation. We develop compact formulations for these approximations and compare their strengths.

When the decision variables are pure binary (i.e., $S\subseteq\{0,1\}^{n}$), we first show that the nonlinear constraints in the reformulation can be recast as submodular knapsack constraints. Then, by exploiting the polyhedral properties of submodular functions, we propose a new big-M free mixed integer linear reformulation, which can be effectively solved by a branch and cut algorithm.

The remainder of the paper is organized as follows. Section 2 presents exact reformulations of DRCC set $Z$. Section 3 provides inner and outer approximations of set $Z$ and compares their strengths. Section 4 studies binary DRCCP (i.e., $S\subseteq\{0,1\}^{n}$), develops a big-M free formulation. Section 5 numerically illustrates the proposed methods. Section 6 concludes the paper.

Exact Reformulations

In this section, we will show that DRCC set $Z$ admits a conditional-value-at-risk (${\bf{CVaR}}$) interpretation and is mixed integer representable. This reformulation also allows us to derive tight inner and outer approximations in next section.

In this subsection, we will reformulate the set $Z$ into its deterministic counterpart with respect to empirical distribution. The main idea of this reformulation is that we first use the strong duality result from to formulate the worst-case chance constraint into its dual form, and then break down the indicator function according to its definition.

and ${\mathcal{I}}(\bm{x})=\emptyset$ if $\bm{a}(\bm{x})\neq 0$ and ${\mathcal{I}}(\bm{x})=[I]$, otherwise, and characteristic function $\chi_{\mathcal{R}}(\bm{x})=\infty$ if $\bm{x}\notin\mathcal{R}$ and 0, otherwise.

We separate the proof into three steps, where the first step is to apply strong duality result for distributionally robust optimization, the second step is to break down the indicator function, and the third step is to replace the dual variable with its reciprocal.

Next, we break down the indicator function in the infimum of (6b) by discussing the conditions under which it is equal to zero or one and reformulate it as below.

For given $\lambda\geq 0$ and $\bm{\zeta}\in{\mathcal{Z}}$, we have

Therefore, we only need to show that for any $i\in[I]$,

If $\bm{a}(\bm{x})^{\top}\bm{\zeta}_{i}>b_{i}(\bm{x})$, then in the left-hand side of (6d), the infimum is equal to $-1$ by letting $\bm{\xi}:=\bm{\zeta}$, which equals the right-hand side since the infimum is also achieved by $\bm{\xi}:=\bm{\zeta}$.

If $\bm{a}(\bm{x})^{\top}\bm{\zeta}_{i}\leq b_{i}(\bm{x})$, then for any $\bm{\xi}\in\Xi$, we either have $\bm{a}(\bm{x})^{\top}\bm{\xi}_{i}>b_{i}(\bm{x})$ or $\bm{a}(\bm{x})^{\top}\bm{\xi}_{i}\leq b_{i}(\bm{x})$. Hence, the left-hand side of (6d) is equivalent to

where $\inf_{\bm{a}(\bm{x})^{\top}\bm{\xi}_{i}\leq b_{i}(\bm{x})}\left[{\|\bm{\xi}-\bm{\zeta}\|}\right]=0$ by letting $\bm{\xi}:=\bm{\zeta}$.

Finally, let $Z^{\prime}$ denote the set in the right-hand side of (4) , we only need to show that $Z=Z^{\prime}$.

Given $\bm{x}\in Z$, there exists $\lambda\geq 0$ such that $(\bm{x},\lambda)$ satisfies (4). If $\lambda>0$, then let $\gamma=\frac{1}{\lambda}$. Then it is easy to see that $(\bm{x},\gamma)$ satisfies (4) . Hence, $\bm{x}\in Z^{\prime}$.

Now suppose that $\lambda=0$, then in (4), we have

Similarly, given $\bm{x}\in Z^{\prime}$, there exists $\gamma\geq 0$ such that $(\bm{x},\gamma)$ satisfies (4) . If $\gamma>0$, then let $\lambda=\frac{1}{\gamma}$. Then it is easy to see that $(\bm{x},\lambda)$ satisfies (4). Hence, $\bm{x}\in Z$.

Now suppose that $\gamma=0$, then in (4) , we have

for each $j\in[N]$. Thus, (4) reduces to $\delta\leq 0$ contradicting that $\delta>0$.∎

Please note that in the proof, we use the fact that $\delta>0$ from Assumption (A1), and the formulation (4) does not hold if $\delta=0$.

while its $(1-\epsilon)$-conditional value-at-risk (CVaR) is defined as

With the definitions above, we observe that set $Z$ in (4) has a ${\bf{CVaR}}$ interpretation.

First, we observe that the constraint in (4) directly implies $\gamma\geq 0$, thus the nonnegativity constraint of $\gamma$ can be dropped, i.e., equivalently, we have

Next, in the above formulation, letting $\gamma^{\prime}:=-\gamma$ and replacing the existence of $\gamma^{\prime}$ by finding the best $\gamma^{\prime}$ such that the constraint still holds, we arrive at

In the following sections, we will derive the inner and outer approximations mainly based upon ${\bf{CVaR}}$ formulation in Corollary 1.

2 Exact Mixed Integer Program Reformulation

In this subsection, we show that set $Z$ is mixed integer representable. To do so, we first observe that the reformulation of set $Z$ in Theorem 1 can be further simplified as a disjunction of a nonconvex set and a convex set.

We need to show that $Z_{1}\cup Z_{2}\subseteq Z$ and $Z\subseteq Z_{1}\cup Z_{2}$.

Given $\bm{x}\in Z_{2}$, we have ${\mathcal{I}}(\bm{x})=[I]$, thus $f(\bm{x},\bm{\zeta})$ (defined in (5)) is $\infty$. Thus, let $\gamma=\frac{\delta}{\epsilon}$. Clearly, $(\gamma,\bm{x})$ satisfies all the constraints in (4), i.e., $\bm{x}\in Z$. Hence, $Z_{2}\subseteq Z$.

Given $\bm{x}\in Z_{1}$, there exists $(\gamma,\nu,\bm{z},\bm{x})$ which satisfies constraints in (8). Suppose that ${\mathcal{I}}(\bm{x})=[I]$, then we have $\bm{a}(\bm{x})=\bm{0}$. Hence, for each $i\in{\mathcal{I}}(\bm{x})$, we have (8a) and (8b) imply that

That is, $b_{i}(\bm{x})>0$. Thus, $\bm{x}\in Z_{2}\subseteq Z$.

Now we suppose that ${\mathcal{I}}(\bm{x})=\emptyset$. For each $i\in[I]$, (8a) and (8b) along with $\nu>0$ imply that

where the second inequality is due to (8d). Then according to (8a), we have

i.e., $(\gamma/\nu,\bm{x})$ satisfies the constraints in (4), i.e., $\bm{x}\in Z$. Thus, $Z_{1}\subseteq Z$.

Similarly, given $\bm{x}\in Z$, there exists $(\gamma,\bm{x})$ which satisfies constraints in (4). Suppose that $\bm{a}(\bm{x})=\bm{0}$, then we must have $b_{i}(\bm{x})\geq 0$ for all $i\in[I]$, otherwise, we have $f(\bm{x},\bm{\zeta}^{j})=0$ for all $j\in[I]$. Then (4a) is equivalent to

a contradiction that $\gamma\geq 0,\epsilon\in(0,1)$. Hence, we must $\bm{x}\in Z_{2}$.

From now on, we assume that $\bm{a}(\bm{x})\neq 0$. Let us define $\widehat{\gamma}=\gamma\|\bm{a}(\bm{x})\|_{*},\nu=\|\bm{a}(\bm{x})\|_{*}$, and $z_{j}=\min_{i\in[I]}(\max\{b_{i}(\bm{x})-\bm{a}(\bm{x})^{\top}\bm{\zeta}_{i}^{j},0\}-\widehat{\gamma},0)$ for each $j\in[N]$. Clearly, $(\widehat{\gamma},\nu,\bm{z},\bm{x})$ satisfies constraints in (8), i.e., $\bm{x}\in Z_{1}$. ∎

We make the following remarks about the disjunctive formulation of set $Z$.

For DRCCP with left-hand uncertainty (i.e., $\eta_{1}=1,\eta_{2}=0$), we have

For DRCCP with right-hand uncertainty or two-side uncertainty (i.e., $\eta_{1}\in\{0,1\},\eta_{2}=1$), we have $Z_{2}=\emptyset$.

According to Lemma 2 , the feasible region induced by a chance constraint is closed, so is set $Z$. However, set $Z_{1}$ might not be closed due to $\nu>0$ in (8e). In practice, one can find a lower bound $0<\underline{\nu}$ such that

or let $\underline{\nu}$ be a sufficiently small number. Then replace the constraint $\nu>0$ in (8e) by $\nu\geq\underline{\nu}$.

We observe that set $Z_{1}$ can be formulated as a mixed integer set when it is bounded, i.e., we can use binary variables to represent the nonlinear constraints (8b) as mixed integer linear ones. This result has been observed independently by (see their Proposition 1) for single DRCCP.

for all $j\in[N]$. Then $Z_{1}$ is mixed integer representable, i.e.,

We first observe that the constraints (8b) are equivalent to

Above, the outer maximum in the right-hand side can be linearized by using a binary variable $y_{j}$, a continuous variable $s_{j}$, and big-M coefficient $M_{j}$ for each $j\in[N]$. By doing so, we arrive at (10). ∎∎

Usually, we can derive the big-M coefficients by inspection; for example, suppose that $\bm{x}\in[\bm{L},\bm{U}]$, then for each $j\in[N]$, we can find $M_{j}$ in the following way: (i) rewrite $b_{i}(\bm{x})-\bm{a}(\bm{x})^{\top}\bm{\zeta}_{i}^{j}=\sum_{\tau\in[n]}(B_{i\tau}-\eta_{1}\zeta_{i\tau}^{j})x_{\tau}+b^{i}-\eta_{2}\zeta_{i(n+1)}^{j}$ for each $i\in[I]$, (ii) define sets $\widehat{S}_{+}=\{\tau\in[n]:B_{i\tau}-\eta_{1}\zeta_{i\tau}^{j}>0\}$ and $\widehat{S}_{-}=[n]\setminus\widehat{S}_{+}$, and (iii) let $M_{j}$ be

There are various methods introduced in literature to further tighten big-M coefficients.

Formulation (10) involves $N$ binary variables and big-M coefficients. In Section 4, we will show that for binary DRCCP, set $Z_{1}$ can be reformulated as a big-M free formulation without introducing additional binary variables.

3 A Special Case: DRCCP with Right-hand Uncertainty

In this subsection, we consider DRCCP with right-hand uncertainty, i.e., $\eta_{1}=0,\eta_{2}=1,\bm{a}(\bm{x})=\bm{e}_{n+1}$. We first observe that when $\bm{a}(\bm{x})=\bm{e}_{n+1}\neq\bm{0}$, in Theorem 1, set $Z$ of DRCCP with right-hand uncertainty has a more compact representation.

For DRCCP with right-hand uncertainty (i.e., $\eta_{1}=0,\eta_{2}=1,\bm{a}(\bm{x})=\bm{e}_{n+1}$), set $Z$ is equivalent to the following mathematical program:

The result directly follows from Theorem 1.∎∎

The differences between this result and the one in Proposition 1 are: (i) for DRCCP with right-hand uncertainty, we do not need to reformulate set $Z$ as a disjunction of two sets, and (ii) compared to set $Z_{1}$ in (8), there is no need to introduce additional positive variable $\nu$ in the formulation (11).

Following the similar derivation in Theorem 2, we can also reformulate the set $Z$ in (11) as a mixed integer program as below. This result has been observed independently by (see their Proposition 2).

for all $j\in[N],i\in[I]$. Then set $Z$ is mixed integer representable, i.e.,

The proof is similar as that of Theorem 2, thus is omitted. ∎∎

Similar to Theorem 2, suppose that $\bm{x}\in[\bm{L},\bm{U}]$, then for each $j\in[N]$, one possible $M_{j}$ can be derived as below:

where $\widehat{S}_{+}=\{\tau\in[n]:B_{i\tau}>0\}$ and $\widehat{S}_{-}=[n]\setminus\widehat{S}_{+}$.

Outer and Inner Approximations

In this section, we will introduce one outer approximation and three different inner approximations by exploiting the exact reformulations in the previous section. The outer approximation can provide a lower bound for DRCCP, while inner approximations can provide good-quality feasible solutions. Our numerical study in Section 5 will demonstrate that together these approximations, we can obtain better solutions than those from the exact mixed integer programming model in the previous section, in particular, for large-sized instances.

Therefore, in Corollary 1, if we replace ${\bf{CVaR}}_{1-\epsilon}\left(\cdot\right)$ by ${\bf{VaR}}_{1-\epsilon}\left(\cdot\right)$, then we have the following outer approximation of set $Z$.

According to Corollary 1 and the well-known result in that

and ${\mathcal{I}}(\bm{x})=\emptyset$ if $\bm{a}(\bm{x})\neq 0$, otherwise, ${\mathcal{I}}(\bm{x})=[I]$. Thus we further have

Using the fact that $\frac{\delta}{\epsilon}>0$, we arrive at (13).∎∎

We make the following remarks about outer approximation $Z_{{\bf{VaR}}}$.

A particular interpretation of formulation (13) is that in order to enforce the robustness, we further penalize the left-hand side of uncertain constraints by the dual norm $\|\bm{a}(\bm{x})\|_{*}$; and

Recently, there are several works on distributionally robust optimization with $\infty-$Wasserstein ambiguity set, and set $Z_{{\bf{VaR}}}$ is in fact equal to the feasible region induced by DRCC with $\infty-$Wasserstein ambiguity set.

Consider $\infty-$Wasserstein ambiguity set ${\mathcal{P}}^{W}$ defined as

where $\infty-$Wasserstein distance is defined as

Then set $Z_{{\bf{VaR}}}$ is equivalent to

This result demonstrates that set $Z_{{\bf{VaR}}}$ indeed can be viewed as a deterministic counterpart of DRCCP with $\infty-$Wasserstein ambiguity set. Thus, in practice, it can serve as an alternative for the set $Z$.

For the completeness of this paper, we present the mixed integer program formulation of outer approximation set $Z_{{\bf{VaR}}}$. The proof is omitted as it directly follows the proof of Theorem 2.

for all $j\in[N]$. Then $Z_{{\bf{VaR}}}$ is mixed integer representable, i.e.,

Similar to Theorem 2, suppose that $\bm{x}\in[\bm{L},\bm{U}]$, then for each $j\in[N]$, one possible $M_{j}$ can be derived as below:

where $\widehat{S}_{+}=\{\tau\in[n]:B_{i\tau}-\eta_{1}\zeta_{i\tau}^{j}>0\}$ and $\widehat{S}_{-}=[n]\setminus\widehat{S}_{+}$.

2 Inner Approximation I- Robust Scenario Approximation

Using the definition of $f(\bm{x},\bm{\zeta})$ and the fact that $\frac{\delta}{\epsilon}>0$, we arrive at (16). ∎∎

We remark that set $Z_{R}$ in (16) is very similar to scenario approach to regular chance constrained program . That is, we generate $N$ i.i.d. samples $\{\bm{\zeta}^{j}\}_{j\in[N]}$ and enforce all the sampled constraints to hold. It has been shown in that if $N$ is larger than a threshold, it guarantees with high probability that the solution of scenario approach is feasible to the regular chance constrained program. Different from scenario approach, in formulation (16), we add a penalty $\frac{\delta}{\epsilon}\|\bm{a}(\bm{x})\|_{*}$ to the sampled constraints, which can be viewed as a “robust” scenario approach to the regular chance constrained problem. That is, if the sample size $N$ is not sufficiently large (i.e., $N$ is smaller than the threshold given by ), one might want to add a penalty $\frac{\delta}{\epsilon}\|\bm{a}(\bm{x})\|_{*}$ to enforce that set $Z_{R}$ is indeed a subset of the feasible region induced by a regular chance constraint.

3 Inner Approximation II- An Inner Chance Constrained Programming Approximation

Next we propose an inner chance constrained programming approximation of set $Z$ by constructing a feasible $\gamma$ in (4).

According to the definition of $f(\bm{x},\bm{\zeta})$ in (5) and the fact $\frac{\delta}{\epsilon}>0$, set $Z_{I}$ is equivalent to

For any $\bm{x}\in Z_{I}$, we need to show that $\bm{x}\in Z$. Since $\bm{x}\in Z_{I}$, there exists an $\alpha$ such that $(\bm{x},\alpha)$ satisfies constraints in (18). Now let us define $\gamma=\frac{\delta}{\epsilon-\alpha}$. It remains to show that $(\gamma,\bm{x})$ satisfies the constraints (4).

where the first inequality is due to $\gamma=\frac{\delta}{\epsilon-\alpha}$ and the second inequality is due to (18). Hence,

where the first inequality is due to $f(\bm{x},\bm{\zeta}^{j})\geq 0$ according to its definition in (5) and the second inequality is due to $|{\mathcal{C}}|\leq N\alpha$.∎∎

We also observe that (i) set $Z_{R}$ is a special case of set $Z_{I}$ by letting $\alpha=0$, thus, we must have $Z_{R}\subseteq Z_{I}$; (ii) there are $\lceil N\epsilon\rceil$ disjoint intervals that $\alpha$ belong to, that is,

The feasible region induced by the above chance constraint increases if we decrease the value of $\alpha$ to $\frac{i-1}{N}$. Therefore, to optimize over set $S\cap Z_{I}$, we only need to enumerate these $\lceil N\epsilon\rceil$ different values of $\alpha$ and choose the one which yields the smallest objective value; (iii) for each given $\alpha$, the chance constraint in (17) is mixed integer program representable. These three results are summarized below.

set $Z_{I}=\cup_{\alpha\in\left\{0,\frac{1}{N},\ldots,\frac{\lceil N\epsilon\rceil-1}{N}\right\}}Z_{I}^{\alpha}$, where set $Z_{I}^{\alpha}$ is defined as

for each $\alpha\in\left\{0,\frac{1}{N},\ldots,\frac{\lceil N\epsilon\rceil-1}{N}\right\}$; and

for all $j\in[N]$, then set $Z_{I}^{\alpha}$ is mixed integer representable, i.e.,

Similar to Theorem 2, suppose that $\bm{x}\in[\bm{L},\bm{U}]$, then for each $j\in[N]$ and $\alpha\in\left\{0,\frac{1}{N},\ldots,\frac{\lceil N\epsilon\rceil-1}{N}\right\}$, one possible $M_{j}^{\alpha}$ in Corollary 5 can be derived as below:

where $\widehat{S}_{+}=\{\tau\in[n]:B_{i\tau}-\eta_{1}\zeta_{i\tau}^{j}>0\}$ and $\widehat{S}_{-}=[n]\setminus\widehat{S}_{+}$.

According to Corollary 5, to solve the inner approximation of DRCCP (i.e., $\min_{\bm{x}\in S\cap Z_{I}}$), we can solve $\min_{\bm{x}\in S\cap Z_{I}^{\alpha}}$ for each $\alpha\in\left\{0,\frac{1}{N},\ldots,\frac{\lceil N\epsilon\rceil-1}{N}\right\}$ and choose the smallest value.

where the inequality is due to $N\epsilon+1\geq\lceil N\epsilon\rceil$ and $\lceil c_{1}N^{1-\frac{c_{2}}{2}}\rceil\geq c_{1}N^{1-\frac{c_{2}}{2}}$. This observation is summarized below.

According to , any light-tail distribution (e.g., Gaussian distribution) satisfies the assumption in above proposition; and

Sets $Z_{{\bf{VaR}}}$ and $Z_{I}$ together build up a hierarchy of regular chance constrained programs, which converges to DRCC set $Z$ as $N\rightarrow\infty$ and preserves the outer and inner approximations, i.e., $Z_{I}\subseteq Z\subseteq Z_{{\bf{VaR}}}$ for all $N$ and $Z_{I}\rightarrow Z$ and $Z_{{\bf{VaR}}}\rightarrow Z$ as $N\rightarrow\infty$.

4 Inner Approximation III- 𝐂𝐕𝐚𝐑𝐂𝐕𝐚𝐑{\bf{CVaR}} Approximation

In this subsection, we will study a well-known convex approximation of a chance constraint, which is to replace the nonconvex chance constraint by a convex constraint defined by ${\bf{CVaR}}$ (cf., ). For DRCC set $Z$, the resulting approximation is

Set $Z_{{\bf{CVaR}}}$ (21) is convex and is an inner approximation of set $Z$. The following results show a reformulation of set $Z_{{\bf{CVaR}}}$. We would like to acknowledge that this result has been independently observed by a recent work in . Thus, the proof is omitted.

Set $Z_{{\bf{CVaR}}}\subseteq Z$ is equivalent to

We remark that one can directly derive the equivalent form (22) of set $Z_{{\bf{CVaR}}}$ based upon formulation (4).

Since $\max\{b_{i}(\bm{x})-\bm{a}(\bm{x})^{\top}\bm{\zeta}_{i},0\}\geq b_{i}(\bm{x})-\bm{a}(\bm{x})^{\top}\bm{\zeta}_{i}$, by replacing $\max\{b_{i}(\bm{x})-\bm{a}(\bm{x})^{\top}\bm{\zeta}_{i},0\}$ with $b_{i}(\bm{x})-\bm{a}(\bm{x})^{\top}\bm{\zeta}_{i}$, then function $f(\bm{x},\bm{\zeta})$ is lower bounded by

Thus, set $Z$ can be inner approximated by the following set

By introducing additional variables $\bm{z}$ to linearize the nonlinear function $\min\left\{\underline{f}(\bm{x},\bm{\zeta})-\gamma,0\right\}$, we arrive at

This set can be proven to be exactly equal to set $Z_{{\bf{CVaR}}}$ by discussing whether $\|\bm{a}(\bm{x})\|_{*}>0$ or not:

if $\|\bm{a}(\bm{x})\|_{*}>0$, then replace $\delta,\{z_{j}\}_{j\in[N]}$ and $\gamma$ by $\delta\|\bm{a}(\bm{x})\|_{*},\{z_{j}\|\bm{a}(\bm{x})\|_{*}\}_{j\in[N]}$ and $\gamma\|\bm{a}(\bm{x})\|_{*}$;

if $\|\bm{a}(\bm{x})\|_{*}=0$, since $\delta\leq\epsilon\gamma+\frac{1}{N}\sum_{j\in[N]}z_{j}$ and $\delta>0,\epsilon>0$, according to the pigeonhole principle, we must have $z_{j_{0}}+\gamma>0$ for some $j_{0}\in[N]$, which implies that $b_{i}(\bm{x})\geq 0$ for all $i\in[I]$.

This observation inspires us that $Z_{{\bf{CVaR}}}=Z$ if $N\epsilon\leq 1$. In fact, if $N\epsilon\leq 1$, then we must have $f(\bm{x},\bm{\zeta}^{j})=\underline{f}(\bm{x},\bm{\zeta}^{j})$ for all $j\in[N]$, which implies that $Z_{{\bf{CVaR}}}=Z$.

Suppose that $\epsilon\in(0,1/N]$, then $Z=Z_{{\bf{CVaR}}}$.

We note that $Z_{{\bf{CVaR}}}\subseteq Z=Z_{1}\cup Z_{2}$, where $Z_{1}$ and $Z_{2}$ are defined in (8) and (9), respectively. We note that set $Z_{2}\subseteq Z_{{\bf{CVaR}}}$. Indeed, suppose that $\bm{x}\in Z_{2}$, i.e., $\bm{a}(\bm{x})=\bm{0},b_{i}(\bm{x})\geq 0$ for each $i\in[I]$, then let $\nu=0,\gamma=0$ and $z_{j}=0$ for each $j\in[N]$. Clearly, $(\nu,\gamma,\bm{z},\bm{x})$ satisfies the constraints in (22). Hence, $\bm{x}\in Z_{{\bf{CVaR}}}$.

Thus, it is sufficient to show that $Z_{1}\subseteq Z_{{\bf{CVaR}}}$. Indeed, given $\bm{x}\in Z_{1}$, there exists $(\nu,\gamma,\bm{z})$ such that $(\nu,\gamma,\bm{z},\bm{x})$ satisfies the constraints in (8). We only need to show that $z_{j}+\gamma>0$ for each $j\in[N]$. Suppose that there exists a $j_{0}\in[N]$ such that $z_{j_{0}}+\gamma\leq 0$. Then according to (8a), we have

where the second inequality is due to $\epsilon N\leq 1$ and $z_{j_{0}}+\gamma\leq 0$, a contradiction that $\delta>0$. Therefore, in (8b), we must have

for each $i\in[I],j\in[N]$. Hence, $(\nu,\gamma,\bm{z},\bm{x})$ satisfies the constraints in (22), i.e., $\bm{x}\in Z_{{\bf{CVaR}}}$. ∎∎

The result in Proposition 5 shows that if the risk parameter $\epsilon$ is small enough (i.e., less than or equal to $\frac{1}{N}$), then set $Z$ is convex and is equivalent to its ${\bf{CVaR}}$ approximation.

5 Formulation Comparisons

First, we would like to compare sets $Z_{R},Z_{{\bf{CVaR}}}$. Indeed, we can show that $Z_{R}\subseteq Z_{{\bf{CVaR}}}$, i.e., set $Z_{R}$ is at least as conservative as ${\bf{CVaR}}$ approximation $Z_{{\bf{CVaR}}}$.

Let $Z_{R},Z_{{\bf{CVaR}}}$ be defined in (16), (22) , respectively. Then $Z_{R}\subseteq Z_{{\bf{CVaR}}}.$

Given $\bm{x}\in Z_{R}$, we only need to show that $\bm{x}\in Z_{{\bf{CVaR}}}$. Indeed, let us consider $\nu=\|\bm{a}(\bm{x})\|_{*}$, $\gamma=\frac{\delta}{\epsilon}\|\bm{a}(\bm{x})\|_{*},z_{j}=0$ for all $j\in[N]$, then we see that $(\nu,\gamma,\bm{z},\bm{x})$ satisfies the constraints in (22), i.e., $\bm{x}\in Z_{{\bf{CVaR}}}$.∎

The following example illustrates sets $Z,Z_{{\bf{VaR}}},Z_{{\bf{CVaR}}},Z_{R},Z_{I}$ and their inclusive relationships.

Suppose $N=3,n=2,I=2,\delta=1/6,\epsilon=2/3$ and $\bm{\zeta}_{1}^{1}=(0,0,\sqrt{2})^{\top},\bm{\zeta}_{2}^{1}=(0,0,3\sqrt{2})^{\top},\bm{\zeta}_{1}^{2}=(0,0,3\sqrt{2})^{\top},\bm{\zeta}_{2}^{2}=(0,0,\sqrt{2})^{\top},\bm{\zeta}_{1}^{3}=(0,0,3\sqrt{2})^{\top},\bm{\zeta}_{2}^{3}=(0,0,2\sqrt{2})^{\top},\bm{a}(\bm{x})=\bm{e}_{3}=\begin{pmatrix}0,0,1\end{pmatrix}^{\top},b_{1}(x)=x_{1},b_{2}(x)=x_{2}$. Then, (2) becomes:

Clearly, we have Z_{R}\subsetneq\left\{\begin{subarray}{c}Z_{{\bf{CVaR}}}\\ \rotatebox[origin={c}]{-90.0}{\not\subseteq}\\ Z_{I}\end{subarray}\right\}\subsetneq Z\subsetneq Z_{{\bf{VaR}}} (see Figure 1 for an illustration).

Finally, the theoretical inclusive relationships of sets $Z,Z_{{\bf{VaR}}},Z_{R},Z_{I},Z_{{\bf{CVaR}}}$ are shown in Figure 2 and their reformulations are summarized in Table 1.

DRCCP with Pure Binary Decision Variables

In this section, we will study DRCCP with pure binary decision variables $\bm{x}\in\{0,1\}^{n}$, i.e., we assume that $S\subseteq\{0,1\}^{n}$. If $S$ is a bounded integer set, we can use binary expansion to reformulate $S$ an equivalent binary set (c.f., ). For binary DRCCP, we will show that the reformulations in the previous section can be improved.

Our main derivation of stronger formulations is based upon some polyhedral results of submodular functions, which will be briefly reviewed in this subsection.

We first briefly introduce the definition of submodularity and interested readers are referred to for more details.

for every $T_{1},T_{2}\subseteq[n]$ with $T_{1}\subseteq T_{2}$ and every $t\in[n]\setminus T_{2}$, we must have $g(T_{1}\cup\{t\})-g(T_{1})\geq g(T_{2}\cup\{t\})-g(T_{2})$.

We first begin with the following lemmas on submodular functions.

Since $\bm{d}_{1}^{\top}\bm{x}+d_{2}$ is a nondecreasing submodular function and $-\max\left(t,d_{3}\right)$ is a nonincreasing concave function, the submodularity of their composition follows by Table 1 in .∎∎

Given $q\geq 1$, function $f(\bm{x})=\|\bm{x}\|_{q}$ with $q\geq 1$ is submodular over the binary hypercube.

This is because $f(\bm{x})=\|\bm{x}\|_{q}=\sqrt[q]{\sum_{l\in[n]}x_{l}}$, and $g(\bm{e}^{\top}\bm{x})$ is a submodular function if $g(\cdot)$ is a concave function (cf., ). ∎∎

Next, we will introduce polyhedral properties of submodular functions. For any given submodular function $f(\bm{x})$ with $\bm{x}\in\{0,1\}^{n}$, let us denote $\Pi_{f}$ to be its epigraph, i.e.,

Then the convex hull of $\Pi_{f}$ is characterized by the system of “extended polymatroid inequalities” (EPI) , i.e.,

where $\Omega$ denotes a collection of all permutations of set $[n]$ and $\rho_{\sigma_{l}}=f(\bm{e}_{A_{l}^{\sigma}})-f(\bm{e}_{A_{l-1}^{\sigma}})$ for each $l\in[n]$ with $A_{0}^{\sigma}=\emptyset,A_{l}^{\sigma}=\{\sigma_{1},\ldots,\sigma_{l}\}$ and $(\bm{e}_{T})_{\tau}=\begin{cases}1,&\text{if }\tau\in T\\ 0,&\text{if }\tau\in[n]\setminus T\end{cases}$.

In addition, although there are $n!$ number of inequalities in (24), these inequalities can be easily separated by a greedy procedure.

2 Reformulating a Binary DRCCP by Submodular Knapsack Constraints: Big-M free

In this section, we will replace the nonlinear constraints defining the feasible region of a binary DRCCP (i.e., set $S\cap Z$) by submodular knapsack constraints. These constraints can be equivalently described by the system of EPI in (24). Therefore we obtain a big-M free mixed integer representation of set $S\cap Z$.

First, we introduce $n$ auxiliary variables complementing binary variables $\bm{x}$, denoted by $\bm{w}$, i.e., $w_{l}+x_{l}=1$ for each $l\in[n]$. With these $n$ additional variables, we can reformulate function $b_{i}(\bm{x})-\bm{a}(\bm{x})^{\top}\bm{\zeta}_{i}^{j}$ as

Thus, from above discussion, we can formulate $S\cap Z$ (recall that set $Z=Z_{1}\cup Z_{2}$ according to Proposition 1) as the following mixed integer set with submodular knapsack constraints.

Suppose that $S\subseteq\{0,1\}^{n}$. Then $S\cap Z=(S\cap\widehat{Z}_{1})\cup(S\cap Z_{2})$, where

According to Proposition 1, equalities (25) and the fact that $\bm{a}(\bm{x})=\begin{pmatrix}\eta_{1}\bm{x}\\ \eta_{2}\end{pmatrix}$ with constant $\eta_{1},\eta_{2}\in\{0,1\}$, constraints (8b) and (8d) are equivalent to (26b) and (26d). Thus, we only need to show that $S\cap Z_{1}\subseteq(S\cap\widehat{Z}_{1})\cup(S\cap Z_{2})$. There are two cases.

If $\eta_{2}=1$, then we must have $\|\bm{a}(\bm{x})\|_{*}=\left\|\begin{pmatrix}\eta_{1}\bm{x}\\ \eta_{2}\end{pmatrix}\right\|_{*}\geq 1$, then $S\cap Z_{1}=S\cap\widehat{Z}_{1}$. We are done.

If $\eta_{2}=0$, then we must have $\eta_{1}=1$. For any $\bm{x}\in S\cap Z_{1}$, we need to show that $\bm{x}\in(S\cap\widehat{Z}_{1})\cup(S\cap Z_{2})$. If $\bm{x}=0$, then the constraints (8) become

Since $\nu>0,\delta>0,1>\epsilon>0$, thus by the pigeonhole principle, we must have $z_{j_{0}}+\gamma>0$ for some $j_{0}\in[N]$. This implies that $b_{i}(\bm{x})>0$ for each $i\in[I]$. Together with $\bm{a}(\bm{x})=\bm{0}$, we must have $\bm{x}=0\in S\cap Z_{2}$.

Now suppose that $\bm{x}\neq 0$. Note that $S\cap Z_{1}\subseteq\{0,1\}^{n}$, therefore, $\bm{x}\neq 0$ implies that $\|\bm{x}\|_{*}\geq 1$, thus, $v\geq\left\|\begin{pmatrix}\eta_{1}\bm{x}\\ \eta_{2}\end{pmatrix}\right\|_{*}=\|\bm{x}\|_{*}\geq 1$. Thus, $\bm{x}\in S\cap\widehat{Z}_{1}$.∎

From the proof of Theorem 7, we note that if $b^{i}\geq\frac{\delta}{\epsilon}$ for each $i\in[I]$, then we have $S\cap Z_{2}\subseteq S\cap\widehat{Z}_{1}$. Thus, $S\cap Z=S\cap\widehat{Z}_{1}$.

Suppose that $S\subseteq\{0,1\}^{n}$ and $b^{i}\geq\frac{\delta}{\epsilon}$ for each $i\in[I]$. Then $S\cap Z=S\cap\widehat{Z}_{1}$.

From the proof of Theorem 7, we only need to show that $\bm{x}=\bm{0}\in S\cap\widehat{Z}_{1}$. In this case, we have $\bm{w}=\bm{e}-\bm{x}=\bm{e}$. Then according to (25), we have $\bm{r}_{ij}^{\top}\bm{x}+\bm{t}_{ij}^{\top}\bm{w}+u_{ij}=b_{i}(\bm{x})-\bm{a}(\bm{x})^{\top}\bm{\zeta}_{i}^{j}=b^{i}$. Let us set $\nu=1,\gamma=\frac{\delta}{\epsilon},\bm{z}=\bm{0}$. Then it is easy to see that $(\bm{x},\bm{w},\bm{z},\gamma,\nu)$ satisfies the constraints in (26), i.e., $\bm{0}\in S\cap\widehat{Z}_{1}$. ∎∎

We note that the left-hand sides of constraints (26b) and (26d) are submodular functions according to Lemma 1 and Lemma 2. Therefore, equivalently, we can replace these constraints with the convex hulls of epigraphs of their associated submodular functions. Thus, we arrive at the following equivalent representation of set $S\cap\widehat{Z}_{1}$.

Suppose that $S\subseteq\{0,1\}^{n}$ and $\|\cdot\|$ is $L_{p}$ norm with $p\geq 1$. Then

Note that the optimization problem $\min_{\bm{x}\in S\cap Z_{1}}\bm{c}^{\top}\bm{x}$ can be solved by a branch and cut algorithm. In particular, at each branch and bound node, denoted as $(\widehat{\bm{x}},\widehat{\bm{w}},\widehat{\bm{z}},\widehat{\gamma},\widehat{\nu})$, there might be too many (i.e., $N\times I+1$) valid inequalities to add, since in (28b) and (28d), there are $N\times I+1$ convex hulls of epigraphs (i.e., $\left\{\operatorname{conv}(\Pi_{ij})\right\}_{i\in[I],j\in[N]},\operatorname{conv}(\Pi_{0})$) to be separated from. Therefore, instead, we can first check and find the epigraphs of $\kappa$ (e.g., $\kappa=10$ in our numerical study) most violated constraints in (26b) and (26d), i.e., find the epigraphs corresponding to the $\kappa$ largest values in the following set

Finally, we can generate and add valid inequalities by separating $(\widehat{\bm{x}},\widehat{\bm{w}},\widehat{\bm{z}},\widehat{\gamma},\widehat{\nu})$ from the convex hulls of these $\kappa$ epigraphs according to Lemma 3.

Numerical Demonstration

In this section, we present a series of numerical studies to demonstrate the effectiveness of the proposed formulations and also show how to use cross validation to choose a proper Wasserstein radius $\delta$.

where the chance constraint here is to guarantee that the worst-case probability that each knapsack’s capacity should be satisfied is at least $1-\epsilon$.

In the following subsections, we generated different random instances to test the proposed formulations. All the instances were executed on a MacBook Pro with a 2.80 GHz processor and 16GB RAM with a call of the commercial solver Gurobi (version 7.5, with default settings). We set the time limit of solving each instance to be 3600 seconds.

The numerical results with sample size $N=100$ are displayed in Table 2, where we use BigM Model, ${\bf{VaR}}$ Model, ${\bf{CVaR}}$ Model and ICCP Model denote exact formulation in Theorem 2, outer approximation in Corollary 4, ${\bf{CVaR}}$ approximation in Theorem 6 and inner chance constrained programming approximation in Corollary 5, respectively. We also use “Opt.Val” to denote the optimal value $v^{*}$, “Value” to denote the best objective value output from an approximation model and “Time” to denote the computational time in seconds. Additionally, since we can solve exact BigM Model to the optimality, we use GAP denote the optimality gap of an approximation model, which is computed as

We also let $\alpha^{*}$ denote the best $\alpha$ found in ICCP Model. In BigM Model (10), we chose a lower bound of $\nu$ as $\underline{\nu}=1$. We chosen the big-M coefficients in BigM Model, ${\bf{VaR}}$ Model, and ICCP Model according to the remarks after Theorem 2, Corollary 4, and Corollary 5, respectively, where $\bm{L}=\bm{0},\bm{U}=\bm{e}$. From Table 2, we see that all the models can be solved to the optimality within 2 minutes, where BigM Model and ICCP Model often take the longest time to solve, and for each instance, ${\bf{CVaR}}$ Model can be solved within a second. This might be because (i) ${\bf{CVaR}}$ Model is a second order conic program and does not involve any binary variables; (ii) on the contrary, the BigM Model not only has binary variables but also involves the most number of auxiliary variables, while to solve ICCP Model, one needs to solve $\lceil N\epsilon\rceil$ regular chance constrained programs. In terms of approximation accuracy, we see that VaR Model is usually 2-3% away from the true optimality, ${\bf{CVaR}}$ Model is 1-2% away from the true optimality, while ICCP Model nearly finds the true optimal solution. This demonstrates that all of the proposed approximation models can find near-optimal solutions.

The numerical results with sample size $N=1000$ are displayed in Table 3, where similarly, we use BigM Model, ${\bf{VaR}}$ Model, ${\bf{CVaR}}$ Model and ICCP Model denote exact formulation in Theorem 2, outer approximation in Corollary 4, ${\bf{CVaR}}$ approximation in Theorem 6 and inner chance constrained programming approximation in Corollary 5, respectively. We use “UB” to denote the best upper bound found by BigM Model or ${\bf{VaR}}$ Model, “LB” to denote the best lower bound found by BigM Model, ${\bf{CVaR}}$ Model, or ICCP Model, and “Time” to denote the computational time in seconds. Additionally, since we cannot solve the BigM Model to optimality, we use GAP denote its optimality gap, which is computed as

To evaluate the effectiveness of approximation models, we use Improvement to denote the percentage of differences between the bounds of approximation models and bounds of BigM Model, i.e., for the ${\bf{VaR}}$ Model,

while for the ${\bf{CVaR}}$ Model or ICCP Model,

where Approximation Model here is either ${\bf{CVaR}}$ Model or ICCP Model. We found that ICCP Model is difficult to solve these instances to optimality, and thus we chose a particular $\alpha=\frac{\epsilon}{2}$ in ICCP Model. Similarly, in BigM Model (10), we chose a lower bound of $\nu$ as $\underline{\nu}=1$, and the big-M coefficients in BigM Model, ${\bf{VaR}}$ Model, and ICCP Model were computed according to the remarks after Theorem 2, Corollary 4, and Corollary 5, respectively, where $\bm{L}=\bm{0},\bm{U}=\bm{e}$. From Table 3, we see that ${\bf{CVaR}}$ Model can be solved to optimality within 2 seconds, while all the other models cannot be solved within the time limit. In terms of approximation accuracy, we see that VaR Model consistently provides better upper bounds and can help close more 10% optimality gap on average compared to BigM Model, ${\bf{CVaR}}$ Model often provides slightly better feasible solutions than BigM Model, while, ICCP Model yields the best feasible solutions. This demonstrates that all of the proposed approximation models are useful, to some extent, to improve the exact bigM model. In particular, ${\bf{VaR}}$ Model provides a better upper bound, which helps evaluate the solution quality more accurately, ${\bf{CVaR}}$ Model and ICCP Model often provide better feasible solutions. Also, we notice that mixed integer programs ${\bf{VaR}}$ Model and ICCP Model outperform BigM Model, which might be because (i) the BigM Model requires more auxiliary variables than ${\bf{VaR}}$ Model or ICCP Model; (ii) the naive big-M coefficients of the BigM Model are typically larger than the other two. In practice, it is worthy of trying all the BigM Model, ICCP Model, and CVaR Model first, then choose the best solution from three models and use the outer approximation- VaR Model to provide a numerical optimality guarantee on how good the solution quality is.

2 Choosing a Wasserstein Radius using Cross Validation

for each $i\in[I]$ and $j\in[N]$, and $\rho\in\{0,0.1,\ldots,1\}$. For each $l\in[n]$, we independently generated $c_{l}$ from the uniform distribution on the interval $$, while for each$i\in[I]$, we set$b^{i}:=50$. We also suppose that the possible Wasserstein radii are from$\delta\in\{0.01,0.02,\ldots,0.1\}$.

For the comparison purpose, we also solve a regular chance constrained programming counterpart of DRMKP (30) with respect to the empirical samples $\{{\bm{\zeta}}^{j}\}_{j\in[N]}$ and used the same procedure to compute its 90-percentile violation. The numerical results are displayed in Table 4, where we use “CCP Model” to denote the chance constrained programming counterpart of DRMKP, and “Opt.Val” to denote the optimal value of a corresponding model.

3 Binary DRMKP: Strength of Big-M Free Formulation

and a lower bound of $\nu$ as $\underline{\nu}=1$. In the branch and cut implementation described in the end of Section 4, each time we added $\kappa=10$ EPI inequalities.

The results are displayed in Table 5. We use BigM Model and BigM-free Model to denote the big-M formulation in Theorem 2 and big-M free formulation in Corollary 7, respectively. In addition, we use UB, LB, GAP, Opt.Val and Time to denote the best upper bound, the best lower bound, optimality gap, the optimal objective value, and the total running time in seconds, respectively.

From Table 5, we observe that the overall running time of BigM-free Model significantly outperforms that of BigM Model, i.e., almost all of the instances of BigM-free Model can be solved within 10 minutes, while the majority of the instances of BigM Model reach the time limit. The main reasons are two-fold: (i) BigM Model involves $\mathcal{O}(N+n)$ binary variables and $\mathcal{O}(N\times I)$ continuous variables, while BigM-free Model only involves $\mathcal{O}(n)$ binary variables and $\mathcal{O}(N)$ continuous variables; and (ii) BigM Model contains big-M coefficients, while BigM-free Model does not. We also observe that, as the risk parameter $\epsilon$ increases or Wasserstein radius $\delta$ decreases, both formulations take longer time to solve, but BigM-free Model still significantly outperforms BigM Model. These results demonstrate the effectiveness of our proposed BigM-free Model.

Conclusion

In this paper, we studied a distributionally robust chance constrained problem (DRCCP) with Wasserstein ambiguity set. We showed that a DRCCP could be formulated as a conditional value-at-risk constrained optimization, thus admits tight inner and outer approximations. Once the feasible region is bounded, we showed that a DRCCP could be mixed integer representable with big-M coefficients and additional binary variables, i.e., a DRCCP can be formulated as a mixed integer conic program. We also compared various inner and outer approximations and proved their corresponding inclusive relations. We further proposed a big-M free formulation for a binary DRCCP and a branch and cut solution algorithm. The numerical studies demonstrated that the proposed formulations are quite promising.

Acknowledgments

The author would like to thank Professor Shabbir Ahmed (Georgia Tech) for his helpful comments on an earlier version of the paper. Valuable comments from the editors and three anonymous reviewers are gratefully acknowledged.

References