Alternate second order conic program reformulations for hub location under stochastic demand and congestion

Abstract

In this paper, we study the single allocation hub location problem with capacity selection in the presence of congestion at hubs. Accounting for congestion at hubs leads to a non-linear mixed integer program, for which we propose 18 alternate mixed integer second order conic program (MISOCP) based reformulations. Based on our computational studies, we identify the best MISOCP-based reformulation, which turns out to be 20–60 times faster than the state-of-the-art. Using the best MISOCP-based reformulation, we are able to exactly solve instances up to 50 nodes in less than half-an-hour. We also theoretically examine the dimensionality of the second order cones associated with different formulations, based on which their computational performances can be predicted. Our computational results corroborate our theoretical findings. Such insights can be helpful in the generation of efficient MISOCPs for similar classes of problems.

Appendices

Appendix

SHLPCC for the SK-based model

For SK-based model, the flow variables $x_{ikm}$ are replaced by path variables $x_{ijkm}$, where $x_{ijkm} = {\left\{ \begin{array}{ll} 1, &{}\text {if flows from }i\text { to }j\text { are routed via hub }k\text { and }m \\ 0, &{}\text {otherwise}. \end{array}\right. }$

Definition of other variables and parameters remain same.

1.1 Two-subscripted capacity allocation variable

$$\begin{aligned} \mathbf \small [SK-2s] \quad \text {min} \quad&{\sum _{i}\sum _{j}\sum _{k}\sum _{m}} F_{ijkm}x_{ijkm} + \sum _{k}\sum _{l} Q^l_{k}y_{kl} + \theta \sum _{k} 1/2 {{E}[N_{k}(y,z)]} \nonumber \\ \text {s.t.} \quad&(3){-}(5), (7){-}(10), (12), (14){-}(21)\nonumber \\&\sum _{m}x_{ijkm} =z_{ik} \qquad \qquad \forall i,j,k \end{aligned}$$

(76)

$$\begin{aligned}&\sum _{k}x_{ijkm} =z_{jm} \qquad \qquad \forall i,j,m \end{aligned}$$

(77)

$$\begin{aligned}&x_{ijkm} \in \{0,1\} \qquad \qquad \forall i,j,k,m,l \end{aligned}$$

(78)

Here, $F_{ijkm}=W_{ij}( \chi d_{ik}+\alpha d_{km}+ \delta d_{mj})$ is the total flow through path $i-j-k-m$. (76) and (77) connect the assignment variables and path variables.

[SK-MISOCP1] (3)–(5), (7), (9)–(10), (12), (14)-(21), (37)–(39), (76)–(78)

[SK-MISOCP2] (3)–(5), (7), (9)–(10), (12), (14)–(21), (41)–(44), (76)–(78).

[SK-MISOCP3] (3)–(5), (7), (9)–(10), (12), (14)–(21), (41)–(43), (46), (76)–(78).

[SK-MISOCP4] (3)–(5), (7), (9)–(10), (12), (14)–(21), (41)–(43), (47), (49), (50), (76)–(78).

[SK-MISOCP5] (3)–(5), (7)–(10), (12), (14)-(21), (52)–(54), (56), (76)–(78).

Table 4 Computation time for N = 25 for CAB dataset

1.2 Three-subscripted capacity allocation variable

$$\begin{aligned} \mathbf \small [SK-3s] \quad \text {min} \quad&{\sum _{i}\sum _{j}\sum _{k}\sum _{m}} F_{ijkm}x_{ijkm} + \sum _{k}\sum _{l} Q^l_{k}t^{l}_{kk} + \theta {{E}[N_{k}]} \nonumber \\ \text {s.t.} \quad&(26)-(28), (30), (32),(34) \nonumber \\&\sum _{m}x_{ijkm} =\sum _{l}t^{l}_{ik} \qquad \qquad \forall i,j,k \end{aligned}$$

(79)

$$\begin{aligned}&\sum _{k}x_{ijkm} =\sum _{l}t^{l}_{jm} \qquad \qquad \forall i,j,m \end{aligned}$$

(80)

$$\begin{aligned}&x_{ijkm} \in \{0,1\} \qquad \qquad \forall i,j,k,m,l \end{aligned}$$

(81)

[SK-MISOCP6] :

(26)–(28), (30), (32),(34), (58)–(60), (79)–(81).

[SK-MISOCP7] :

(26)–(28), (30), (32),(34), (62)–(64), (79)–(81).

[SK-MISOCP8] :

(26)–(28), (30), (32),(34),(66)–(70), (79)–(81).

[SK-MISOCP9] :

(26)–(28), (30), (32),(34), (71)–(73), (75), (79)–(81).

OA method

1.1 EK-OA-2s: OA-based method for EK-2s

The auxillary variable $L_{kl}$ and $\rho _{k}$ which were defined as $L_{kl} = \rho _{k} y_{kl}$ and $\rho _{k}=\frac{ s_{k}}{1+s{k}}$, imply

$$\begin{aligned} L_{kl} = {\left\{ \begin{array}{ll} 0, &{} \text { if }y_{kl} =0 \\ \frac{s_{k}}{1+s_{k}}, &{} \text { if }y_{kl}=1. \end{array}\right. } \end{aligned}$$

(82)

Also, earlier results, $\sum _{l} L_{kl} =\rho _{k} \quad \forall k$ , $L_{kl} \le y_{kl} \quad \forall k,l$ and ${\sum _{i}\sum _{j} W_{ij} z_{ik}} = \sum _{l} \gamma ^l_{k} L_{kl}$, remain. The function, $L_{kl} = \frac{s_{k}}{(1+s_{k})}$ is a concave function which can be approximated with piecewise linear functions that are tangent to the function $L_{kl}$ . The method chooses the minimum of these tangents at points $s^h_{k} \quad \forall h \in H , {k} \, \in {N} , \,{{l}} \in {L}$. The cuts are given by

$$\begin{aligned}&L_{kl} = \text {min}_{h \in H} \Bigg \{ \frac{1}{(1+s^h_{k})^2} s_{k} + \bigg ( \frac{s^h_{k}}{1+s^h_{k}} \bigg )^2 \Bigg \} \nonumber \\&\quad \Longleftrightarrow L_{kl} \le \frac{1}{(1+s^h_{k})^2} s_{k} + \bigg ( \frac{s^h_{k}}{1+s^h_{k}} \bigg )^2 \quad \forall k \in N , l \in L, h \in H \end{aligned}$$

(83)

In the OA-based method proposed by Elhedhli and Hu (2005), for every $k-l$ pair, the non linear congestion term is approximated with tangents (cuts), given by (82), at points $s_{k}^h$ (set at $h_{0}$ initially). At every iteration a relaxed mixed integer linear problem is solved. The solution of which not only gives the lower bound but also supplies information for the next cut. Also, this solution is feasible for the main problem thus giving the upper bound. The algorithm terminates when both upper and lower bounds are $\epsilon$ (or less) away from each other where $\epsilon \ge 0$. Günlük and Linderoth (2012) proposed perspective counterpart for (83) as

$$\begin{aligned} L_{kl} = \frac{s_{k}}{1+s_{k}/y_{kl}}{k} \, \in {N} , \,{{l}} \in {L} \end{aligned}$$

and the corresponding perspective cut at $s_{k}^h$ as

$$\begin{aligned} L_{kl} \le \frac{1}{(1+s^h_{k})^2} s_{k} + \bigg ( \frac{s^h_{k}}{1+s^h_{k}} \bigg )^2 y_{kl} , \quad \forall k, l \end{aligned}$$

(84)

The formulation for the EK-OA-2s is as follows:

$$\begin{aligned}&\mathbf [EK-OA-2s] \\&\text {min} \quad {\sum _{i}\sum _{k}} C_{ik}(\chi O_{i}+\delta D_{i})z_{ik}+\sum _{i}\sum _{k}\sum _{m} \alpha C_{km} x_{ikm} \\&\quad + \sum _{k}\sum _{l} Q^l_{k}y_{kl} + \theta /2 \sum _{k} \left\{ s_{k} + \rho _{k} +\sum _{l}c^2_{kl} \big ( V_{kl} -L_{kl} \big )\right\} \\&\quad \text {s.t.} \quad (3){-}(21), (84) \end{aligned}$$

Algorithm for the above discussed OA-based method is as follows:

1.2 EK-OA-3s: OA-based method for EK-3s

For formulations with $t^{l}_{ik}$, we had objective function as

$$\begin{aligned}&\text {min} \quad \left( {\sum _{i}\sum _{k}} C_{ik}(\chi O_{i}+\delta D_{i}) \sum _{l}t^l_{ik}\right) +\sum _{i}\sum _{k}\sum _{m} \alpha C_{km} x_{ikm}+ \sum _{k}\sum _{l} Q^l_{k}t^l_{kk} \\&\quad + \theta \sum _{k} \sum _{l} 1/2 \left\{ \left( 1+c^2_{kl} \right) \frac{\sum _{i}\sum _{j} W_{ij} t^l_{ik}}{\left( \gamma ^l_{k} - \sum _{i}\sum _{j} W_{ij} t^l_{ik}\right) }+ \left( 1-c^2_{kl} \right) \frac{\sum _{i}\sum _{j} W_{ij} t^l_{ik}}{\gamma ^l_{k}}\right\} \end{aligned}$$

We introduce variable $\rho _{kl}$ and $s_{kl}$ such that

$$\begin{aligned}&\frac{\sum _{i}\sum _{j} W_{ij} t^l_{ik}}{\gamma ^l_{k}} \le \rho _{kl} \qquad \quad \forall k,l \end{aligned}$$

(85)

$$\begin{aligned}&\frac{\sum _{i}\sum _{j} W_{ij} t^l_{ik}}{\gamma ^l_{k} - \sum _{i}\sum _{j} W_{ij} t^l_{ik}} \le s_{kl} \quad \forall k,l \end{aligned}$$

(86)

$\rho _{kl}$ and $s_{kl}$ are related as

$$\begin{aligned} \rho _{kl} = \frac{s_{kl}}{1+s_{kl}} \quad \forall k,l \end{aligned}$$

which is non-linear concave function and can be approximated using tangent hyperplanes as discussed in the previous section. For perspective reformulation, we introduce the following constraint

$$\begin{aligned} \rho _{kl} \le t^{l}_{kk} \end{aligned}$$

(87)

to the following equivalent perspective form

$$\begin{aligned} \rho _{kl} = \frac{s_{kl}}{1+(s_{kl}/t^{l}_{kk})} \end{aligned}$$

(88)

By following logic similar to (83), we have perspective cuts as

$$\begin{aligned} \rho _{kl} \le \displaystyle \frac{1}{(1+s^h_{kl})^2} s_{kl} + \displaystyle \frac{(s^h_{kl})^2}{(1+s^h_{kl})^2} t^{l}_{kk} , \quad \forall k,l \end{aligned}$$

(89)

The overall formulation for the approximation method is as follows:

$$\begin{aligned}&\mathbf [EK-OA-3s] \\&\text {min} {\sum _{i}\sum _{k}} C_{ik}(\chi O_{i}+\delta D_{i}) \sum _{l}t^l_{ik})+\sum _{i}\sum _{k}\sum _{m} \alpha C_{km} x_{ikm}+ \sum _{k}\sum _{l} Q^l_{k}t^l_{kk} \\&\quad + \theta \sum _{k} \sum _{l} 1/2 \Bigg \{ \Big (1+c^2_{kl} \Big ) s_{kl} + \Big (1-c^2_{kl} \Big ) \rho _{kl} \Bigg \}\\&\quad \text {s.t.} \quad (26)-(34), (87), (89) \end{aligned}$$

Algorithm for the above discussed OA-based method is as follows:

Table 5 Comparison of EK-MISOCPs against EK-OA-2s and EK-OA-3s for N = 10 (CAB dataset)

Table 6 Comparison of EK-MISOCPs against EK-OA-2s and EK-OA-3s for N = 15 (CAB dataset)

Table 7 Comparison of EK-MISOCPs against EK-OA-2s and EK-OA-3s for N = 20 (CAB dataset)

Table 8 Comparison of EK-MISOCPs against EK-OA-2s and EK-OA-3s for N = 25 (CAB dataset)

Performance profile for coefficient of variation (c), and unit congestion cost $\theta$ for CAB dataset

Fig. 2

Performance profile of EK-MISOCPs and EK-OA-based method for $\theta$= $\{20, 50\}$

Fig. 3

Performance profile of EK-MISOCPs and EK-OA-based method for $\{c= 0, 1, 2\}$

Table 9 Comparison among EK-MISOCP4, 5, 6 and 9 for $|\mathbf {N}|=\mathbf {25}$ (AP dataset)

Table 10 Comparison among EK-MISOCP4, 5, 6 and 9 for $|\mathbf {N}|=\mathbf {50}$ (AP dataset)