A distributionally robust chance-constrained model for humanitarian relief network design

Abstract

We propose a novel two-stage distributionally robust joint chance-constrained (DRJCC) model to design a resilient humanitarian relief network with uncertainties in demand and unit allocation cost of relief items in the post-disaster environment. This model determines the locations of the supply facilities with pre-positioning inventory levels and the transportation plans. We investigate the problem under two types of ambiguity sets: moment-based ambiguity and Wasserstein ambiguity. For moment-based ambiguity, we reformulate the problem into a mixed-integer conic program and solve it via a sequential optimization procedure by optimizing scaling parameters iteratively. For Wasserstein ambiguity, we reformulate the problem into a mixed-integer linear program. We conduct comprehensive numerical experiments to assess the computational efficiency of the proposed reformulation and algorithmic framework, and evaluate the reliability of the generated network by the proposed model. Through a case study in the Gulf Coast area, we demonstrate that the DRJCC model under Wasserstein ambiguity achieves a better trade-off between cost and network reliability in out-of-sample tests than the moment-based DRJCC model and the classical stochastic programming model.

Appendices

Appendix A: Proof of Theorem 1

Proof

Following Theorem 3 of (Zhang et al. 2016), let $\varvec{q}^*$ be the optimal solution for second-stage problem for given a feasible first-stage solution of decision variables $\varvec{y}$ and $\varvec{r}$. Denote $\Psi (\varvec{y},\varvec{r})$ as a feasible set of $\varvec{X}$ for given $\varvec{y},\varvec{r}$. Let set $\Lambda =\left\{ \varvec{x} \in {\mathbb {R}}^{|A|}: \;\; \begin{aligned} (\varvec{x}-\varvec{\mu }_{\zeta })^\intercal \varvec{\Sigma }_{\zeta }^{-1} (\varvec{x}-\varvec{\mu }_{\zeta }) \le \varphi ^2 \end{aligned} \right\} \subseteq {\mathbb {R}}^{|A|+}$. Then, we have

$$\begin{aligned} \sup _{F \in {\mathcal {F}}^M} {\mathbb {E}}_{F} \left[ \varvec{q}^\intercal \varvec{\zeta } \right] = \sup _{F \in {\mathcal {F}}^M} \; \min _{\varvec{q} \in \Psi (\varvec{y},\varvec{r})} \varvec{q}^\intercal {\mathbb {E}}_{F} [\varvec{\zeta }] = \max _{{\mathbb {E}} \left[ \varvec{\zeta }\right] \in \Lambda } \; \min _{\varvec{q} \in \Psi (\varvec{y},\varvec{r})} \varvec{q}^\intercal {\mathbb {E}}_{F} [\varvec{\zeta }]. \end{aligned}$$

(10a)

Note that $\min _{\varvec{q} \in \Psi (\varvec{y},\varvec{r})} \varvec{q}^\intercal {\mathbb {E}}_{F} [\varvec{\zeta }] = {(\varvec{q}^*)}^\intercal {\mathbb {E}}_{F} [\varvec{\zeta }]$, by taking the maximum over ${\mathbb {E}}_{F} [\varvec{\zeta }] \in \Lambda$ on both sides of the equation, we have

$$\begin{aligned} \max _{{\mathbb {E}} \left[ \varvec{\zeta }\right] \in \Lambda } \; \min _{\varvec{q} \in \Psi (\varvec{y},\varvec{r})} \varvec{q}^\intercal {\mathbb {E}}_{F} [\varvec{\zeta }] =\max _{{\mathbb {E}} \left[ \varvec{\zeta }\right] \in \Lambda } \; {(\varvec{q}^*)}^\intercal {\mathbb {E}}_{F} [\varvec{\zeta }]. \end{aligned}$$

(10b)

We know that $\varvec{q} \in \Psi (\varvec{y},\varvec{r})$, the next inequality follows:

$$\begin{aligned} \max _{{\mathbb {E}} \left[ \varvec{\zeta }\right] \in \Lambda } \; {(\varvec{q}^*)}^\intercal {\mathbb {E}}_{F} [\varvec{\zeta }] \ge \min _{\varvec{q} \in \Psi (\varvec{y},\varvec{r})} \max _{{\mathbb {E}} \left[ \varvec{\zeta }\right] \in \Lambda } \; \varvec{q}^\intercal {\mathbb {E}}_{F} [\varvec{\zeta }] \end{aligned}$$

(10c)

Combining (10b) and (10c), we can get

$$\begin{aligned} \max _{{\mathbb {E}} \left[ \varvec{\zeta }\right] \in \Lambda } \; \min _{\varvec{q} \in \Psi (\varvec{y},\varvec{r})} \varvec{q}^\intercal {\mathbb {E}}_{F} [\varvec{\zeta }] \ge \min _{\varvec{q} \in \Psi (\varvec{y},\varvec{r})} \max _{{\mathbb {E}} \left[ \varvec{\zeta }\right] \in \Lambda } \; \varvec{q}^\intercal {\mathbb {E}}_{F} [\varvec{\zeta }] \end{aligned}$$

(10d)

According to minimax weak duality (Bertsekas 2009), we can easily get

$$\begin{aligned} \max _{{\mathbb {E}} \left[ \varvec{\zeta }\right] \in \Lambda } \; \min _{\varvec{q} \in \Psi (\varvec{y},\varvec{r})} \varvec{q}^\intercal {\mathbb {E}}_{F} [\varvec{\zeta }] \le \min _{\varvec{q} \in \Psi (\varvec{y},\varvec{r})} \max _{{\mathbb {E}} \left[ \varvec{\zeta }\right] \in \Lambda } \; \varvec{q}^\intercal {\mathbb {E}}_{F} [\varvec{\zeta }] \end{aligned}$$

(10e)

Combining (10d) and (10e), we have

$$\begin{aligned} \max _{{\mathbb {E}} \left[ \varvec{\zeta }\right] \in \Lambda } \; \min _{\varvec{q} \in \Psi (\varvec{y},\varvec{r})} \varvec{q}^\intercal {\mathbb {E}}_{F} [\varvec{\zeta }] = \min _{\varvec{q} \in \Psi (\varvec{y},\varvec{r})} \max _{{\mathbb {E}} \left[ \varvec{\zeta }\right] \in \Lambda } \; \varvec{q}^\intercal {\mathbb {E}}_{F} [\varvec{\zeta }] \end{aligned}$$

(10f)

Because the optimal objective value of $\max _{{\mathbb {E}} \left[ \varvec{\zeta }\right] \in \Lambda } \varvec{q}^\intercal {\mathbb {E}}_{F} [\varvec{\zeta }]$ is $\varphi \sqrt{\varvec{q}^\intercal \varvec{\Sigma }_{\zeta }\varvec{q}} + \varvec{\mu }_{\zeta }^\intercal \varvec{q}$, we can get the reformulation in Theorem 1. $\square$

Appendix B: Proof of Theorem 2

Proof

Denote $g_i = r_{i} + \sum _{(j,i)\in A} x_{ji} - \sum _{(i,j)\in A} x_{ij}$. According to the definition of the worst-case CVaR, constraint (4) can be expressed as follows:

$$\begin{aligned}&\sup _{G \in {\mathcal {G}}^M}\textrm{CVaR}_{\epsilon } \left( \max _{i = 1...|I|} \left\{ \alpha _i \left( \xi _i - g_i \right) \right\} \right) \nonumber \\ =&\inf _{\beta \in {\mathbb {R}}} \left\{ \beta + \frac{1}{\epsilon } \sup _{G \in {\mathcal {G}}^M} {\mathbb {E}}_G \left( \left[ \max _{i = 1...|I|} \left\{ \alpha _i \left( \xi _i - g_i \right) \right\} -\beta \right] ^+ \right) \right\} . \end{aligned}$$

(11a)

To solve the worst-case CVaR problem, the following worst-case expectation (11b) should be solved at the beginning.

$$\begin{aligned} \sup _{G \in {\mathcal {G}}^M}{\mathbb {E}}_G \left( \left[ \max _{i = 1...|I|} \alpha _i\left( \xi _i - g_i \right) -\beta \right] ^+ \right) \end{aligned}$$

(11b)

For fixed first-stage decision variables $r_i, \; x_{ij} \in {\mathbb {R}}^+$ and $\varvec{\alpha } \in {\mathcal {A}} = \{\varvec{\alpha } \in \varvec{R}^{|I|}:\varvec{\alpha }>\varvec{0} \}$, the worst-case expectation can be expressed as the following problem (12)

$$\begin{aligned}&\sup _G \int \left( \left[ \max _{i = 1...|I|} \alpha _i\left( \xi _i - g_i \right) -\beta \right] ^+ \right) \mathrm{{d}} G(\varvec{\xi }) \end{aligned}$$

(12a)

$$\begin{aligned}&\mathrm{s.t.} \int \mathrm{{d}} G(\varvec{\xi }) = 1 \end{aligned}$$

(12b)

$$\begin{aligned}&\quad \; \int \varvec{\xi } \mathrm{{d}} G(\varvec{\xi }) = \varvec{\mu }_{\zeta } \end{aligned}$$

(12c)

$$\begin{aligned}&\quad \; \int \varvec{\xi }\varvec{\xi }^\intercal G F(\varvec{\xi }) = \varvec{\Sigma }_{\zeta } + \varvec{\mu }_{\zeta } \varvec{\mu }_{\zeta }^\intercal \end{aligned}$$

(12d)

where F is probability measure on ${\mathbb {R}}^{|I|}$. Assign dual variables $\delta _1 \in {\mathbb {R}},\varvec{\delta }_2 \in {\mathbb {R}}^{|I|}, \varvec{\Delta }_3 \in {\mathbb {Q}}^{|I|}$ to constraints (12b), (12c), and (12d), respectively. Then, we have the following dual problem (13)

$$\begin{aligned} \inf \quad&\delta _1 + \varvec{\delta }_2^\intercal \varvec{\mu }_{\zeta } + \langle \varvec{\Delta }_3, \varvec{\Sigma }_{\zeta } + \varvec{\mu }_{\zeta } \varvec{\mu }_{\zeta }^\intercal \rangle \end{aligned}$$

(13a)

$$\begin{aligned} \mathrm{s.t.} \quad&\delta _1 \in {\mathbb {R}},\varvec{\delta }_2 \in {\mathbb {R}}^{|I|}, \varvec{\Delta }_3 \in {\mathbb {R}}^{{|I| \times |I|}} \end{aligned}$$

(13b)

$$\begin{aligned}&\delta _1 + \varvec{\delta }_2^\intercal \varvec{\xi } + \langle \varvec{\Delta }_3, \varvec{\xi } \varvec{\xi }^\intercal \rangle \ge \left[ \max _{i = 1...|I|} \alpha _i\left( \xi _i - g_i \right) -\beta \right] ^+, \quad \forall \; \varvec{\xi } \in {\mathbb {R}}^{|I|} \end{aligned}$$

(13c)

Because $\varvec{\Sigma }_{\zeta } \succ 0$, the strong duality holds between problems (12) and (13) (Isii et al. 1960). We assign an indicate vector $\varvec{1}_i \in {\mathbb {R}}^{|I|}$ to represent a vector that have the ith element is 1, and the other elements are zero. The indicate vector $\varvec{1}_i$ is used to represent the right-hand side of constraint (13c). Then, we have:

$$\begin{aligned} \delta _1 + \varvec{\delta }_2^\intercal \varvec{\xi } + \langle \varvec{\Delta }_3, \varvec{\xi } \varvec{\xi }^\intercal \rangle \ge \alpha _i \left( \varvec{\xi }^\intercal \varvec{1}_i - g_i \right) -\beta , \quad \quad \forall \; \varvec{\xi } \in {\mathbb {R}}^{|I|}, \; i \in I \end{aligned}$$

(14)

Denoting a combined decision variable $\varvec{S} = \begin{bmatrix} &{}\varvec{\Delta }_3 \quad &{} \quad \varvec{\delta }_2 \\ &{}\varvec{\delta }_2 \quad &{} \quad \delta _1 \end{bmatrix}$, constraint (14) can be reformulated as:

$$\begin{aligned} \varvec{S} - \begin{bmatrix} \varvec{0} \;&{}\; \frac{1}{2} \alpha _i\varvec{1}_i \\ \frac{1}{2} \alpha _i\varvec{1}_i^\intercal \;&{}\; -\alpha _i g_i -\beta \end{bmatrix} \succeq 0, \quad i \in I \end{aligned}$$

(15)

where $\varvec{0} \in {\mathbb {Q}}^{|I|}$ is a all zero matrix. By defining $\varvec{\Omega } = \begin{bmatrix} \varvec{\Sigma }_{\zeta }+\varvec{\mu }_{\zeta } \varvec{\mu }_{\zeta }^\intercal \;&{}\; \varvec{\mu }_{\zeta } \\ \varvec{\mu }_{\zeta }^\intercal \;&{}\; 1 \end{bmatrix}$, we have the tractable reformulation for problem (11a). Therefore, Theorem 2 is proved. $\square$

Appendix C: Proof of Theorem 3

Proof

Because $f(\varvec{y},\varvec{r},\varvec{\zeta },\varvec{\xi }) > - \infty$, according to Proposition 5.5.1 proposed by Bertsekas (2009), we have $\sup _{F \in {\mathcal {F}}} \min _{\varvec{x}} {\mathbb {E}}_{F} [f(\varvec{y},\varvec{r},\varvec{\zeta },\varvec{\xi }) ] = \min _{\varvec{x}} \sup _{F \in {\mathcal {F}}} {\mathbb {E}}_{F} [f(\varvec{y},\varvec{r},\varvec{\zeta },\varvec{\xi }) ]$. For a feasible second-stage auxiliary decision vector $\bar{\varvec{q}}$, the worst-case expectation $\sup _{F \in {\mathcal {F}}} {\mathbb {E}}_{F} [ h(\bar{\varvec{q}},\varvec{\zeta }) ]$ can be acquired by solving following problem:

$$\begin{aligned}&\sup \int _{\Xi } \sum _{m \in M} h (\bar{\varvec{q}},\varvec{\zeta }) Q(\mathrm{{d}} \varvec{\zeta }, \hat{\varvec{\zeta }}_m) \end{aligned}$$

(16a)

$$\begin{aligned}&\mathrm{s.t.} \; \int _{\Xi } Q(\mathrm{{d}} \varvec{\zeta }, \hat{\varvec{\zeta }}_m) = \frac{1}{M} ,\; \forall m \in M \end{aligned}$$

(16b)

$$\begin{aligned}&\quad \; \; \int _{\Xi } \sum _{m \in M} \mathrm{{d}}(\varvec{\zeta }, \hat{\varvec{\zeta }}_m) Q (\mathrm{{d}} \varvec{\zeta }, \hat{\varvec{\zeta }}_m) \le \theta _F \end{aligned}$$

(16c)

Developing Lagrangian function for problem (16), we get

$$\begin{aligned}&L(\varvec{\zeta }, \lambda , \varvec{s}) = \int _{\Xi } \sum _{m \in M} h (\bar{\varvec{q}},\varvec{\zeta }) Q(\mathrm{{d}} \zeta , \hat{\varvec{\zeta }}_m) - \int _{\Xi } \sum _{k=1}^N s_m Q(\mathrm{{d}} \varvec{\zeta }, \hat{\varvec{\zeta }}_m) \nonumber \\&\quad - \int _{\Xi } \sum _{m \in M} \lambda \mathrm{{d}}(\varvec{\zeta }, \hat{\varvec{\zeta }}_m) Q (\mathrm{{d}} \varvec{\zeta }, \hat{\varvec{\zeta }}_m) \nonumber \\&\quad + \frac{1}{M} \sum _{m \in M} s_m + \lambda \theta _F \end{aligned}$$

(17a)

Then, we have the dual problem of problem (16) as follows:

$$\begin{aligned}&\min _{\varvec{s}, \lambda } \; \frac{1}{M} \sum _{m \in M} s_m + \lambda \theta _F \end{aligned}$$

(17b)

$$\begin{aligned}&\; \mathrm{s.t.} \;\;\; h (\bar{\varvec{q}},\varvec{\zeta }) - \lambda \mathrm{{d}}(\varvec{\zeta }, \hat{\varvec{\zeta }}_m) \le s_m ,\; \forall \varvec{\zeta } \in \Xi , m \in M \end{aligned}$$

(17c)

$$\begin{aligned}&\quad \quad \; \lambda \ge 0 \end{aligned}$$

(17d)

Constraint (17c) can be expressed as

$$\begin{aligned} \sup _{\varvec{\zeta } \in \Xi _F} \left\{ \varvec{\zeta }^\intercal \bar{\varvec{q}} - \lambda \Vert \varvec{\zeta }-\hat{\varvec{\zeta }}_m\Vert _p \right\} \le s_m ,\; \forall m \in M \end{aligned}$$

(17e)

Because the allocation cost cannot be negative, $\Xi _F = [0, \infty ]$. Let $\Delta \nu _m = \varvec{\zeta }-\hat{\varvec{\zeta }}_m$, then we have the Lagrangian function of (17e) as follows:

$$\begin{aligned} L (\Delta \nu _m, \varvec{\gamma }_m ) = ( \bar{\varvec{q}} + \varvec{\gamma }_m)^\intercal (\Delta \nu _m + \varvec{\zeta }_m) - \lambda \Vert \Delta \nu _m \Vert _p \end{aligned}$$

(17f)

Then, we have the Lagrangian dual function of (17f)

$$\begin{aligned} \sup _{\Delta \nu _m} L (\Delta \nu _m, \varvec{\gamma }_m )&= \sup _{\Delta \nu _m} ( \bar{\varvec{q}} + \varvec{\gamma }_m)^\intercal (\Delta \nu _m + \varvec{\zeta }_m) - \lambda \Vert \Delta \nu _m \Vert _p \end{aligned}$$

(17g)

$$\begin{aligned}&= \sup _{\Delta \nu _m} \left\{ (\Vert \bar{\varvec{q}} + \varvec{\gamma }_m \Vert _p^* - \lambda ) \Vert \Delta \nu _m \Vert _p \right\} \end{aligned}$$

(17h)

where $\Vert \cdot \Vert _p^*$ is the dual norm of $l_p$ norm. Equality (17h) holds because of Lemma 1 in Zhang et al. (2017). If $\Vert \bar{\varvec{q}} + \varvec{\gamma }_m \Vert _p^* - \lambda > 0$, the problem $\sup _{\Delta \nu _m} L (\Delta \nu _m, \varvec{\gamma }_m )$ is unbounded. If $\Vert \bar{\varvec{q}} + \varvec{\gamma }_m \Vert _p^* - \lambda \le 0$, the problem $\sup _{\Delta \nu _m} L (\Delta \nu _m, \varvec{\gamma }_m )$ is equivalent with the following:

$$\begin{aligned}&\min _{\varvec{\gamma }_m} \; ( \bar{\varvec{q}} + \varvec{\gamma }_m)^\intercal \varvec{\zeta }_m \end{aligned}$$

(18a)

$$\begin{aligned}&\; \mathrm{s.t.} \;\;\; \Vert \bar{\varvec{q}} + \varvec{\gamma }_m \Vert _p^* \le \lambda \end{aligned}$$

(18b)

$$\begin{aligned}&\quad \quad \; \varvec{\gamma }_m \ge 0 \end{aligned}$$

(18c)

Bringing problem (18) back to (17c), we have a tractable reformulation of the problem DR-HRN. $\square$

Appendix D: Pseudo-code of Algorithm 1

Appendix E: Proof of Theorem 4

Proof

Based on Theorem 3 of Chen et al. (2022), the chance constraint (2b) can be reformulated as

$$\begin{aligned}&t \ge 0, \varvec{e} \ge 0 \end{aligned}$$

(19a)

$$\begin{aligned}&\epsilon t \ge \theta _G + \frac{1}{N} \sum _{n \in N} e_n \end{aligned}$$

(19b)

$$\begin{aligned}&\textbf{dist}(\varvec{\xi }_n, \varvec{g}) \ge t - e_n, \quad \forall n \in N \end{aligned}$$

(19c)

where $\textbf{dist}(\hat{\varvec{\xi }}_n, \varvec{g})$ is the distance between data point $\hat{\varvec{\xi }}_n$ of the empirical distribution ${\widehat{G}}$ and the lower bound of feasible region. Therefore, we have

$$\begin{aligned} \textbf{dist}(\varvec{\xi }_n, \varvec{g}) = \max \left\{ 0, \min _{i \in I} \{g_i - \xi _i\} \right\} \end{aligned}$$

(19d)

We denote the binary variable $z_n$. Based on Proposition 2 of Chen et al. (2022), we have the reformulation for constraint (19c) as follows:

$$\begin{aligned}&\varvec{z} \in \{ 0,1 \}^{N} \end{aligned}$$

(19e)

$$\begin{aligned}&{\widetilde{M}} (1-\varvec{z}_n) \ge t - e_n, \quad \forall n \in N \end{aligned}$$

(19f)

$$\begin{aligned}&g_i -\xi _{in} + {\widetilde{M}} z_n \ge t - e_n, \quad \forall i \in I ,\; n \in N \end{aligned}$$

(19g)

Then, according to Theorem 3 of Ho-Nguyen et al. (2021), the following reformulation is developed for constraint (19g) to improve the computational efficiency.

$$\begin{aligned}&\epsilon t \ge \theta _G + \frac{1}{N} \sum _{n \in N} e_n \end{aligned}$$

(19h)

$$\begin{aligned}&\sum _{n \in N} z_n \le \lfloor \epsilon N \rfloor \end{aligned}$$

(19i)

$$\begin{aligned}&g_i - q_i \ge t, \;&\forall i \in I \end{aligned}$$

(19j)

$$\begin{aligned}&g_i -\xi _{in} + (\xi _{in} - q_i) z_n \ge t - e_n,\;&\forall n \in N_{i}^*, \; \forall i \in I \end{aligned}$$

(19k)

where $N_{i}^*:= \{ n \in N: \xi _{in} > q_i \},\; \forall i \in I$ and $q_i$ equal to the $\lfloor \epsilon N \rfloor$th largest value of $\xi _{i}$. $\square$

Appendix F: Stochastic Programming Model

The uncertain parameter of scenario-based approach is given by a set of scenarios. The random variables $\varvec{\zeta }$ and $\varvec{\xi }$ are assumed to follow a discrete distribution independently. Let $\zeta _{ijw}$ be the observation of the unit allocation cost between location i and j in scenario w and $\xi _{iw}$ denote the demand at location i under the scenario w. We also assume that the probability of each scenario equal.

We reformulate the JCC to linear inequalities by defining a binary decision variable $b_w$ and a large number M. The objective function of the SP model is formulated as first-stage cost plus the expected second-stage cost over all scenarios. Denote the scenario set W indexed by w. The model shows the following:

$$\begin{aligned} \min \quad&\sum _{i \in I} C^F y_{i} +\sum _{i\in I} C^P r_{i} +\frac{1}{|W|} \sum _{w \in W} \left[ \sum _{(i,j)\in A} D^T_{ij} x_{ijw} \zeta _{ijw} \right] \end{aligned}$$

(20a)

$$\begin{aligned} \text {s.t.} \quad&O^S r_{i} \le C^H y_{i}, \quad&\forall i \in I \end{aligned}$$

(20b)

$$\begin{aligned}&\sum _{w \in W} b_w \le |W| (1-\epsilon ) \end{aligned}$$

(20c)

$$\begin{aligned}&r_{i} + \sum _{(j,i)\in A} x_{jiw} - \sum _{(i,j)\in A} x_{ijw} \ge \xi _{iw} - M (1-b_w) , \quad&\forall i \in I,\quad \forall w \in W \end{aligned}$$

(20d)

$$\begin{aligned}&y_{i} \in \{0,1\} , \quad&\forall i \in I \end{aligned}$$

(20e)

$$\begin{aligned}&r_{i} \ge 0, \quad&\forall i \in I \end{aligned}$$

(20f)

$$\begin{aligned}&b_w \in \{0,1\} , \quad&\forall w \in W \end{aligned}$$

(20g)

Constraints (20c) and (20d) ensure that the probability of the demand satisfaction for all scenarios is greater than $1-\epsilon$. If $b_w=1$, constraint (20d) equal to $r_{i} + \sum _{(j,i)\in A} x_{jiw} - \sum _{(i,j)\in A} x_{ijw} \ge \xi _{iw}$, if $b_w=0$, constraint (20d) is redundant.