A Branch-Price-and-Cut Procedure for the Discrete Ordered Median Problem

INFORMS JOURNAL ON COMPUTING Articles in Advance, pp. 1–18 http://pubsonline.informs.org/journal/ijoc ISSN 1091-9856 (print), ISSN 1526-5528 (online) A Branch-Price-and-Cut Procedure for the Discrete Ordered Median Problem Samuel Deleplanque, a Martine Labbé, b Diego Ponce, c,d,e Justo Puerto c a Ifsttar, COSYS, ESTAS, Université Lille Nord de France, 59000 Lille, France; b Départament d’Informatique, Faculté des Sciences, Université Libre de Bruxelles, 1050 Bruxelles, Belgium; c Instituto de Matem´ aticas de la Universidad de Sevilla (IMUS), 41012 Sevilla, Spain; d Department of Mechanical, Industrial and Aerospace Engineering, Concordia University, Montreal, Quebec H3G 1M8, Canada; e Centre interuniversitaire de recherche sur les réseaux d’enterprise, la logistique et le transport (CIRRELT), Montreal, Quebec H3T 1J4, Canada Contact: [email protected],http://orcid.org/0000-0003-4119-6006 (SD); [email protected], http://orcid.org/0000-0001-7471-2308 (ML); [email protected],http://orcid.org/0000-0003-3380-6601 (DP); [email protected], http://orcid.org/0000-0003-4079-8419 (JP) Received: February 9, 2018 Revised: December 7, 2018; March 11, 2019 Accepted: June 2, 2019 Published Online in Articles in Advance: January 7, 2020 https://doi.org/10.1287/ijoc.2019.0915 Copyright: © 2020 INFORMS Abstract. The discrete ordered median problem (DOMP) is formulated as a set-partitioning problem using an exponential number of variables. Each variable corresponds to a set of demand points allocated to the same facility with the information of the sorting position of their corresponding costs. We develop a column generation approach to solve the continuous relaxation of this model. Then we apply a branch-price-and-cut algorithm to solve smallto large-sized instances of DOMP in competitive computational time. Funding: This research was supported by Spanish/FEDER [Grant MTM2016-74983-C02-01-R]. The research of the second and third authors was supported by the Interuniversity Attraction Poles Programme initiated by the Belgian Science Policy Office. History: Accepted by Andrea Lodi, Area Editor for Design & Analysis of Algorithms—Discrete. Keywords: discrete optimization •location theory •branch and price •ordered median problems 1. Introduction Broadly speaking, a facility-location problem consists of locating one or several facilities to minimize an objective function of the assignment costs of clients to facilities. Median and center location problems constitute the most popular ones. The ordered median location problem is a flexible model that provides a common framework to cast most popular location problems. The objectivefunctiontobeminimizedisaweightedfunction of the allocation costs in which the weights are assigned to the ordered values of the costs rather than to specific costs. Ordered median location problems were first introduced in networks and continuous spaces by Nickel and Puerto (1999) and Puerto and Fern´ andez (2000), respectively. Later, they were extended to the discrete setting by Nickel (2001)andBolandetal.(2006). Given a set of clients and a set of candidate locations and assuming that the allocation costs of clients to facilities are known, the discrete ordered median problem (DOMP) consists of choosing pfacility locations and assigning each client to a chosen facility with the smallest allocation cost to minimize the ordered weighted average of these costs. The ordered weighted average sorts the allocation costs in a nondecreasing sequence and then it performs the scalar product of this soobtained sorted cost vector with a given vector of weights. DOMP has been widely studied since the 1990s, and there are a number of different formulations, solution approaches, and applications available in the literature. To cite a few, DOMP has been applied to discrete facility location in Boland et al. (2006), Mar´ ınetal.(2009, 2010), Nickel (2001), Puerto (2008), and Puerto et al. (2009); to location on networks in Nickel and Puerto (1999); to hub network design problems in Puerto et al. (2011,2013,2016); to determine values in cooperative game theory in Perea and Puerto (2013); to combinatorial optimization problems with ordering in Fern´ andez et al. (2013,2014,2017); and to voting problems in Ponce et al. (2018), etc. The reader is referred to the monographs by Nickel and Puerto (2005) and Puerto and Rodr´ ıguez-Ch´ ıa(2015)forsomeother applications. There exist several valid formulations for DOMP that exploit specific features of the problem, for instance, free self-service, ties in the matrix of costs, or null elements in the vector of weights (see, e.g., Boland et al. (2006), Mar´ ınetal.(2009,2010), Puerto et al. (2013), Labbé et al. (2017), and the references therein). In Labbé et al. (2017), a new formulation for DOMP has been proposed based on a set-packing approach that is valid for general cost coefficients. This formulation gives rise to rather tight integrality gaps and was shown to be reasonably efficient to solve medium-sized instances when embedded in a branch-and-cut (B&C) scheme. For general cost coefficients (with no ties), all these formulations have a very large (cubic) number of 1 binary variables, and therefore, already for instances with 100 clients, they fail to be loaded by the solvers. In this paper, we explore a different paradigm for solving DOMP based on column generation that avoids considering explicitly all the variables of the problem. Moreover, we introduce a new extended formulation using an exponential number of variables that corresponds to a set-partitioning model and provides even better linear relaxation lower bounds. Each variable represents a set of couples (client, position). In each element of the partition, its clients are served by the same facility, and their positions indicate the place of their allocation costs in the sorted list of allocation costsfortheentireconsidered solution. To handle the exponential number of variables, we use a columngeneration approach that is embedded in a branchprice-and-cut (B&P&C) algorithm. Branch-and-price algorithms to solve location problems have been proposed by du Merle and Vial (2002), Lorena and Senne (2004), Senne et al. (2005), Ceselli and Righini (2005), Avella et al. (2006), and Contreras et al. (2011) to cite a few. In most cases, the considered problem is the p-median although there are some exceptions,forinstance,Contrerasetal.(2011)for capacitated hub location or Doulabi et al. (2016)fora different problem. A branch-and-price approach has never been applied to DOMP, and even more, our approach based on a set partition for couples is fully new. These two facts open new avenues of research in the field of location analysis. Therefore, the contributions of this paper are twofold: (1) methodological, to propose a new perspective in the resolution of DOMP based on formulations with an exponential number of variables and to develop an efficient B&P&C algorithm to handle them, and (2) numerical, to provide solutions for large instances of DOMP. Moreover, in those cases in which optimality of a solution cannot be certified, our approach provides, at least, valid lower bounds that can be used to measure the quality of feasible solutions of DOMP given either by heuristic algorithms (Dom´ ınguez-Mar´ ınetal.2005, Stanimirovic et al. 2007, Puerto et al. 2014, Olender and Ogryczak 2018). The remainder of this paper is organized as follows. Notation, models, and algorithms are presented in Section 2.Section2.2 introduces a new set-partitioning formulation for DOMP. This formulation uses an exponential number of variables in which each element of the partition is a set of clients that are assigned to the same facility together with their sorted positions. This formulation is theoretically compared with another valid formulation described in Section 2.1 and borrowed from Labbé et al. (2017). Section 2.3 describes the column generation algorithm that we have designed to overcome the large number of variables in themodel.Weprovethatthepricingsubproblemis solvable efficiently in polynomial time by using an ad hoc dynamic programming algorithm. We devote Section 3to the implementation details of our B&P&C algorithm. We develop a GRASP heuristic in Section 3.1 that is used to generate both a promising initial solution and a pool of variables to initialize the columngeneration routine. We also develop a stabilization routine based on Pessoa et al. (2010)thatreduces considerably the number of iterations of the columngeneration approach in Section 3.2.Inaddition,Section 3.3 is devoted to an additional improvement, namely a preprocessing. The next two sections, 3.4 and 3.5, present our branching strategies and some families of valid inequalities that are added to the branch-and-price algorithm. In the last section, namely Section 4,wereportonthefinal computational experiments. We evaluate the performance of the B&P&C algorithm and compare it to the compact formulation in Section 2.1. The paper ends with some concluding remarks. 2. Problem Definition and Formulations Let Ibe a set of npoints that, at the same time, represent clients and potential uncapacitated facility locations, and let cij denote the cost for serving client i’s demand from facility j. Given a set Jof popen facilities, let ci(J)represent thecheapestcostforallocatingclientito a facility in J,thatis,ci(J):minj∈Jcij. Now let us sort the costs ci(J),i∈Iby nondecreasing order of their values. The elements of the resulting vector of ordered costs are denoted by c(k)(J)and satisfy c(1)(J)≤···≤c(n)(J).We denote the set of all possible positions 1,...,nin this ordered vector by K. Given the vector λ(λk)k∈Ksatisfying λk≥0,k∈K, the objective function of DOMP is defined as z(J):∑ k∈K λkc(k)(J).(1) Recall that this objective function provides a very general paradigm to encompass standard and new location models. For instance, if λ1...λn1, we obtain the median objective; if λ1λ2... λn−10,λ n1, we obtain the center objective; if λ1 λ2...λn−1α, λn1, where α[0,1],weobtain a convex combination of median and center objectives (centdian), etc. The p-facility discrete ordered median problem looks for the subset Jof pfacilities to open to minimize the ordered median function: min J⊆I:|J|pz(J).(DOMP) Several formulations of DOMP have been proposed in the literature using different types of variables. Among Deleplanque et al.: A Branch-Price-and-Cut Procedure for the DOMP 2INFORMS Journal on Computing, Articles in Advance, pp. 1–18, © 2020 INFORMS them, we mention those based on a combination of the p-median and permutation polytopes (Boland et al. 2006) or on covering approaches by using radius variables (Puerto 2008;Mar ´ ınetal.2009,2010). 2.1. An Explicit Formulation for DOMP: The Weak Order Constraints In the following, we recall the weak order constraints formulation, WOC,introducedinLabbéetal.(2017), and that is the starting point for the developments presented in this paper. This formulation uses two types of binary variables. Variable yjassumes value one if facility j∈Iis open (i.e., j∈J) and zero otherwise. Variable xk ij is equal to one if client i∈Iis allocated to facility j∈Iand the corresponding cost occupies position k∈Kin the allocation cost ranking (i.e., c(k)(J) cij). The choice of this formulation is motivated by its good performance in terms of integrality gap (see Labbé et al. (2017)). However, it requests important memory space because it needs O(n3)binary variables, which may become prohibitive for moderate n. We denote the rank of the allocation cost cij by rij, that is, rij if cij is the th element in the list of the costs cij for all i,j∈I, sorted by order of nondecreasing values and for which ties are broken arbitrarily. For the sake of readability, the reader is referred to Example 1in Section 2.3.Theformulationis as follows: (WOC)min ∑ i∈I∑ j∈I∑ k∈K λkcijxk ij (2) s.t. ∑ j∈I∑ k∈K xk ij 1i∈I(3) ∑ i∈I∑ j∈I xk ij 1k∈K(4) ∑ k∈K xk ij ≤yji,j∈I(5) ∑ j∈I yjp(6) ∑ i∈I∑ j∈I∑ i∈I∑ j∈I: rij≤rij xk ij+∑ i∈I∑ j∈I: rij≥rij xk−1 ij ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠≤n2 k∈K,k 1(7) xk ij,yj∈{0,1}i,j∈I,k∈K.(8) By means of (3), we ensure that each location is served by exactly one facility. In the same way, in each position there must be exactly one allocation cost (4). Constraints (5) translate the fact that a client can be allocated to a facility only if this facility is open and that the allocation cost of a client to a facility can be placed in at most one position. The equality constraint (6) implies that there are exactly popen facilities. Constraints (7), called weak order constraints,ensure that, if client iis allocated to facility jand the corresponding costs cij occupy the kth position in the cost ranking of the solution, then, in the (k−1)th position, there must be a smaller allocation cost. This property is enforced by the coefficients of each variable in the inequality. In each constraint, there are two different positions, kand k−1sothat,by(4), only two variables must take value one, and all the others are equal to zero. If we do not take into account the variables assuming the value zero and we assume that the variables with value one for positions kand k−1 correspond to allocation pairs in sorted position sand t,respectively, the inequality reduces to the following expression: (n2−(s−1))xk isjs+txk−1 itjt≤n2,which is valid if and only if t<s. Finally, the variables are binary, see (8). WOC can be reinforced by adding the following valid inequalities: ∑ i∈I∑ j∈I: rij≤rij xk ij+∑ i∈I∑ j∈I: rij≥rij xk−1 ij≤1,i,j∈I,k∈K,k 1.(9) Observe that constraints (7) are the aggregation over i,j∈Iof inequalities (9). These inequalities are the socalled strong order constraints;seeLabbéetal.(2017) for a detailed explanation. 2.2. A Set-Partitioning Formulation From a linear programming relaxation point of view, the preceding formulation is not the strongest one, but it provides a good compromise between the number of required constraints and the quality of its linear relaxation bound; see Labbé et al. (2017). Further, it allows solving to optimality problems of moderate size. Oneofitsdrawbacksistheuseofacubicnumberof variables, which can be prohibitive for large n.A second important problem of most known formulations for DOMP is their high degree of symmetry in case of allocation costs cij or weights (λk)withmany ties. These reasons motivate the introduction of a new formulation based on a different rationale. We observe that a solution for DOMP is a partition of the clients together with their positions in the sorted vector ofcostssothateachsubsetofclientsinthepartitionis allocated to the same facility. Let us consider sets of couples (i,k),wherethefirst component refers to a client iand the second to a position k,namelyS{(i,k):for some i∈I,k∈K}.Further, we denote by 3(I×K)the family of all sets Sfor which all first (respectively, second) coordinates of its couples are different. Associated with each set Sand facility j,wedefine a variable yj Sequal to one if the set Sis part of a feasible solution ((i,k)∈Siff xk ij 1) and zero otherwise. Deleplanque et al.: A Branch-Price-and-Cut Procedure for the DOMP INFORMS Journal on Computing, Articles in Advance, pp. 1–18, © 2020 INFORMS 3 Let Sbe the set of couples whose first coordinate corresponds to the clients allocated to a given facility j of a feasible solution. The positions of these clients in the solution, that is, the second coordinates of couples in Smust be compatible with the ranking of all the allocation costs involved in the solution. Hence, they must, in particular, be compatible with the ranking of the costs cij of the clients iallocated to j. This implies that, for facility j∈J, we only need to consider subsets of couples Sbelonging to 6(j){S∈3(I×K):cij ≤cij for all (i,k),(i,k)∈S,andk<k}. Because, in any feasible solution, each client imust be allocated to a unique facility jand its allocation cost must occupy a unique position kin the sorted list, the following relationship holds: xk ij ∑ S∈6(j):(i,k)∈S yj S,i,j∈I,k∈K.(10) Next, we can evaluate the cost cj Sinduced by the set S provided that its clients are assigned to facility jin a feasible solution: cj S∑ (i,k)∈S λkcij.(11) To simplify the presentation in the following we denote by (i,·) any couple whose first entry is iregardless of the value of the second entry. Analogously, (·,k) denotes any couple whose second entry is kregardless of the value of the first entry. The following valid formulation uses variables yj S and constitutes our master problem (MP): (MP)min ∑ j∈I∑ S∈6(j) cj Syj S(12) s.t.∑ j∈I∑ S∈6(j): (i,·)∈S yj S1i∈I(13) ∑ j∈I∑ S∈6(j): (·,k)∈S yj S1k∈K(14) ∑ S∈6(j) yj S≤1j∈I(15) ∑ j∈I∑ S∈6(j) yj S≤p(16) ∑ n i1∑ n j1∑ S∈6(j): (i,k)∈S rij≤rij yj S+∑ S∈6(j): (i,k−1)∈S rij≥rij yj S ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ≤n2k∈K,k 1(17) yj S∈{0,1}S∈6(j),j∈I. (18) The objective function (12) accounts for the sorted weighted cost of any feasible solution. Constraints (13) ensure that each client appears in exactly one set S. Constraints (14) ensure that each position is taken by exactly one client appearing in one set S.Constraints (15) guarantees that each facility jserves at most one set Sof clients. Inequality (16) states that at most pfacilities will be opened. By the following family of inequalities (17), we enforce the correct sorting of the costs in any feasible solution. Finally, the variables are binary. One can relate MP and WOC. First, remark that, for a given facility j, there is at most one cost cij that occupies a given position k. Hence, the following constraints are valid for WOC: ∑ i∈I xk ij ≤yj,j∈I,k∈K.(19) Let WOC+denote the formulation given by (2)–(8) and (19) and consider the Dantzig–Wolfe reformulation of WOC+in which constraints (5), (19), and (8) constitute the subproblem. The subproblem can be decomposed by facility. On the one hand, the feasible points of the subproblem of a facility jcorrespond one to one to the sets S∈3(I×K).Hence,thisDantzig–Wolfe reformulation of WOC+is given by the master problem in which we consider variables yj Sfor all S∈3(I×K) (instead of only S∈6(j)). More precisely, the variables of WOC+are related to the variables yj Sthrough the following two equations: xk ij ∑ S∈3(I×K):(i,k)∈S yj Si,j∈I,k∈K and yj∑ S∈3(I×K) yj Sj∈I. Moreover, constraints (13) correspond to constraints (3), constraints (14)to(4), constraints (16)to(6), and constraints (17)to(7). Finally, constraints (15)constitute the “convexity”constraints for the subproblems. On the other hand, it is easy to see that the polyhedron of each subproblem, defined by constraints (5) and (19) together with xk ij ≥0andyj≤1, is integer. This implies that the linear relaxations of WOC+and MP in which all sets S∈3(I×K)are considered provide the same bound. By restricting the subsets Sto be considered for each facility jto belong to 6(j), our formulation MP provides, thus, a stronger model. The computational experiments presented in Section 4.3 show that there exist instances for which the linear relaxation of MP provides a strictly better (higher) lower bound than the linear relaxation of WOC. Deleplanque et al.: A Branch-Price-and-Cut Procedure for the DOMP 4INFORMS Journal on Computing, Articles in Advance, pp. 1–18, © 2020 INFORMS Formulation MP can be strengthened by adding valid inequalities borrowed from WOC. Indeed, one can translate valid inequalities (9)intermsoftheyj Svariables so that they can be used in the set partition formulation of DOMP. The translation of (9)resultsin ∑ S∈6(j): (i,k)∈S rij≤rij yj S+∑ S∈6(j): (i,k−1)∈S rij≥rij yj S≤1,i,j∈I,k∈K,k 1. (20) 2.3. Column Generation to Solve the Linear Relaxation of MP (LRMP) Because the number of variables in MP is too large to be handled directly, in this section, we describe a column generation approach to solve it. Let (α, β, γ, δ, ) be the dual variables associated, respectively, to constraints (13)–(17). The dual problem DP of LRMP is (DP)max ∑ i∈I αi+∑ k∈K βk−∑ j∈I γj−pδ−∑ k∈K: k1 n2k (21) s.t.∑ i∈I: (i,·)∈S αi+∑ k∈K: (·,k)∈S βk−γj−δ −∑ i∈I∑ j∈I∑ (i,k)∈S: rij≥rij k1 k+∑ (i,k)∈S: rij≤rij kn k+1 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ≤cj S j∈I,S∈6(j) (22) δ, γj, k≥0j∈I,k∈K,k 1.(23) To apply the column generation procedure, let us assume that we are givena set of columns that definea restricted master problem and denote its linear relaxation by ReLRMP. This problem is solved to optimality, and (α∗,β ∗,γ ∗,δ ∗, ∗) represents its optimal dual solution. See Example 1.Thereducedcost,cj S,of column yj S,namelycj Scj S−zj S,isgivenby cj Scj S+γ∗ j+δ∗+∑ i∈I∑ j∈I∑ (i,k)∈S: rij≥rij k1 ∗ k+∑ (i,k)∈S: rij≤rij kn ∗ k+1 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ −∑ i∈I: (i,·)∈S α∗ i−∑ k∈K: (·,k)∈S β∗ k.(24) If cj S≥0forallj,S∈6(j),thecurrentsolutionofReLRMP is also optimal for the LRMP, and the column generation procedure stops. Otherwise, one has identified one (some) new column(s) to be added to the current reduced master problem to proceed further. In each iteration, ReLRMP and its reduced costs provide lower and upper bounds for the LRMP. Indeed it holds that (Desrosiers and Lübecke 2005) zReLRMP +p·min j∈I,S∈6(j)cj S≤zLRMP ≤zReLRMP,(25) zReLRMP +∑ j∈I min S∈6(j)cj S≤zLRMP ≤zReLRMP,(26) where zReLRMP and zLRMP denote the optimal value of ReLRMP and LRMP, respectively. Example 1. Consider the following vector λ(4,2,1), cost matrix C, and precedence matrix R: C 136 318 681 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠,R 146 528 793 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠. For n3, there are 33 different sets of couples (i,k) in 6. S1{(1,1)} S2{(1,2)} S3{(1,3)} S4{(2,1)} S5{(2,2)} S6{(2,3)} S7{(3,1)} S8{(3,2)} S9{(3,3)} S10 {(1,1),(2,2)} S11 {(1,1),(2,3)} S12 {(1,1),(3,2)} S13 {(1,1),(3,3)} S14 {(1,2),(2,1)} S15 {(1,2),(2,3)} S16 {(1,2),(3,1)} S17 {(1,2),(3,3)} S18 {(1,3),(2,1)} S19 {(1,3),(2,2)} S20 {(1,3),(3,1)} S21 {(1,3),(3,2)} S22 {(2,1),(3,2)} S23 {(2,1),(3,3)} S24 {(2,2),(3,1)} S25 {(2,2),(3,3)} S26 {(2,3),(3,1)} S27 {(2,3),(3,2)} S28 {(1,1),(2,2),(3,3)} S29 {(1,1),(2,3),(3,2)} S30 {(1,2),(2,1),(3,3)} S31 {(1,2),(2,3),(3,1)} S32 {(1,3),(2,1),(3,2)} S33 {(1,3),(2,1),(3,2)}. The sets 6(j)are the following: 6(1){S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S12, S13,S15,S17,S22,S23,S25,S28}, 6(2){S1,S2,S3,S4,S5,S6,S7,S8,S9,S12,S13,S14, S17,S18,S19,S22,S23,S25,S30}, 6(3){S1,S2,S3,S4,S5,S6,S7,S8,S9,S10,S11,S15, S16,S20,S21,S24,S26,S27,S31}. Deleplanque et al.: A Branch-Price-and-Cut Procedure for the DOMP INFORMS Journal on Computing, Articles in Advance, pp. 1–18, © 2020 INFORMS 5 We consider as initial pool of columns the variables y1 18 and y3 8. With this set of variables, the ReLRMP is ( ReLRMP ) min+2y2 5+10y1 13 s.t+y1 13 ≥1i1 +y2 5≥1i2 +y1 13 ≥1i3 +y1 13 ≥1k1 +y2 5≥1k2 +y1 13 ≥1k3 −y1 13 ≥−1j1 −y2 5≥−1j2 ≥−1j3 −y2 5−y1 13 ≥−2 −8y2 5−y1 13 ≥−9k2 −2y2 5−3y1 13 ≥−9k3 y≥0. Actually, we are interested in its dual problem: ( DP ) max +α1+α2+α3+β1+β2+β3−γ1−γ2−γ3 −2δ−92−93 s.t.+α2+β2−γ2−δ−82−23≤2y2 5 () +α1+α3+β1+β3−γ1−δ−2−33≤10 y1 13 () α, β, γ, δ,  ≥0. Solving (DP), the solution is α22,β 310,α 1α3 β1β2δ230 and the value of the objective function is f12. 2.4. Solving the Pricing Subproblem Although any column yj Swith negative reduced cost may be added to ReLRMP, we follow a strategy that identifies the most negative reduced cost for each facility j. This approach may give rise to several candidate columns (multiple pricing; see Chv´ atal 1983), which is advantageous for this procedure. To do that, for each facility j∈I,wesolveasubproblem to find the column yj S,S∈6(j), with minimum reduced cost. This set Smust be such that there is at most one couple (i,·) for each client iand one couple (·,k)for each position k.Furthermore,thesetS must enjoy that the allocation costs if its couples are compatible. We solve this problem by the following dynamic programming algorithm. Let dk ij be the contribution of the pair (i,k)to the reduced cost of any column yj Ssuch that (i,k)∈S. Depending on the values of k,dk ij is given by dk ij  λkcij +∑ i∈I∑ j∈I: rij≤rij k+1−αi−βkif k1, λkcij +∑ i∈I∑ n j∈I: rij≥rij k+∑ i∈I∑ j∈I: rij≤rij k+1−αi−βk if k2,...,n−1, λkcij +∑ i∈I∑ j∈I: rij≥rij k−αi−βk,if kn. ⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩(27) Then, for a facility j,theproblemoffinding the variable yk Swith minimum reduced costs can be formulated as min S∈6(j)cj Sγ∗ j+δ∗+∑ (i,k)∈S dk ij.(28) Now, for each facility j,wedefine a matrix Djas follows: Dj d1 i1jd2 i1j··· dn i1j d1 i2j . . ... . d1 injdn inj ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ,(29) where i1,i2,...,inis a permutation of the indices i 1,...,nsuch that ci1j≤ci2j≤···≤cinj. Example 2 (Continuing from Example 1).We illustrate the procedure that computes the elements dk ij for all i,k1,...,nof the matrix D1(j= 1). d1 11 λ1c11 +r112−α1−β14 d2 11 λ2c11 +n2−r11 +1 () 2+r113−α1−β22 d3 11 λ3c11 ++ n2−r11 +1 () 3−α1−β3−9 d1 21 λ1c21 +r212−α2−β110 d2 21 λ2c21 +n2−r21 +1 () 2+r213−α2−β24 d3 21 λ3c21 ++ n2−r21 +1 () 3−α2−β3−9 d1 31 λ1c31 +r312−α3−β124 d2 31 λ2c31 +n2−r31 +1 () 2+r213−α3−β212 d3 31 λ3c31 ++ n2−r31 +1 () 3−α3−β3−4 Deleplanque et al.: A Branch-Price-and-Cut Procedure for the DOMP 6INFORMS Journal on Computing, Articles in Advance, pp. 1–18, © 2020 INFORMS Because r11 <r21 <r31, the valid permutation is (1,2,3). This implies that D1 42−9 10 4 −9 24 12 −4 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ i1 i2 i3. By using Dj,weobtainthatasetSbelongs to 6(j)if and only if every (i1,k1)and (i2,k2)∈Ssuch that i1<i2: k1<k2. Our dynamic programming algorithm to obtain the minimum reduced cost for each j∈Jbuilds upon this observation by constructing a solution to a reduced version of (28)inwhichonlythefirst ilclients and the first kpositions are considered. For each couple (il,k),wedefine a function gj(il,k)min {cj S:S∈6j ()and for all i l,k () ∈S:i l≤iland k≤k}(30) and we denote an optimal solution of this restricted optimization problem by Sj(il,k). Hence, the optimal value of problem (28)isequalto gj(in,n)+δ+γjand a corresponding optimal solution by Sj(in,n). Our recursive procedure computes gj(il,k)and Sj(il,k)for increasing values of land kand exploits the following feasibility conditions on S: i. For each client i(position k), at most one couple containing i(position k)belongstoS. ii. If (il1,k1)and (il2,k2)∈Sand k1<k2,thenril1j<ril2j. More precisely, if (il,k)belongs to Sj(il,k),then,from (i), it follows that gj(il,k)gj(il−1,k−1)+dk ilj.Otherwise, Sj(il,k)may contain a couple (il,k)withk≤k−1 or a couple (il,k)with l≤l−1 but not both; otherwise, condition (ii) would be violated. Hence, in this case, gj(il,k)min{gj(il−1,k−1),gj(il,k−1),gj(il−1,k)}. Combining the two cases, we obtain the following recurrence relation for l,k2,...,n: gj(il,k)min{gj(il−1,k−1)+dk ilj,gj(il−1,k−1), gj(il,k−1),gj(il−1,k)}.(31) Algorithm 1 (Pricing Subproblem Algorithm) 1: gj(i1,1)min{0,d1 i1j}; 2: if gj(i1,1)d1 i1j<0, then 3: Sj(i1,1){(i1,1)}; 4: else 5: Sj(i1,1)∅; 6: end if 7: for k2,...,n,do 8: gj(i1,k)min{dk i1j,gj(i1,k−1)}; 9: if gj(i1,k)gj(i1,k−1),then 10: Sj(i1,k)Sj(i1,k−1); 11: else 12: Sj(i1,k){(i1,k)}; 13: end if 14: end for 15: for l2,...,n,do 16: gj(il,1)min{d1 ilj,gj(il−1,1)}; 17: if gj(il,1)gj(il−1,1),then 18: Sj(il,1)Sj(il−1,k); 19: else 20: Sj(il,1){(il,1)}; 21: end if 22: end for 23: for k,l2,...,n,do 24: gj(il,k)min{gj(il−1,k−1)+dk ilj,gj(il−1,k−1), gj(il,k−1),gj(il−1,k)}; 25: if gj(il,k)gj(il−1,k−1),then 26: Sj(il,k)Sj(il−1,k−1); 27: else if gj(il,k)gj(il,k−1),then 28: Sj(il,k)Sj(il,k−1); 29: else if gj(il,k)gj(il−1,k),then 30: Sj(il,k)Sj(il−1,k); 31: else 32: Sj(il,k)Sj(il−1,k−1)∪{(il,k)}; 33: end if 34: end for Obviously, if, at the end of the procedure, gj(in,n)+ δ+γjis negative, the variable yj Sj(in,n)is a good candidate to be chosen in the next iteration of the columngeneration scheme. If we solve this problem for all j,wegetcj RminScj S, and if cj R<0, we can activate (at least) yj R.Next,we solve a new ReLRMP with this (these) new activated variable(s). Remark 1. Computing each matrix Djcan be done in O(n2). Next, obtaining gj(in,n)requires the evaluation of the function gj(i,k)for all i∈Iand k∈K. According to the algorithm, the evaluation of each gj(i,k)is done in constant time. Solving the pricing subproblem amounts to evaluating gj(in,n)for all j∈I. Therefore, the entire pricing subproblem can be solved in O(n3)time. Example 3 (Continuing from Example 2).We show the computation of the gj(in,n)and Sj(in,n)for j1. g1(i1,1)min{0,4}0,S1(i1,1)∅ g1(i1,2)min{2,0}0,S1(i1,2)∅ g1(i1,3)min{−9,0}−9,S1(i1,3){(1,3)} g1(i2,1)min{10,0}0,S1(i2,1)∅ g1(i3,1)min{24,0}0,S1(i3,1)∅ g1(i2,2)min{0+4,0,0,0},S1(i2,2)∅ g1(i3,2)min{0+12,0,0,0},S1(i3,2)∅ g1(i2,3)min{0−9,0,−9,0},S1(i2,3){(1,3)} g1(i3,3)min{0−4,0,−9,0},S1(i3,3){(1,3)} Deleplanque et al.: A Branch-Price-and-Cut Procedure for the DOMP INFORMS Journal on Computing, Articles in Advance, pp. 1–18, © 2020 INFORMS 7 We have obtained g1(i3,3)−9, and S1(i3,3)S3is the potential set to be used because its reduced cost is negative. The corresponding reduced cost c1 3g1(i3,3)+ δ+γ1−9+0+0−9<0. Hence, we active variable y1 3. Next, the process continues with the following facilities, that is, j2,3. In this example, the optimal solution can be certifiedafterfourcompleteiterations of the preceding process. 2.5. Dealing with Infeasibility One important issue when implementing a columngeneration procedure to solve a linear optimization problem is how to deal with infeasibility. This is especially crucial if the procedure is used within a branch-and-bound scheme to solve the linear relaxation of the problem at every node of the branching tree. To handle this, we resort to the so-called Farkas pricing. This method was used previously, to the best of our knowledge, in Günlük et al. (2005) and Ceselli et al. (2008). The term “Farkas pricing”was coined in Gamrath (2010). According to Farkas’lemma (Farkas 1894), a reduced master problem is infeasible if its associated dual problem is unbounded. Thus, to recover feasibility in the ReLRMP, we have to revoke the certificate of unboundedness in the dual problem. This can be done by adding constraints to it. Because we are only interested in recovering feasibility in ReLRMP, one can proceed in the same way as for the usual pricing but with null coefficients in the objective function of the primal. In this way, the Farkas dual problem is max ∑ i∈I αi+∑ k∈K βk−∑ j∈I γj−pδ−∑ k∈K: k1 n2k(32) s.t.∑ i∈I: (i,·)∈S αi+∑ n k∈K: (·,k)∈S βk−γj−δ −∑ i∈I∑ j∈I∑ (i,k)∈S: rij≥rij k1 k+∑ (i,k)∈S: rij≤rij kn k+1 ⎛ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎞ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ≤0j∈I,S∈6(j) (33) δ, γjk≥0j∈I,k∈K, k 1.(34) To identify new variables that make the reduced master problem feasible, we use our dynamic programming approach in which we replace the column costs cj Sby zeros. Farkas pricing is an important element in our approach because it allows starting the column generation algorithm with an empty pool of columns although this is not advisable. Furthermore, Farkas pricing is crucial in the branching phase to recover feasibility (whenever possible) in those nodes of the branching tree where it is lost after fixing variables. 3. A Branch-Price-and-Cut Implementation In this section, we precise several components of the implementation of our set-partitioning formulation based on a column generation approach. B&P&C is a branch-and-cut scheme that solves the linear relaxation at each node of the branching tree with the column generation algorithm previously described and may apply cuts to improve the obtained lower bound. (The reader is referred to Doulabi et al. (2016) for another recent implementation of a B&P&C.) Unless otherwise specified, to calibrate the best choice of the different parameters used in our B&P&C, we have performed a preliminary computational study based on a set of 60 instances with sizes n20,30 and with a time limit of 1,800 seconds. Those are the smallest instances that we eventually use in Section 4. 3.1. Upper Bound for the Master Problem: A GRASP Heuristic and an Initialization Stage A heuristic algorithm that generates a good feasible solution for MP provides a promising pool of initial columnsaswellasagoodupperbound. GRASP (Feo and Resende 1989,1995) is a wellknown heuristic technique that usually exhibits good performance in short computing time. In our case, it consists of a multistart greedy algorithm to construct a set of pfacilities from a randomly generated set of facilities with smaller cardinality. Following Puerto et al. (2014), we have chosen, in a greedy manner, an initial set of p/2facilities. Next, we improve this initial solution by performing a fixed number of iterations of a local search procedure. The greedy algorithm adds iteratively a new facility to the current set of open facilities, choosing the one with the maximum improvement of the objective value. The local search consists of an interchange heuristic between open and closed facilities. The pseudocode of the GRASP used to solve the problem is described in Algorithm 2. Algorithm 2 (GRASP for DOMP) 1: Input(n,p,C,λ,n1,n2,q); 2: for n1replications, do 3: PartialSolution ←ConstructRandomizedPartialSolution(q); 4: Solution ←ConstructGreedySolution(PartialSolution); 5: for n2iterations, do 6: Solution ←LocalSearch(Solution); 7: BestSolution ←UpdateSolution(Solution, BestSolution); 8: end for 9: end for Deleplanque et al.: A Branch-Price-and-Cut Procedure for the DOMP 8INFORMS Journal on Computing, Articles in Advance, pp. 1–18, © 2020 INFORMS First of all, we would like to point out the remarkable behavior of the GRASP heuristic for this problem. To illustrate the appropriateness of our heuristic, we have solved to optimality a number of instances of the problem using the mixed integer programming (MIP) formulation to be compared with those given by our GRASP. In all instances, up to a size of n400, the solution provided by GRASP is always as good as the one obtained by the any of our MIP formulations with a CPU time limit of 7,200 seconds; see Section 4. Moreover, it is not only advisable to use the GRASP heuristic because it provides a very good upper bound, thus, helping the exploration of the searching tree by pruning many branches of the branch-and-bound tree, but in addition, the construction phase of the heuristic also provides a very promising pool of initial columns for the B&P&C, in combination with the technique described in the following. Because we are solving the LRMP without generatingitsentiresetofvariables,usingtheprimalsimplex algorithm, the goal of the initialization phase is to find an initial set of columns that allows solving the MP by performing a small number of iterations in the column generation routine. We create variables using amodification of the local search routine of the GRASP algorithm. Every time that we find a promising feasible solution in the heuristic, we create the variables that define that solution (CreateSetVariables(J)). Algorithm 3presents the pseudocode of this process. Function CreateSetVariables(J) determines the costs involved in the solution, that is, the minimum for each client among the open facilities. Then those costs are ordered to determine the position of each client. Once we know the couples (i,k)assigned for each open facility, the corresponding variables are added to the pool. Example 4 (Continuing from Example 1).We illustrate the use of the function CreateSetVariables(J) with the following set J{1,3}(open facilities). The allocation costs for this set Jof open facilities are c11 1,c21 3, c33 1. According to R, the ranks of these costs are r11 1<r33 3<r21 5. Thus, we get the couples (1,1), (3,2)and (2,3). This means that client 1 goes to facility 1 in position 1, client 3 goes to facility 3 in position 2, and client 2 goes to facility 1 in position 3. Therefore, the variables y1 {(1,1),(2,3)} and y3 {(3,2)} are added to the pool. Algorithm 3 (Initial Columns) 1: Input(|J|p); 2: ¯ zz(J); CreateSetVariables(J); 3: for n2iterations, j1∈J,j2∈¯ J,do 4: if z((J\{j1}) ∪ {j2}) <¯ z,then 5: ¯ zz((J\{j1}) ∪ {j2});J(J\{j1}) ∪ {j2}; CreateSetVariables(J); 6: end if 7: end for To test the usefulness of GRASP in solving problem instances, Table 1reports results for the 60 instances of sizes n20,30, enabling or not the use of the GRASP. It shows average results of CPU time (Time(s)), percentage gap at termination, i.e., 100(zUB − zLB)/zLB (Gap(%)), and number of unsolved problems (in parentheses), number of nodes |Nodes|, and number of variables (|Vars|). According to Table 1, it is clearly advisable to use the upper bound provided by the GRASP heuristic: it reduces the number of nodes, thus improving the size of the branch-and-bound tree. In Table 2, the same information as in Table 1is reported but only for the instances solved to optimality within the time limit. One can observe from this table that enabling the use of GRASP reduces the CPU time and number of nodes of the B&B tree and, at the same time, reduces the overall number of variables required by the B&P&C. In addition, by using the GRASP heuristic, B&P&C is able to solve six more instances. For those instances for which B&P&C does not certify optimality, GRASP provides an upper bound that leads to an average gap of 0.89%. Finally, without the use of GRASP, in many cases, no feasible solutions are found within the time limit, and thus, no percentage gap can be reported. Our results show that, by using the GRASP heuristic, 2.03% of the final number of variables are generated when applying Algorithm 3.Thecombination of the incumbent solution (given by GRASP) and that initial pool of variables leads to solving the considered instances faster, requiring a fewer number of nodes and variables to certify optimality. Figure 1reports the performance profile of GAP versus number of solved instances within a time limit of 1,800 seconds for the 60 instances with sizes n 20,30. The dashed line reports results using GRASP and the solid one without it. It is interesting to point out that, when GRASP is enabled, the B&P&C is able to solve to optimality 30 instances, and the GAP Table 2. CPU Time, Number of Nodes, and Number of Variables with and without GRASP Heuristic for n20,30 GRASP Time(s) |Nodes| |Vars| Disabled 386.19 272 10,962 Enabled 147.77 79 6,062 Note. Summary of solved instances. Table 1. CPU Time, Number of Nodes, and Number of Variables with and without GRASP Heuristic for n20,30 GRASP Time(s) Gap(%) |Nodes| |Vars| Disabled 1,107.21 —(36) 158 12,850 Enabled 965.31 0.89 (30) 88 9,907 Deleplanque et al.: A Branch-Price-and-Cut Procedure for the DOMP INFORMS Journal on Computing, Articles in Advance, pp. 1–18, © 2020 INFORMS 9 column #unsolved by Memory(MB). This new column shows the average required memory to solve the corresponding set of instances. In that table, extensive computational experiments are reported for instances up to 400 points. We would like to remark that the increase of the complexity with respect to the instance sizes of |Vars|,|Cuts|,Memory(MB),andGap(%)is moderate (almost linear), which allows one to handle DOMP problemsoflargersize.MoreovertheGap(%)are similar to those reported in Table 5. To conclude, the results show that the overall performance of B&P&C(MP) in solving DOMP is systematically better than the branch-and-cut formulation B&C(WOC) for instances of n≥70. In addition, it is worth noting that B&C(WOC) is not even able to solve the linear relaxation of DOMP problems of sizes n≥100. This fact shows the usefulness of our new approach. 5. Conclusions This paper presents a first branch-price-and-cut, B&P&C(MP), algorithm for solving DOMP. This approach is based on an extended formulation using an exponential number of variables coming from a setpartitioning model. Elements in the partitionsare couples containing information about a client and its sorted position in the sorted sequence of allocation costs. To address the solution of this formulation, we develop a column generation algorithm, and we prove that the pricing routine is polynomially solvable by a dynamic programming algorithm. We embed the column generation algorithm within a branch-and-price framework. Furthermore, we adapt preprocessing and incorporate families of valid inequalities that improve its performance. Extensive computational results compare the performanceofourB&P&C(MP)againstthemostrecent algorithm in the literature for DOMP, B&C(WOC), showing that, for the largest considered instances, B&P&C(MP) performs better, and it requires less memory to upload and run the models. The methodology presented in this paper is able to solve sized instances for DOMP that had never been solved in the literature. Table 6. Numerical Results for B&P&C(MP) for Bigger Instances n p Time(s)|Vars||Nodes||Cuts|Memory(MB)Gap(%) 100 50 7,201.10 23,057 2 5,221 311 2.43 120 60 7,200.54 29,209 1 6,992 424 2.54 140 70 7,202.11 34,222 2 6,805 449 2.87 160 80 7,201.08 41,811 1 8,069 574 3.13 180 90 7,201.24 50,523 1 7,961 656 3.48 200 100 7,201.61 59,931 1 9,050 805 3.96 220 110 7,201.39 57,685 1 11,097 806 4.39 240 120 7,200.94 63,710 1 11,487 874 4.43 260 130 7,201.83 73,143 1 9,772 910 5.05 280 140 7,200.98 83,892 1 9,592 1,037 6.19 300 150 7,201.74 87,444 1 10,811 1,076 6.83 100 33 7,201.10 29,974 2 5,864 562 3.81 120 40 7,201.12 37,204 5 5,433 621 4.80 140 46 7,200.98 48,286 2 7,953 894 4.47 160 53 7,201.55 57,060 3 7,159 926 5.01 180 60 7,203.85 70,077 1 8,997 1,188 5.99 200 66 7,200.55 67,417 1 11,416 1,199 6.16 220 73 7,201.41 71,731 2 9,356 1,073 6.88 240 80 7,202.02 84,466 2 10,093 1,324 9.80 260 86 7,200.74 92,401 1 10,505 1,390 9.54 280 93 7,201.68 105,535 1 13,265 1,728 9.01 300 100 7,200.28 114,427 1 15,595 2,034 9.41 100 25 7,202.69 31,218 1 7,034 801 4.73 120 30 7,200.98 40,074 1 6,818 922 5.91 140 35 7,203.38 49,539 1 8,466 1,162 5.83 160 40 7,201.35 58,464 1 10,688 1,465 7.29 180 45 7,203.43 72,146 1 10,329 1,674 8.70 200 50 7,201.88 62,568 1 10,252 1,364 11.27 220 55 7,203.56 74,359 1 9,985 1,515 10.34 240 60 7,200.19 80,228 1 9,745 1,422 11.57 260 65 7,200.54 99,628 1 7,944 1,346 11.77 280 70 7,203.63 109,544 1 4,187 1,247 11.32 300 75 7,200.79 128,462 1 3,844 1,261 12.26 400 200 86,400.42 158,287 1 9,308 1,183 7.07 400 133 86,401.23 178,265 1 13,652 2,764 11.26 400 100 86,401.11 236,973 1 7,582 3,125 10.29 Note. For instances with n = 400, the time limit was set to 24 hours. Deleplanque et al.: A Branch-Price-and-Cut Procedure for the DOMP 16 INFORMS Journal on Computing, Articles in Advance, pp. 1–18, © 2020 INFORMS Acknowledgments The authors thank the SCIP team (Gamrath et al. 2016) for the helpful advice. They also thank the helpful reports from two anonymous referees that led to improving the quality of the paper. References Achterberg T, Koch T, Martin A (2005) Branching rules revisited. Oper. Res. Lett. 33(1):42–54. Applegate D, Bixby R, Chv´ atal V, Cook W (1995) Finding cuts in the TSP (a preliminary report). DIMACS Technical Report 95-05, DIMACS, Rutgers University, New Brunswick, NJ. Avella P, Sassano A, Vasil’ev I (2006) Computational study of largescale p-median problems. Math. Programming 109(1):89–114. Barnhart C, Johnson E, Nemhauser G, Savelsbergh M, Vance P (1998) Branch-and-price: Column generation for solving huge integer programs. Oper. Res. 46(3):316–329. Beale E, Tomlin J (1970) Special facilities in a general mathematical programming system for non-convex problems using ordered sets of variables. Lawrence J, ed. Proc. 5th Internat. Conf. Oper. Res. (Tavistock Publications, London), 447–454. Benichou M, Gauthier J, Girodet P, Hentges G, Ribiere G, Vincent O (1971) Experiments in mixed-integer programming. Math. Programming 1(1):71–94. Boland N, Dom´ ınguez-Mar´ ın P, Nickel S, Puerto J (2006) Exact procedures for solving the discrete ordered median problem. Comput. Oper. Res. 33(11):3270–3300. Ceselli A, Righini G (2005) A branch-and-price algorithm for the capacitated p-median problem. Networks 45(3):125–142. Ceselli A, Liberatore F, Righini G (2008) A computational evaluation of a general branch-and-price framework for capacitated network location problems. Ann. Oper. Res. 167(1):209–251. Chv´ atal V (1983) Linear Programming (W. H. Freeman and Company, New York). Contreras I, D´ ıaz J, Fern´ andez E (2011) Branch and price for largescale capacitated hub location problems with single assignment. INFORMS J. Comput. 23(1):41–55. Deleplanque S, Labbé M, Ponce D, Puerto J (2018) An extended version of a branch-price-and-cut procedure for the discrete ordered median problem. Preprint, submitted February 9, https:// arxiv.org/abs/1802.03191. Desrosiers J, Lübecke M (2005) A primer in column generation. Desaulniers G, Desrosiers J, Salomon MM, eds. Column Generation (Springer, Boston), 1–32. Dom´ ınguez-Mar´ ın P, Nickel S, Hansen P, Mladenovic N (2005) Heuristic procedures for solving the discrete ordered median problem. Ann. Oper. Res. 136(1):145–173. Doulabi SHH, Rousseau LM, Pesant G (2016) A constraint-programmingbased branch-and-price-and-cut approach operating room planning and scheduling. INFORMS J. Comput. 28(3):432–448. du Merle O, Vial J (2002) Proximal ACCPM, a cutting plane method for column generation and Lagrangean relaxation: Application to the p-median problem. Technical report, HEC Gen` eve, Université de Gen` eve, Geneva, Switzerland. du Merle O, Villenueve D, Desrosiers J, Hansen P (1999) Stabilized column generation. Discrete Math. 194(1–3):229–237. Farkas G (1894) A Fourier-féle mechanikai elv alkalmaz´ asai. Mathematikai és Természettudom´ anyi ´ Erstesit¨ o12:457–472. Feo TA, Resende MGC (1989) A probabilistic heuristic for a computationally difficult set covering problem. Oper. Res. Lett. 8(2):67–71. Feo TA, Resende MGC (1995) Greedy randomized adaptive search procedures. J. Global Optim. 6(2):109–133. Fern´ andez E, Pozo MA, Puerto J (2014) Ordered weighted average combinatorial optimization: Formulations and their properties. Discrete Appl. Math. 169(31):97–118. Fern´ andez E, Puerto J, Rodr´ ıguez-Ch´ ıa AM (2013) On discrete optimization with ordering. Ann. Oper. Res. 207(1):83–96. Fern´ andez E, Pozo MA, Puerto J, Scozzari A (2017) Ordered weighted average optimization in multiobjective spanning tree problems. Eur. J. Oper. Res. 260(31):886–903. Gamrath G (2010) Generic branch-cut-and-price. Master’s thesis, Institut für Mathematik, Technische Universit¨ at Berlin, Berlin. Gamrath G, Fischer T, Gally T, Gleixner AM, Hendel G, Koch T, Maher SJ, Miltenberger M, Müller B, Pfetsch ME, et al. (2016) The SCIP optimization suite 3.2. Technical report 15-60, ZIB, Berlin. Günlük O, Lad´ anyi L, de Vries S (2005) A branch-and-price algorithm and new test problems for spectrum auctions. Management Sci. 51(3):391–406. Johnson EL (1989) Modeling and strong linear programs for mixed integer programming. Wallace S, ed. Algorithms and Model Formulations in Mathematical Programming,NATOASISeries,vol.51 (Springer, Berlin Heidelberg), 1–43. Labbé M, Ponce D, Puerto J (2017) A comparative study of formulations and solution methods for the discrete ordered p-median problem. Comput. Oper. Res. 78:230–242. Lorena L, Senne E (2004) A column generation approach to capacitated p-median problems. Comput. Oper. Res. 31(6):863–876. Mar´ ın A, Nickel S, Velten S (2010) An extended covering model for flexible discrete and equity location problems. Math. Methods Oper. Res. 71(1):125–163. Mar´ ın A, Nickel S, Puerto J, Velten S (2009) A flexible model and efficient solution strategies for discrete location problems. Discrete Appl. Math. 157(5):1128–1145. Nickel S (2001) Discrete ordered Weber problems. Oper. Res. Proc. 2000: 71–76. Nickel S, Puerto J (1999) A unified approach to network location problems. Networks 34(4):283–290. Nickel S, Puerto J (2005) Location Theory: A Unified Approach (Springer, Berlin, Heidelberg). Olender P, Ogryczak W (2018) A revised variable neighborhood search for the discrete ordered median problem. Eur. J. Oper. Res. 274(2):445–465. Perea F, Puerto J (2013) Finding the nucleolus of any n–person cooperative game by a single linear program. Comput. Oper. Res. 40(10):2308–2313. Pessoa A, Uchoa E, Aragão MP, Rodrigues R (2010) Exact algorithm over an arctime-indexed formulation for parallel machine scheduling problems. Math. Programming Comput. 2(3–4): 259–290. Ponce D, Puerto J, Ricca F, Scozzari A (2018) Mathematical programming formulations for the efficient solution of the k-sum approval voting problem. Comput. Oper. Res. 98:127–136. Puerto J (2008) A new formulation of the capacitated discrete ordered median problem with {0,1}-assignment. Kalcsics J, Nickel S, eds. Oper. Res. Proc. 2007 (Springer, Berlin Heidelberg), 165–170. Puerto J, Fern´ andez F (2000) Geometrical properties of the symmetrical single facility location problem. J. Nonlinear Convex Anal. 1(3):321–342. Puerto J, Rodr´ ıguez-Ch´ ıa AM (2015) Ordered median location problems. Laporte G, Nickel S, Saldanha da Gama F, eds. Location Science (Springer) Heidelberg: 249–288. Puerto J, Pérez-Brito D, Garc´ ıa-Gonz´ alez CG (2014) A modified variable neighborhood search for the discrete ordered median problem. Eur. J. Oper. Res. 234(1):61–76. Puerto J, Ramos AB, Rodr´ ıguez-Ch´ ıa AM (2011) Single-allocation ordered median hub location problems. Comput. Oper. Res. 38(2): 559–570. Puerto J, Ramos AB, Rodr´ ıguez-Ch´ ıa AM (2013) A specialized branch & bound & cut for single-allocation ordered median hub location problems. Discrete Appl. Math. 161(16–17):2624–2646. Deleplanque et al.: A Branch-Price-and-Cut Procedure for the DOMP INFORMS Journal on Computing, Articles in Advance, pp. 1–18, © 2020 INFORMS 17 Puerto J, Rodr´ ıguez-Ch´ ıa AM, Tamir A (2009) Minimax regret single-facility ordered median location problems on networks. INFORMS J. Comput. 21(1):77–87. Puerto J, Ramos AB, Rodr´ ıguez-Ch´ ıa AM, S´ anchez-Gil MC (2016) Ordered median hub location problems with capacity constraints. Transportation Res., Part C Emerging Tech. 70:142–156. Ryan DM, Foster A (1981) An integer programming approach to scheduling. Wren A, ed. Computer Scheduling of Public Transport: Urban Passenger Vehicle and Crew Scheduling (North-Holland, Amsterdam), 269–280. Senne E, Lorena L, Pereira MA (2005) A branch-and-price approach to p-median location problems. Comput. Oper. Res. 32(6):1655–1664. Stanimirovic Z, Kratica J, Dugosija D (2007) Genetic algorithms for solving the discrete ordered median problem. Eur. J. Oper. Res. 182(3):983–1001. Wolsey LA (1998) Integer Programming (John Wiley & Sons, New York). Deleplanque et al.: A Branch-Price-and-Cut Procedure for the DOMP 18 INFORMS Journal on Computing, Articles in Advance, pp. 1–18, © 2020 INFORMS