The set covering problem (SCP) is one of representative combinatorial optimization problems. We are given a set of elements , subsets () and their costs () for . We say that is a cover of if holds. The goal of SCP is to find a minimum cost cover of . The SCP is formulated as a 0-1 integer programming (0-1 IP) problem as follows:
where if holds and otherwise, and if and otherwise. That is, a column of matrix represents the corresponding subset by . For notational convenience, for each , let be the index set of subsets that contains the element .
The SCP is known to be NP-hard in the strong sense, and there is no polynomial time approximation scheme (PTAS) unless P = NP. However, the worst-case performance analysis does not necessarily reflect the experimental performance in practice. The continuous development of mathematical programming has much improved the performance of heuristic algorithms accompanied by advances in computing machinery [2, 3]. For example, Beasley  presented a number of greedy algorithms based on Lagrangian relaxation called the Lagrangian heuristics, and Caprara et al.  introduced pricing techniques into a Lagrangian heuristic algorithm to reduce the size of instances. Several efficient heuristic algorithms based on Lagrangian heuristics have been developed to solve very large-scale instances with up to 5000 constraints and 1,000,000 variables with deviation within about 1% from the optimum in a reasonable computing time [5, 6, 7, 8].
The SCP has important real applications such as crew scheduling , vehicle routing , facility location [10, 11], and logical analysis of data . However, it is often difficult to formulate problems in real applications as SCP, because they often have additional side constraints in practice. Most practitioners accordingly formulate them as general mixed integer programming (MIP) problems and apply general purpose solvers, which are usually less efficient compared with solvers specially tailored to SCP.
In this paper, we consider an extension of SCP introducing (i) multicover and (ii) generalized upper bound (GUB) constraints, which arise in many real applications of SCP such as vehicle routing [13, 14], crew scheduling , staff scheduling [16, 17] and logical analysis of data . The multicover constraint is a generalization of covering constraint [19, 20], in which each element must be covered at least ( is the set of nonnegative integers) times. The GUB constraint is defined as follows. We are given a partition of (, , ). For each block (), the number of selected subsets from the block (i.e., ) is constrained to be at most (). We call the resulting problem the set multicover problem with GUB constraints (SMCP-GUB), which is formulated as a 0-1 IP problem as follows:
This generalization of SCP substantially extends the variety of its applications. However, GUB constraints often make the pricing method less effective, because they often prevent solutions from containing highly evaluated variables together. To overcome this problem, we develop heuristic algorithms to reduce the size of instances, in which new evaluation schemes of variables are introduced taking account of GUB constraints. We also develop an efficient implementation of a 2-flip neighborhood local search algorithm that reduces the number of candidates in the neighborhood without sacrificing the solution quality. In order to guide the search to visit a wide variety of good solutions, we also introduce an evolutionary approach called the path relinking method  that generates new solutions by combining two or more solutions obtained so far.
The SMCP-GUB is NP-hard, and the (supposedly) simpler problem of judging the existence of a feasible solution is NP-complete, since the satisfiability (SAT) problem can be reduced to this decision problem. We accordingly allow the search to visit infeasible solutions violating multicover constraints and evaluate their quality by the following penalized objective function. Note that throughout the remainder of the paper, we do not consider solutions that violate the GUB constraints, and the search only visits solutions that satisfy the GUB constraints. Let (
is the set of nonnegative real values) be a penalty weight vector. A solutionis evaluated by
If the penalty weights are sufficiently large (e.g., holds for all ), then we can conclude SMCP-GUB to be infeasible when an optimal solution under the penalized objective function violates at least one multicover constraint. In our algorithm, the initial penalty weights () are set to for all . Starting from the initial penalty weight vector , the penalty weight vector is adaptively controlled to guide the search to visit better solutions.
We present the outline of the proposed algorithm for SMCP-GUB. The first set of initial solutions are generated by applying a randomized greedy algorithm several times. The algorithm then solves a Lagrangian dual problem to obtain a near optimal Lagrangian multiplier vector through a subgradient method (Section 2), which is applied only once in the entire algorithm. Then, the algorithm applies the following procedures in this order: (i) heuristic algorithms to reduce the size of instances (Section 5), (ii) a 2-flip neighborhood local search algorithm (Section 3), (iii) an adaptive control of penalty weights (Section 4), and (iv) a path relinking method to generate initial solutions (Section 6). These procedures are iteratively applied until a given time limit has run out.
2 Lagrangian relaxation and subgradient method
For a given vector , called a Lagrangian multiplier vector, we consider the following Lagrangian relaxation problem of SMCP-GUB:
We refer to as the Lagrangian cost associated with column . For any , gives a lower bound on the optimal value of SMCP-GUB (when it is feasible, i.e., there exists a feasible solution to SMCP-GUB).
The problem of finding a Lagrangian multiplier vector that maximizes is called the Lagrangian dual problem (LRD):
For a given , we can easily compute an optimal solution to as follows. For each block (), if the number of columns satisfying is equal to or less, then set for variables satisfying and for the other variables; otherwise, set for variables with the lowest Lagrangian costs and for the other variables.
The Lagrangian relaxation problem has integrality property. That is, an optimal solution to
is also optimal to its linear programming (LP) relaxation problem obtained by replacingin (4) with for all . In this case, any optimal solution to the dual of the LP relaxation problem of SMCP-GUB is also optimal to LRD, and the optimal value of the LP relaxation problem of SMCP-GUB is equal to .
A common approach to compute a near optimal Lagrangian multiplier vector is the subgradient method. It uses the subgradient vector , associated with a given , defined by
This method generates a sequence of nonnegative Lagrangian multiplier vectors , where is a given initial vector and is updated from by the following formula:
where is the best solution obtained so far under the penalized objective function with the initial penalty weight vector , and is a parameter called the step size.
When huge instances of SCP are solved, the computing time spent on the subgradient method becomes very large if a naive implementation is used. Caprara et al.  developed a variant of the pricing method on the subgradient method. They define a core problem consisting of a small subset of columns (), chosen among those having low Lagrangian costs (). Their algorithm iteratively updates the core problem in a similar fashion that is used for solving large-scale LP problems . In order to solve huge instances of SMCP-GUB, we also introduce a pricing method into the basic subgradient method (BSM) described, e.g., in .
3 2-flip neighborhood local search
The local search (LS) starts from an initial solution and repeats replacing with a better solution in its neighborhood until no better solution is found in . For a positive integer , the -flip neighborhood is defined by , where is the Hamming distance between and . In other words, is the set of solutions obtainable from by flipping at most variables. We first develop a 2-flip neighborhood local search algorithm (2-FNLS) as a basic component of the proposed algorithm. In order to improve efficiency, 2-FNLS searches first, and then .
We first describe the algorithm to search , called the 1-flip neighborhood search. Let
denote the increase of by flipping and , respectively, where and . The algorithm first searches for an improved solution obtainable by flipping by searching for satisfying and for the block containing , where . If an improved solution exists, it chooses with the minimum ; otherwise, it searches for an improved solution obtainable by flipping by searching for satisfying .
We next describe the algorithm to search , called 2-flip neighborhood search. Yagiura el al.  developed a 3-flip neighborhood local search algorithm for SCP. They derived conditions that reduce the number of candidates in and without sacrificing the solution quality. However, those conditions are not directly applicable to the 2-flip neighborhood search for SMCP-GUB because of GUB constraints. Below we derive the following three lemmas that reduce the number of candidates in by taking account of GUB constraints.
Let denote the increase of by flipping the values of and simultaneously.
Suppose that a solution is locally optimal with respect to . Then holds, only if .
See A. ∎∎
This lemma indicates that in searching for improved solutions in , it is not necessary to consider the simultaneous flip of variables and such that or . Based on this, we consider only the set of solutions obtainable by flipping and simultaneously. We assume that holds for the block containing or and are in the same block , because otherwise the move is infeasible. Let
denote the increase of in this case.
Suppose that a solution is locally optimal with respect to , and . Then holds, only if at least one of the following two conditions holds.
Both and belong to the same block satisfying .
See B. ∎∎
Suppose that a solution is locally optimal with respect to , and for a block and a pair of indices with and , , and hold. Let and . Then we have .
See C. ∎∎
Note that the condition of Lemma 3 implies that the condition (i) of Lemma 2 is satisfied. We can conclude that to find an improved solution satisfying the condition (i), it suffices to check only one pair for each block satisfying , instead of checking all pairs with , and (provided that the algorithm also checks the solutions satisfying the condition (ii) of Lemma 2).
The algorithm first searches for an improved solution satisfying the condition (i). For each block () satisfying , it checks the solution obtained by flipping and for and in with the minimum and , respectively. The algorithm then searches for an improved solution satisfying the condition (ii). Let denote the subset of obtainable by flipping . The algorithm searches for each in the ascending order of . If an improved solution is found, it chooses a pair and with the minimum among those in .
Algorithm 2-FNLS searches first, and then . The algorithm is formally described as follows.
A solution and a penalty weight vector .
A solution .
- Step 1:
If holds, choose with the minimum , set and return to Step 1.
- Step 2:
If holds, choose with the minimum , set and return to Step 2.
- Step 3:
For each block satisfying (), if holds for and with the minimum and (), respectively, set and . If the current solution has been updated at least once in Step 3, return to Step 3.
- Step 4:
For each in the ascending order of , if holds, choose with the minimum and set and . If the current solution has been updated at least once in Step 4, return to Step 1; otherwise output and exit.
We note that 2-FNLS does not necessarily output a locally optimal solution with respect to , because the solution is not necessarily locally optimal with respect to in Steps 3 and 4. Though it is easy to keep the solution locally optimal with respect to in Steps 3 and 4 by returning to Step 1 whenever an improved solution is obtained in Steps 2 or 3, we did not adopt this option because it consumes much computing time just to conclude that the current solution is locally optimal with respect to in most cases. We also note that the phase to search in the algorithm (i.e., Steps 1 and 2) always finishes with the search for an improved solution obtainable by flipping to prevent this phase from stopping at solutions having redundant columns.
Let one-round be the computation needed to find an improved solution in the neighborhood or to conclude that the current solution is locally optimal, including the time to update relevant data structures and/or memory [23, 24]. If implemented naively, 2-FNLS requires and one-round time for and , respectively, where . In order to improve computational efficiency, we keep the following auxiliary data
in memory to compute each and in time. We also keep the values of () in memory to update the values of and for in time when is changed, where (see D).
We first consider the one-round time for . In Steps 1 and 2, the algorithm finds and with the minimum and in time, respectively, by using the auxiliary data whose update requires time. Thus, the one-round time is reduced to for .
We next consider the one-round time for . In Step 3, the algorithm first finds and with the minimum and (), respectively, in time. The algorithm then evaluates in time by using (9), where . In Step 4, the algorithm first flips and temporarily updates the values of (), and (, ) in time so that the memory corresponding to these keeps the values of , and for the obtained from by flipping . Then, for searching , the algorithm evaluates
(in time for each pair of and ) only for each such that the value of has been changed to () during the temporary update. Note that the number of such candidates that satisfy is . When an improved solution was not found in , the updated memory values are restored in time to the original values , and before we try another candidate of . The time to search for each is therefore . Thus, the one-round time is reduced to for , where .
Because , , , , , and always hold, these orders are not worse than those of naive implementation, and they are much better if , and hold, which are the case for many instances. We also note that the computation time for updating the auxiliary data has little effect on the total computation time of 2-FNLS, because, in most cases, the number of solutions actually visited is much less than that of evaluated neighbor solutions.
4 Adaptive control of penalty weights
We observed that 2-FNLS tends to be attracted to locally optimal solutions of insufficient quality when the penalty weights are large. We accordingly incorporate a mechanism to adaptively control the values of () [8, 25, 26]; the algorithm iteratively applies 2-FNLS, updating the penalty weight vector after each call to 2-FNLS. We call such a sequence of calls to 2-FNLS the weighting local search (WLS) according to [27, 28].
Let denote the solution at which the previous 2-FNLS stops. Algorithm WLS resumes 2-FNLS from after updating the penalty weight vector . Starting from an initial penalty weight vector , where we set for all , the penalty weight vector is updated as follows. Let denote the best solution obtained in the current call to WLS with respect to the penalized objective function with the initial penalty weight vector . If holds, WLS uniformly decreases the penalty weights for all , where the parameter is decided so that for 15% of variables satisfying , the new value of becomes negative. Otherwise, WLS increases the penalty weights by
where is the amount of violation of the th multicover constraint, and is a parameter that is set to 0.2 in our computational experiments. Algorithm WLS iteratively applies 2-FNLS, updating the penalty weight vector after each call to 2-FNLS, until the best solution with respect to obtained in the current call to WLS has not improved in the last 50 iterations.
A solution .
A solution and the best solution with respect to obtained in the current call to WLS.
- Step 1:
Set , , and .
- Step 2:
Apply to obtain an improved solution and then set . Let be the best solution with respect to obtained during the call to .
- Step 3:
If holds, then set and ; otherwise, set . If holds, output and and halt.
- Step 4:
If holds, then uniformly decrease the penalty weights for all by ; otherwise, increase the penalty weights for all by (12). Return to Step 2.
5 Heuristic algorithms to reduce the size of instances
For a near optimal Lagrangian multiplier vector , the Lagrangian costs give reliable information on the overall utility of selecting columns for SCP. Based on this property, the Lagrangian costs are often utilized to solve huge instances of SCP. Similar to the pricing method for solving the Lagrangian dual problem, several heuristic algorithms successively solve a number of subproblems, also called core problems, consisting of a small subset of columns (), chosen among those having low Lagrangian costs () [5, 6, 7, 8]. The Lagrangian costs are unfortunately unreliable for selecting columns for SMCP-GUB, because GUB constraints often prevent solutions from containing more than variables with the lowest Lagrangian costs . To overcome this problem, we develop two evaluation schemes of columns for SMCP-GUB.
Before updating the core problem for every call to WLS, the algorithm heuristically fixes some variables to reflect the characteristics of the incumbent solution and the current solution . Let be a near optimal Lagrangian multiplier vector, and be an index set from which variables to be fixed are chosen. Let be the maximum value of the Lagrangian cost (). The algorithm randomly chooses a variable (
) with probability
. We note that the uniform distribution is used ifholds for all . The algorithm iteratively chooses and fixes a variable () until holds for 20% of multicover constraints . It then updates the Lagrangian multiplier if holds for , and computes the Lagrangian costs for , where is the index set of the fixed variables. The variable fixing procedure is formally described as follows.
The incumbent solution , the current solution and a near optimal Lagrangian multiplier vector .
A set of fixed variables and a Lagrangian multiplier vector .
- Step 1:
Set , , and .
- Step 2:
If holds, then set for each satisfying , output and , and halt.
- Step 3:
Randomly choose a column with probability defined by (13), and set . Return to Step 2.
Subsequent to the variable fixing procedure, the algorithm updates the instance to be considered by setting , (), and ().
The first evaluation scheme modifies the Lagrangian costs to reduce the number of redundant columns resulting from GUB constraints. For each block (), let be the value of the st lowest Lagrangian cost among those for columns in if holds and otherwise. We then define a score for , called the normalized Lagrangian score, by if holds, and otherwise.
The second evaluation scheme modifies the Lagrangian costs by replacing the Lagrangian multiplier vector with the adaptively controlled penalty weight vector . We define another score for , called the pseudo-Lagrangian score, by . Intuitive meaning of this score is that we consider a column to be promising if it covers many frequently violated constraints in the recent search. The variable fixing procedure for the second evaluation scheme is described in a similar fashion to that of the first evaluation scheme by replacing the Lagrangian multiplier vectors and with penalty weight vectors and , respectively.
Given a score vector