Towards large-scale probabilistic set covering problem: an efficient Benders decomposition approach

In this paper, we consider the probabilistic set covering problem (PSCP) of the form:

where $A$ is a 0-1 matrix of dimension $m\times n$, $c\in\mathbb{R}^{n}_{+}$ is the cost of columns, $\epsilon$ is a confidence parameter chosen by the decision maker, typically near zero, and $\xi$ is a random vector taking values in $\{0,1\}^{m}$. (PSCP) generalizes the well-known set covering problem (SCP) [29] where the random vector $\xi$ is replaced by the all ones vector and has obvious applications in location problems [8]. It also serves as a building block in various applications; see [4, 26].

(PSCP) falls in the class of probabilistic programs (also known as chance-constrained programs), which has been extensively investigated in the literatures. We refer to the surveys [3, 16, 15, 24] for a detailed discussion of probabilistic programs and the various algorithms. Beraldi and Ruszczyński [8] first investigated (PSCP) and developed a specialized branch-and-bound algorithm, which involves enumerating the so-called $(1-\epsilon)$-efficient points introduced by [23] and solving a deterministic SCP for each $(1-\epsilon)$-efficient point. Saxena et al. [26] studied the special case in which the random vector $\xi$ can be decomposed into $T$ blocks/subvectors $\xi^{1},\ldots,\xi^{T}$ formed by subsets $\mathcal{M}^{t}\subseteq[m]$, and $\xi^{t_{1}}$ and $\xi^{t_{2}}$ are independent for any two distinct $t_{1},t_{2}\in[T]$ (throughout the paper, for a nonnegative integer $\tau$, we denote $[\tau]=\{1,\ldots,\tau\}$ where $[\tau]=\varnothing$ if $\tau=0$). Problem (PSCP) in this special case can be rewritten as

where $A^{t}$ is the submatrix formed by the rows in $\mathcal{M}^{t}$, $t\in[T]$. Saxena et al. [26] developed an equivalent mixed integer programming (MIP) formulation for (PSCP’) based on the $(1-\epsilon)$-efficient and -inefficient points of the distribution of $\xi^{t}$ and solved it by a customized branch-and-cut (B&C) algorithm. For set covering models with uncertainty on the constraint matrix $A$, we refer to [2, 7, 12, 14, 20, 27, 28, 30] and the references therein.

One weakness of the above two approaches for solving (PSCP) is that it requires an enumeration of the $(1-\epsilon)$-efficient (and -inefficient) points of the distribution of $\xi$ or $\xi^{t}$, which could be computationally demanding, especially when the dimensions are large. In addition, the MIP formulation of Saxena et al. [26] may grow exponentially with the dimensions of $\xi^{t}$, making it only capable of solving (PSCP) with small dimensions of $\xi^{t}$. Luedtke and Ahmed [19] addressed this difficulty by using the sample average approximation (SAA) approach [19, 22], where (PSCP) is replaced by an approximation problem based on Monte Carlo samples of the random vector $\xi$. The approximation problem is also a form of (PSCP) but has a finite discrete distribution of $\xi$. As such, it can be reformulated as a deterministic MIP problem and solved by the state-of-the-art MIP solvers. In general, the selection of the scenario size is crucial for the approximation quality; the larger the scenario size, the smaller the approximation error (see [3, 19, 21] for related discussions). However, since the problem size of the MIP formulation grows linearly with the number of scenarios, this also poses a new challenge for MIP solvers to solve (PSCP) with a huge number of scenarios.

The first contribution of this paper is the proposed customized Benders decomposition (BD) algorithm that is capable of solving (PSCP) with a huge number of scenarios due to the following two reasons. First, the number of variables in the underlying Benders reformulation is independent of the number of scenarios of the random data. Second, the Benders cuts can be separated by an efficient combinatorial algorithm. Our results show that the proposed BD algorithm is much more efficient than a black-box MIP solver’s branch-and-cut and automatic BD algorithms.

The second contribution of this paper is that for (PSCP) with block structure of $\xi$ [26], i.e., problem (PSCP’), we develop a mixed integer nonlinear programming (MINLP) formulation for an alternative approximation problem that explicitly takes the block structure $\xi$ into consideration, and propose a BD algorithm for it. The proposed BD algorithm is much more robust than the state-of-the-art approach in [26] in terms of solving (PSCP’) with different input parameters. Moreover, we show that the new approximation problem can return almost the same solution as the traditional one (when the scenario size is large) but is much more computationally solvable by the proposed BD algorithm.

The rest of the paper is organized as follows. Section 2 presents the mathematical formulations for (PSCP) and (PSCP’), and provides a favorable property for the formulations. Section 3 develops the BD algorithms for solving the two formulations. Finally, Section 4 presents the numerical results.

Consider (PSCP) with (arbitrarily) finite discrete distributions of random vectors $\xi$, that is, there exist finitely many scenarios $\xi_{1},\ldots,\xi_{s}\in\{0,1\}^{m}$ such that $\mathbb{P}\{\xi=\xi_{i}\}=p_{i}$ for $i\in[s]$ and $\sum_{i=1}^{s}p_{i}=1$. We introduce for each $k\in[m]$, a binary variable $v_{k}$, where $v_{k}=1$ guarantees $A_{k}x\geq 1$, and for each $i\in[s]$, a binary variable $z_{i}$, where $z_{i}=1$ guarantees $v\geq\xi_{i}$. Then (PSCP) can be formulated as the following MIP formulation [19]:

where $\mathcal{S}_{k}=\left\{i\in[s]\,:\,\xi_{ik}=1\right\}$, $k\in[m]$. Constraints (1c)–(1e) ensure $\mathbb{P}\{v\geq\xi\}\geq 1-\epsilon$. Note that if $\xi_{i}=\boldsymbol{0}$ for some $i\in[s]$, we can set $z_{i}=1$ and remove variable $z_{i}$ from the formulation. Therefore, we can assume $\xi_{i}\neq\boldsymbol{0}$ for all $i\in[s]$ and thus $\cup_{k\in[m]}\mathcal{S}_{k}=[s]$ in the following.

For problem (PSCP’) (i.e., problem (PSCP) with block structure of $\xi$) with a finite discrete distribution of $\xi^{t}$: $\mathbb{P}\{\xi^{t}=\xi^{t}_{i}\}=p_{it}$ for $i\in[s]$ and $\sum_{i=1}^{s}p_{it}=1$, a similar convex MINLP problem can be presented as follows:

Here, $\mathcal{S}^{t}_{k}=\left\{i\in[s]\,:\,\xi^{t}_{ik}=1\right\}$ (similarly, we assume $\cup_{k\in\mathcal{M}^{t}}\mathcal{S}^{t}_{k}=[s]$); variable $\eta_{t}$ represents the value of $\ln(\mathbb{P}\{v^{t}\geq\xi^{t}\})$ where $v^{t}$ is subvector of $v$ formed by the rows in $\mathcal{M}^{t}$; and constraints (2c)–(2e) ensure that $\prod_{t=1}^{T}\mathbb{P}\{v^{t}\geq\xi^{t}\}\geq 1-\epsilon$.

Both problems (1) and (2) can be seen as approximations of problem (PSCP) in the context of the SAA approach. When the block number $T$ is equal to one, problems (1) and (2) are equivalent. However, when the block number $T$ is larger than one, the two problems are not equivalent in general as the finite discrete distribution of $\xi$, obtained from Monte Carlo sampling, may not have an independent block structure. In Section 4, we will perform computational experiments to compare the performance of the two approximations.

Although problems (1) and (2) can be solved to optimality by the state-of-the-art MIP/MINLP solvers, the potentially huge numbers of scenario variables $z$ and related constraints make this approach only able to solve moderate-sized instances. To alleviate the computational difficulty, the authors in [17, 18, 19] proposed some preprocessing techniques to reduce the numbers of variables and constraints in problem (1). The preprocessing techniques can also be extended to problem (2). Unfortunately, the numbers of variables and constraints after preprocessing can still grow linearly with the number of scenarios, making it still challenging to solve problems with a huge number of scenarios. For problem (2), Saxena et al. [26] proposed an equivalent MIP reformulation in the space of $x$ and $v$. However, the problem size of this MIP reformulation depends on the number of the so-called $(1-\epsilon)$-efficient and -inefficient points, which generally grows exponentially with the block sizes of $\xi^{t}$.

To overcome the above weakness, we will develop an efficient BD approach in Section 3, which projects the scenario variables $z$ out from problems (1) or (2) and adds the Benders cuts on the fly. This approach is based on a favorable property of the two formulations, detailed as follows. Let (MIP’) (respectively, (MINLP’)) be the relaxation of problem (1) (respectively, (2)) obtained by removing the integrality constraints $z\in\{0,1\}^{s}$ (respectively, $z\in\{0,1\}^{s\times T}$). The following proposition shows that this operation does not change the objective values of problems (1) and (2).

It is also easy to see that relaxing the binary variables $v$ into continuous variables in $[0,1]$ in formulations (MIP’) and (MINLP’) does not change their optimal values. However, as noted by [19], leaving the integrality constraints on variables $v$ can render MIP solvers to favor the branching on variables $v$ and effectively improve the overall performance. Our preliminary experiments also justified this statement (for both (MIP’) and (MINLP’)). Due to this, we decide to leave the integrality constraints $v\in\{0,1\}^{m}$ in the two formulations.

In this section, we first review the (generalized) BD approach for generic convex MINLP and then apply the BD approach to formulations (1) and (2).

In this section, we present computational results to illustrate the effectiveness of the proposed BD algorithms for solving large-scale PSCPs. The proposed BD algorithms were implemented in Julia 1.7.3 using CPLEX 20.1.0. We use a branch-and-Benders-cut (B&BC) approach [25] to implement the BD algorithms in Sections 3.2 and 3.3, where the Benders feasibility cuts are added on the fly at each node of the search tree via the cut callback framework. In addition, before starting the B&BC algorithm, we apply the processing technique [17, 18, 19] to simplify formulations (1) and (2) and follow [9, 11] to apply an iterative procedure to solve the LP relaxation of the master problem. The latter attempts to collect effective Benders feasibility cuts for the initialization of the relaxed master problem, and improve the overall performance. At every iteration of the procedure, the LP is solved, and the Benders feasibility cuts violated by the current LP solution are added to the LP. This loop is repeated until the dual bound increases less than or equal to 0.01% for 100 consecutive times. At the end of this procedure, all violated cuts are added as static cuts to the relaxed master problem. The parameters of CPLEX were set to run the code in a single-threaded mode, with a time limit of 7200 seconds and a relative MIP gap of 0%. The computational experiments were conducted on a cluster of computers equipped with Intel(R) Xeon(R) Gold 6140 CPU @ 2.30GHz. Throughout this section, all averages are reported as the shifted geometric mean, with a shift of $1$ [1].

Following [26], we construct PSCP instances using the 60 deterministic SCP instances [5].We construct instances with two block sizes, namely 10 and 20, and instances that do not have a block structure of the random variable $\xi$ (that is, $\xi$ has only a single block). We use two probability distributions of $\xi$ called circular and star distributions [8, 26] to construct scenarios of $\xi_{i}$ or $\xi_{i}^{t}$. For the circular distribution, the $k$-th entry of random vector $\xi$ is set to $\xi_{k}=\max\{Y_{k},Y_{(k\bmod m)+1}\}$, where $Y_{k}$, $k\in[m]$, are independent Bernoulli random variables following the distribution $\Pr\{Y_{k}=1\}=\alpha_{k}$, $\alpha_{k}$ are uniformly chosen from $[0.01,0.0275]$, and $\bmod$ denotes the modulo operator. For the star distribution, the $k$-th entry of random vector $\xi$ is set to $0$ if $V_{k}\leq 1$ and $1$ otherwise, where $V_{k}=Y_{k}+Y_{m+1}$, $Y_{k}$, $k\in[m+1]$, are independent Poisson random variables following the distribution $\Pr(Y_{k}=\tau)=\frac{e^{-\lambda_{k}}(\lambda_{k})^{\tau}}{\tau!}$ ($\tau\in\mathbb{Z}_{+}$), and $\lambda_{k}$ are uniformly chosen from $[0.1,0.2]$. The distributions of random variables $\xi^{t}$ can be described in a similar manner. In our random construction procedure, the scenario size $s$ and the confidence parameter $\epsilon$ are chosen from $\{10^{4},10^{6}\}$ and $\{0.05,0.1\}$, respectively. In total, there are $1440$ PSCP instances of problems (1) and (2) in our testbed.

We first compare the performance of problems (1) and (2) when using them as Monte Carlo sampling approximations to solve PSCPs. We use the BD algorithms in Sections 3.2 and 3.3 to solve the two problems. Table 1 reports the total number of instances solved within the time limit (S), the average CPU time in seconds (T), and the average separation time in seconds (ST). For instances that cannot be solved within the time limit, we report the average relative end gap (G(%)) in percentage returned by CPLEX. To see the difference between the two problems, we report the number of instances with different optimal values (ND) and the gap between the optimal values of these instances (defined by $\texttt{$\rm\Delta$O(\%)}=100\times\frac{|o_{1}-o_{2}|}{\max\{o_{1},o_{2}\}}$, where $o_{1}$ and $o_{2}$ represent the optimal values of the two problems). As shown in Table 1, for most cases, the optimal values of problems (1) and (2) are identical, and even if the optimal values are different in some cases, their difference is usually very small. This shows that the gap between the two problems is indeed very small, especially for the cases with a huge scenario size. As for the computational efficiency, for instances with multiple blocks (or block sizes of 10 and 20), problem (2) is much more solvable by the proposed BD algorithm. This is reasonable as by exploiting the block structure of the random variables $\xi$, up to $T$ cuts can be constructed in a single iteration, making the proposed BD algorithm for problem (2) converge much faster. In contrast, for instances with a single block, solving problem (2) is outperformed by solving problem (1).

Next, we compare the performance of the proposed BD algorithm (called BD) with the CPLEX’s default branch-and-cut and automatic BD algorithms (called CPX and AUTO-BD, respectively) for solving generic PSCPs with arbitrary distribution of $\xi$ (i.e., problem (1)). The results are summarized in Table 2. From Table 2, we see that the performance of BD is much better than CPX. Indeed, for large-scale instances with $10^{6}$ scenarios, CPX was even unable to solve the LP relaxation of problem (1) within the given time limit. Comparing AUTO-BD and BD, we observe that for small-scale instances with $10^{4}$ scenarios, BD is outperformed by AUTO-BD by a factor of $1.4$ (in terms of CPU time); but for large-scale instances with $10^{6}$ scenarios, BD outperforms AUTO-BD by a factor of $3.1$.

Finally, we compare the performance of the proposed BD algorithm with the B&C algorithm of Saxena et al. [26] (denoted by B&C) for solving PSCPs with block structure of $\xi$ (i.e., problem (2)). B&C solves an MIP formulation built on the so-called $(1-\epsilon)$-efficient and -inefficient points (which must be enumerated before starting the algorithm). Table 3 summarizes the performance results on instances with block sizes being $10$ and $20$. We do not report the results on instances with a single block, as the block size is too large to enumerate the $(1-\epsilon)$-efficient and -inefficient points within a reasonable time of period. For setting B&C, we report the average number of (linear) constraints that are used to model the probabilistic constraint (C); and for completeness, we also report the CPU time spent in the enumeration phase (ET), which is not included in the time spent in solving the MIP formulation. As shown in Table 3, when $\epsilon$ and the block size are small or $\xi^{t}$ is constructed by a circular distribution, the number of constraints in the MIP formulation of Saxena et al. [26] is small, and thus B&C performs better than BD. However, for instances with a large $\epsilon$, a large block size, or the star distribution of $\xi^{t}$, due to the large number of constraints in the MIP formulation, using B&C yields worse performance. In contrast, as the proposed BD algorithm does not need to solve a large MIP formulation, it performs much better on these instances. Moreover, it can avoid calling a potentially time-consuming enumeration procedure of $(1-\epsilon)$-efficient and -inefficient points.

In summary, our computational results illustrate the effectiveness of our proposed BD algorithm for solving PSCPs. More specifically, for generic PSCPs with arbitrary distribution of $\xi$, the proposed BD algorithm is much more efficient than the branch-and-cut and automatic BD algorithms of CPLEX; for PSCPs with block structure of $\xi$, the proposed BD algorithm is much more robust than the state-of-the-art approach in [26] in terms of solving (PSCP’) with different input parameters.