1 Introduction
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 01 integer programming (01 IP) problem as follows:
(1) 
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 NPhard in the strong sense, and there is no polynomial time approximation scheme (PTAS) unless P = NP. However, the worstcase 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 [4] presented a number of greedy algorithms based on Lagrangian relaxation called the Lagrangian heuristics, and Caprara et al. [5] 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 largescale 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 [5], vehicle routing [9], facility location [10, 11], and logical analysis of data [12]. 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 [15], staff scheduling [16, 17] and logical analysis of data [18]. 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 (SMCPGUB), which is formulated as a 01 IP problem as follows:
(2) 
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 2flip 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 [21] that generates new solutions by combining two or more solutions obtained so far.
The SMCPGUB is NPhard, and the (supposedly) simpler problem of judging the existence of a feasible solution is NPcomplete, 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 solution
is evaluated by(3) 
If the penalty weights are sufficiently large (e.g., holds for all ), then we can conclude SMCPGUB 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 SMCPGUB. 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 2flip 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 SMCPGUB:
(4) 
We refer to as the Lagrangian cost associated with column . For any , gives a lower bound on the optimal value of SMCPGUB (when it is feasible, i.e., there exists a feasible solution to SMCPGUB).
The problem of finding a Lagrangian multiplier vector that maximizes is called the Lagrangian dual problem (LRD):
(5) 
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 replacing
in (4) with for all . In this case, any optimal solution to the dual of the LP relaxation problem of SMCPGUB is also optimal to LRD, and the optimal value of the LP relaxation problem of SMCPGUB 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
(6) 
This method generates a sequence of nonnegative Lagrangian multiplier vectors , where is a given initial vector and is updated from by the following formula:
(7) 
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. [5] 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 largescale LP problems [22]. In order to solve huge instances of SMCPGUB, we also introduce a pricing method into the basic subgradient method (BSM) described, e.g., in [3].
3 2flip 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 2flip neighborhood local search algorithm (2FNLS) as a basic component of the proposed algorithm. In order to improve efficiency, 2FNLS searches first, and then .
We first describe the algorithm to search , called the 1flip neighborhood search. Let
(8) 
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 2flip neighborhood search. Yagiura el al. [8] developed a 3flip 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 2flip neighborhood search for SMCPGUB 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.
Lemma 1.
Suppose that a solution is locally optimal with respect to . Then holds, only if .
Proof.
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
(9) 
denote the increase of in this case.
Lemma 2.
Suppose that a solution is locally optimal with respect to , and . Then holds, only if at least one of the following two conditions holds.
 (i)

Both and belong to the same block satisfying .
 (ii)

.
Proof.
See B. ∎∎
Lemma 3.
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 .
Proof.
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 2FNLS searches first, and then . The algorithm is formally described as follows.
Algorithm 2FNLS
 Input:

A solution and a penalty weight vector .
 Output:

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 2FNLS 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 oneround 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, 2FNLS requires and oneround time for and , respectively, where . In order to improve computational efficiency, we keep the following auxiliary data
(10) 
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 oneround 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 oneround time is reduced to for .
We next consider the oneround 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
(11) 
(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 oneround 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 2FNLS, 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 2FNLS 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 2FNLS, updating the penalty weight vector after each call to 2FNLS. We call such a sequence of calls to 2FNLS the weighting local search (WLS) according to [27, 28].
Let denote the solution at which the previous 2FNLS stops. Algorithm WLS resumes 2FNLS 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
(12) 
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 2FNLS, updating the penalty weight vector after each call to 2FNLS, until the best solution with respect to obtained in the current call to WLS has not improved in the last 50 iterations.
Algorithm WLS
 Input:

A solution .
 Output:

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 SMCPGUB, 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 SMCPGUB.
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
(13) 
and fixes
. We note that the uniform distribution is used if
holds 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.Algorithm FIX
 Input:

The incumbent solution , the current solution and a near optimal Lagrangian multiplier vector .
 Output:

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 pseudoLagrangian 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