A comparative study of formulations and solution methods for the discrete ordered p-median problem

A Comparative Study of Formulations and Solution Methods for the Discrete Ordered p-Median Problem Martine Labb´ea, Diego Ponceb,c, Justo Puertob,c aComputer Science Department, Universit´e Libre de Bruxelles bIMUS, Instituto de Matem´aticas, Universidad de Sevilla cDepartamento de Estad´ıstica e Investigaci´on Operativa, Universidad de Sevilla Abstract This paper presents several new formulations for the Discrete Ordered Median Problem (DOMP) based on its similarity with some scheduling problems. Some of the new formulations present a considerably smaller number of constraints to define the problem with respect to some previously known formulations. Furthermore, the lower bounds provided by their linear relaxations improve the ones obtained with previous formulations in the literature even when strengthening is not applied. We also present a polyhedral study of the assignment polytope of our tightest formulation showing its proximity to the convex hull of the integer solutions of the problem. Several resolution approaches, among which we mention a branch and cut algorithm, are compared. Extensive computational results on two families of instances, namely randomly generated and taken from Beasley’s OR-library, show the power of our methods for solving DOMP. Keywords: discrete multifacility location, MIP formulations, ordered median problem 1. Introduction Recognizing the need for more flexible logistic models, Nickel [13] proposed the Discrete Ordered Median Location Problem (DOMP) which could be used to model different locations problems, as the p-median or the p-center. It is a flexible formulation based on applying an ordered weighted averaging operator to the costs as they appear in the solution and taking them into account with a suitable n-vector λ. Given a vector of weights, the ordered weighted average of nreal numbers is obtained by first ranking those numbers by nondecreasing order and then computing the scalar product of the ranked allocation cost vector and the weight vector, see e.g. Nickel and Puerto [14]. Consider a set of clients and a set of candidate locations where some facility can be established. Further we are given the costs for allocating clients to facilities. DOMP consists in choosing pfacility locations and assigning each client to a facility with smallest allocation cost in order to minimize a special objective function the so called ordered median function. Given a vector of weights, this function consists in an ordered weighted average of the allocation costs, namely it sorts these costs in non-decreasing sequence and then it performs the scalar product of this so obtained sorted cost vector with the given vector of weights. This objective function has been widely applied in the field of location analysis and distribution models (Mar´ın et al. [10], Kalcsics et al. [9] and Puerto et al. [17]). In addition, it has the potential to yield new models for order statistics embedded within mathematical programming formulations; thus enlarging the applications of optimization tools to the resolution of statistical problems in data analysis. DOMP is known to be NP-hard, see Nickel and Puerto [14]. The first formulation of DOMP, proposed by Nickel [13] consists in an integer nonlinear problem. Then, Boland et al. [3] propose several linearization of Nickel’s model. Email addresses: [email protected] (Martine Labb´e), [email protected] (Diego Ponce), [email protected] (Justo Puerto) Preprint submitted to Elsevier May 18, 2016 Instances with up to 30 clients could be solved to optimality by Boland et al. [3]. Further, if clients and facility locations coincide and if the allocation cost of a client to itself is equal to zero (the so-called free self service), then instances with up to 100 clients could be solved by Mar´ın et al. [10, 11]. We observe that in all previously considered formulations the gaps with respect to the linear programming relaxations of those models are rather large, as mentioned in all those papers. In this paper, we propose new formulations for DOMP and we develop a theoretical comparison of the lower bounds obtained from their LP-relaxations and show that our new formulations are rather tight. Our theoretical results also attempt to shed some light on the polyhedral structure of the new formulation based on scheduling constraints. We conclude with extensive computational experiments to compare the respective efficiency of these formulations. For the ease of presentation, we assume in some cases in this paper that the client and facility location sets coincide. However, it is important to remark that all the models and results presented carry out to the general case in which client and facility location may differ since we do not impose that the cost of allocating a client to itself is equal to zero. The remaining paper is organized as follows: in Section 2 we define the problem and some previous formulations. Next, we present new formulations for the DOMP. In addition, we analyze the relationship between the polytopes of the previously known formulations of this problem and we identify facets for the related assignment polytope in Section 3. Finally, some computational experiments are reported in Section 4. 2. The problem and some formulations Let Ibe a set of npoints which at the same time represent clients and potential facility locations. The cost for serving client i’s demand from facility jis denoted by Cij and a facility can serve as many clients as possible, i.e. facilities are uncapacitated. The Discrete Ordered Median Problem (DOMP ) consists in (i) determining a subset Jof pfacility locations, J⊂I, to open and (ii) assigning clients to closest open facilities in order to minimize the ordered median objective function defined as follows. Given the set Jof popen facilities, let ci(J) represents the cost for allocating client ito some facility in J such that ci(J) = min j∈JCij. Now let us rank the costs ci(J), i∈Iby non-decreasing order of their values. These ordered costs are denoted by ck ≤(J) and verify c1 ≤(J)≤ · · · ≤ cn ≤(J). Then, given a vector λ= (λk)n k=1 satisfying λk≥0, k = 1, . . . , n, the DOMP objective function, also called ordered median function, is defined as n X k=1 λkck ≤(J).(1) Note 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 we obtain a convex combination of median and center objectives (centdian), etcetera. We define the p-facility Discrete Ordered Median Problem as determining the subset J, of pfacilities to open in order to minimize the ordered median function: min J⊂I:|J|=p n X k=1 λkck ≤(J).(DOMP) 2 2.1. Three-index formulation The formulation that we present below, denoted by (DOMP1), was introduced by Boland et al. [3]. It uses three-index variables xk ij such that xk ij = 1, if client iis served by facility jand cost ck(J) = Cij is the k-th smallest in the ordered sequence c≤(J) and xk ij = 0 otherwise. Further, it also uses location variables yjsuch that yj= 1 if j∈Jand yj= 0 otherwise. If xk ij = 1, we say that allocation of client ito facility jis in position k, or that couple ij is in position k. (DOMP1) min n X i=1 n X j=1 n X k=1 λkCijxk ij (2) s.t. n X j=1 n X k=1 xk ij = 1 i= 1, . . . , n (3) n X i=1 n X j=1 xk ij = 1 k= 1, . . . , n (4) n X k=1 xk ij ≤yji, j = 1, . . . , n (5) n X j=1 yj=p(6) n X i=1 n X j=1 Cijxk−1 ij ≤ n X i=1 n X j=1 Cijxk ij k= 2,· · · , n (7) xk ij ∈ {0,1}i, j, k = 1, . . . , n (8) yj∈ {0,1}j= 1, . . . , n (9) 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 (4). We know that a client can be allocated to a facility only if this facility is open, i.e. xk ij ≤yjfor all i, j, k. Furthermore, each allocation of client to facility can be placed in at most one position. Hence, xk ij ≤yjcan be strengthened yielding constraint (5). The equality constraint (6) implies that there are exactly popen facilities. Inequality (7) imposes that the allocation cost in position k−1 cannot be greater than the one in position k. Finally, the variables are binary, see (8) and (9). 2.2. Two-index formulation This formulation(DOMP2) was described for the first time in Puerto [16] and Mar´ın et al. [10] and later applied to a hub problem in Puerto et al. [17]. It considers a vector that contains all the different values in the cost matrix C, augmented with zero if it is not present in matrix C, as it is explained below. Let Cbe a matrix and assume that it contains Gdifferent values such that cij >0. Then, the (G+ 1)- vector c(.)is constructed as follows c(0) = 0 < c(1) < c(2) <· · · < c(G−1) < c(G)= max{Cij :i, j = 1, . . . , n} To formulate the problem we need to define the following binary variables. Variable xij = 1 if client iis served by facility jand 0 otherwise, variable yj= 1 if j∈Jand 0 otherwise and variable ukh = 1 if the k-th smallest allocation cost is greater than c(h−1) and 0 otherwise. Further, we set uk0= 1 and ukG+1 = 0, k = 1, . . . , n. 3 The problem to solve is (DOMP2) min n X k=1 G X h=1 λk(c(h)−c(h−1))ukh (10) s.t. n X j=1 yj=p(11) n X j=1 xij = 1 i= 1, . . . , n (12) xij ≤yji, j = 1, . . . , n (13) ukh ≥ukh+1 k= 1, . . . , n, h = 1, . . . , G −1 (14) uk+1h≥ukh k= 1, . . . , n −1, h = 1, . . . , G (15) n X i=1 n X j=1: Cij >c(h−1) xij = n X k=1 ukh h= 1, . . . , G (16) xij ∈ {0,1}i, j = 1, . . . , n (17) ukh ∈ {0,1}k= 1, . . . , n, h = 1, . . . , G (18) yj∈ {0,1}j= 1, . . . , n. (19) The objective function (10) is equivalent to (1). Assume that the k-th smallest allocation cost is equal to c(hk)for some hkthen by the definition of the variable ukhk G X h=1 (c(h)−c(h−1))ukh = hk X h=1 (c(h)−c(h−1)) = c(hk)−c(0) =c(hk) provided that ukhk= 1 and ukhk+1 = 0. Finally, n X k=1 G X h=1 λk(c(h)−c(h−1))ukh = n X k=1 λk G X h=1 (c(h)−c(h−1))ukh = n X k=1 λkc(hk)= n X i=1 λici(J). The equality constraint (11) ensures that there are exactly pfacilities to be located. Constraints (12) and (13) state that each client is served by one open facility. Equality (16) ensures a good definition of the variable ukh and it relates the sorting (ukh) and design variables (xij). We need to impose some sorting constraints on the ukh variables (15). Constraint (14) is redundant but it is included because, according to Puerto [16], it significantly strengthen the formulation. Furthermore, all variables are binary, (17), (18) and (19). Note that two more two index formulations has been proposed in Mar´ın et al. [10, 11]. However, they are only valid if cii = 0 ∀i= 1, . . . , n (free self-service). 2.3. A new formulation for DOMP There are several formulations for DOMP but they all have large integrality gap, as observed previously in the literature see Boland et al. [3], Mar´ın et al. [10, 11] and Puerto [16]. The main motivation for addressing a new formulation for the DOMP relies on the attempt to reduce this gap. 4 First, we introduce the following notation. ij ≺ˆıˆ ≡            Cij < Cˆıˆ or Cij =Cˆıˆ and i < ˆı or Cij =Cˆıˆ, i = ˆı and j < ˆ ij ˆıˆ ≡            Cij > Cˆıˆ or Cij =Cˆıˆ and i > ˆı or Cij =Cˆıˆ, i = ˆı and j > ˆ ij ˆıˆ ≡   ij ≺ˆıˆ or ij = ˆıˆ ij ˆıˆ ≡   ij ˆıˆ or ij = ˆıˆ (20) 2.3.1. New three index formulation. Our new formulation uses the same variables and constraints as in the three-index formulation, DOMP1. Except that (7) is replaced by (21) that are called order constraints. The resulting formulation is denoted DOMP3. (DOMP3) min n X i=1 n X j=1 n X k=1 λkCijxk ij s.t. (3),(4),(5),(6) n X ˆı=1 n X ˆ=1: ˆıˆij xk−1 ˆıˆ + n X ˆı=1 n X ˆ=1: ˆıˆij xk ˆıˆ ≤1i, j = 1,· · · , n, k = 2,· · · , n (21) xk ij ∈ {0,1}i, j, k = 1, . . . , n yj∈ {0,1}j= 1, . . . , n The rationale behind constraints (21) is illustrated by the following example. Example 1. Consider the following matrix. C=     0274 1055 3602 9410    (22) The order of couples ij by means of the above preference order is 11 ≺22 ≺33 ≺44 ≺21 ≺43 ≺12 ≺34 ≺31 ≺14 ≺42 ≺23 ≺24 ≺32 ≺13 ≺41. The columns of Figure 1 represent the n2possible assignments of clients to facilities whereas its rows represent the npositions in DOMP objective function. Each point (bullet or circle) thus represents a variable xk ij and the bullets correspond to variables which cannot take value 1 simultaneously because, in any feasible solution, the cost of couples assigned to consecutive positions should be non-decreasing. k=1 k=2 k=3 k=4 ij 11 22 33 44 21 43 12 34 31 14 42 23 24 32 13 41 Figure 1: An order constraint for the case n=4 5 So, Figure 1 corresponds to the following inequality: n X i=1 n X j=1: ˆıˆ43 x2 ˆıˆ + n X i=1 n X j=1: ˆıˆ43 x3 ˆıˆ ≤1, or equivalently x3 11 +x3 22 +x3 33 +x3 44 +x3 21 +x3 43 +x2 43 +x2 12 +x2 34 +x2 31 +x2 14 +x2 42 +x2 23 +x2 24 +x2 32 +x2 13 +x2 41 ≤1. Note that the number of order constraints is O(n3). Further they can be seen as cliques of a conflict graph induced by the incopatibility among xk ij variables on three index formulations (Nemhauser and Trotter [12], Fern´andez, Puerto, Rodr´ıguez-Ch´ıa [7]). In a similar way, we can choose several allocations which are sorted and identify a new type of constraints. Let sbe a positive integer with s≤n−1 and let (i1j1),(i2j2),...,(isjs) be scouples of clients and facilities such that irjrir+1jr+1 for all r= 1, . . . , s −1. Then the following family of inequalities, called staircase inequalities is valid for k=s, . . . , n: n X i=1 n X j=1: iji1j1 xk ij + s−1 X r=1 n X i=1 n X j=1: ijirjr ijir+1jr+1 xk−r ij + n X i=1 n X j=1: ijisjs xk−s ij ≤1.(23) Figure 2 provides an example of staircase inequality when n= 5. k=1 k=2 k=3 k=4 i1j1 k=5 i2j2 i3j3i4j4 Figure 2: A staircase constraint for the case n=5 Notice that there exists an exponential number of additional staircase constraints. But the question is whether these new constraints actually strengthen our formulation. The answer is provided by the following result that states that all of them are implied by those with a single step, i.e. the order constraints. Proposition 1. Staircase inequalities (23) with irjrir+1jr+1 are all dominated by (21) and (4). Proof. Let sbe a positive integer such that s≤n−1, then by (21) we obtain that also the following s inequalities hold: n X i=1 n X j=1: ijirjr xk−(r−1) ij + n X i=1 n X j=1: ijirjr xk−r ij ≤1r= 1, . . . , s. (24) Now, since irjrir+1jr+1 using (4), we obtain the following s−1 equations: n X i=1 n X j=1 xk−r ij = 1 r= 1, . . . , s −1.(25) 6 Adding all inequalities (24) we obtain a new inequality: n X i=1 n X j=1: iji1j1 xk ij + n X i=1 n X j=1: iji1j1 xk−1 ij +· · · + n X i=1 n X j=1: ijisjs xk−(s−1) ij + n X i=1 n X j=1: ijisjs xk−s ij ≤s. After rearranging, we get n X i=1 n X j=1: iji1j1 xk ij + s−1 X r=1 n X i=1 n X j=1: ijirjr xk−r ij + s−1 X r=1 n X i=1 n X j=1: ijir+1jr+1 xk−r ij + n X i=1 n X j=1: ijisjs xk−s ij ≤s. Next, we conveniently split some terms of the above inequality to get: n X i=1 n X j=1: iji1j1 xk ij + s−1 X r=1 n X i=1 n X j=1: ijirjr ijir+1jr+1 xk−r ij + s−1 X r=1 n X i=1 n X j=1: ijir+1jr+1 xk−r ij + s−1 X r=1 n X i=1 n X j=1: ijirjr xk−r ij + s−1 X r=1 n X i=1 n X j=1: ijirjr ijir+1jr+1 xk−r ij + n X i=1 n X j=1: ijisjs xk−s ij ≤1+(s−1). (26) On the other hand, using equality (25) we can write n X i=1 n X j=1: ijirjr xk−r ij + n X i=1 n X j=1: ijirjr ijir+1jr+1 xk−r ij + n X i=1 n X j=1: ijir+1jr+1 xk−r ij = 1 r= 1, . . . , s −1. Adding the above equations for all r= 1, . . . , s −1, we obtain s−1 X r=1 n X i=1 n X j=1: ijirjr xk−r ij + s−1 X r=1 n X i=1 n X j=1: ijirjr ijir+1jr+1 xk−r ij + s−1 X r=1 n X i=1 n X j=1: ijir+1jr+1 xk−r ij =s−1. Finally, using the above equation in (26) results in n X i=1 n X j=1: iji1j1 xk ij + s−1 X r=1 n X i=1 n X j=1: ijirjr ijir+1jr+1 xk−r ij + n X i=1 n X j=1: ijisjs xk−s ij ≤1, which is a staircase inequality. 2.3.2. A two index formulation with scheduling constraints. Formulation DOMP3is rather efficient and tight whenever there are few ties in the structure of the allocation costs. This is for instance the case of problems with assignment costs based on flat costs (for instance randomly generated). However, if the number of ties in the allocation costs is large the number of binary variables and constraints is relatively large, as compared with similar numbers in formulation DOMP2. In order to exploit this advantage without losing the usage of scheduling constraints, that relate the formulation with the stable set problem, we will develop in the following another formulation. The new formulation, called DOMP3C, is an extension of DOMP3using the rationale of DOMP2. Moreover, it provides a compact form for representing ties in the allocation costs. 7 Consider a new set of binary variables, vkh such that vkh = 1 if the k-th smallest allocation cost is c(h) and 0 otherwise. Next, DOMP3Cis a new valid formulation for DOMP : (DOMP3C) min n X k=1 G X h=0 λkc(h)vkh (27) s.t. (11),(12),(13) G X h=0 vkh = 1 k= 1, . . . , n (28) n X i=1 n X j=1: Cij =c(h) xij = n X k=1 vkh h= 0, . . . , G (29) X h0<h vk+1h0+X h0≥h vkh0≤1k= 1, . . . , n −1 h= 1, . . . , G (30) xij, yj, vkh ∈ {0,1}i, j, k = 1, . . . , n h= 0, . . . , G. (31) Clearly, the objective function (27) accounts for the ordered weighted sum of the allocation costs. Constraints (29) state that the number of allocations that are attained at the value c(h)regardless of the level kthat they occupy (vkh variables) must be equal to the number of allocations of clients ito facility jwith ij such that Cij is equal to c(h)(see Figure 3). Finally, constraints (30) are scheduling constraints based on costs values rather than in couples ij of client-facility (see Figure 4). k=1 k=2 k=3 k=4 ij 11 22 33 44 21 43 12 34 31 14 42 23 24 32 13 41 0 cij 0 0 0 1 1 2 2 3 4 4 5 5 6 7 9 c(h) 0 1 2 3 4 5 6 7 9 v12+v22+v32+v42 = x12+x34 Figure 3: The rationale of Constraints (29) k=1 k=2 k=3 k=4 h 0 2 3 4 5 6 7 81 v30+v31+v32+v23+v24+v25+v26+v27+v28 ≤ 1 Figure 4: The rationale of Constraints (30) 2.4. An aggregated formulation Here, we introduce another formulation (DOMP4) based on the aggregation of order constraints from DOMP3corresponding to the same position. It therefore requires a smaller number of constraints. (DOMP4) min n X i=1 n X j=1 n X k=1 λkCijxk ij s.t. (3),(4),(5),(6) n X i n X j     n X i0=1 n X j0=1: i0j0ij xk i0j0+ n X i0=1 n X j0=1: i0j0ij xk−1 i0j0    ≤n2k= 2,· · · , n (32) xk ij ∈ {0,1}i, j, k = 1, . . . , n yj∈ {0,1}j= 1, . . . , n 8 The new constraints (32), that we call weak order constraints, ensure that if a couple ij occupies the k-th position then in (k−1)-th position there must be a more preferred allocation. It is due to the coefficients of each variable in the inequality. In each constraint there are two different positions, kand k−1, so that, by (4), two variables must take value one and all the others will be equal to zero. If we do not take into account the variables taking the value zero and we assume that the variables with value one for positions k and k−1 are in position sand t, respectively, we have the following: (n2−(s−1))xk isjs+txk−1 itjt≤n2, which is valid if and only if t<s. Aggregating the scheduling constraints in DOMP3Cfor the different values of the costs, namely in h= 1, . . . , G, results in a new valid model. This model is the aggregated version of DOMP3Cthat we denote as DOMP4C: (DOMP4C) min n X k=1 G X h=0 λkc(h)vkh s.t. (11),(12),(13),(28),(29) G X h=1  X h0<h vkh0+X h0≥h vk−1h0 ≤G k = 2, . . . , n (33) xij, vkh ∈ {0,1}k= 1, . . . , n h= 0, . . . , G. Using the rationale of (7), since there is only one binary variable vkh in each position kby (28), the following constraints are valid inequalities for both formulations DOMP3Cand DOMP4C: G X h=1 c(h)vk−1h− G X h=1 c(h)vkh ≤0k= 2, . . . , n. (34) In fact, these constraints can also define another valid formulation for DOMP replacing (33) by (34) in the above formulation. In our experiments, we will not use this last possibility and instead, we shall use (34) as valid inequalities to strengthen DOMP3Cand DOMP4C. 3. Theoretical results In this section, we provide a theoretical comparison of the four formulations presented in Section 2 and some polyhedral results regarding our formulation DOMP3. Our goal is to state the formal relationships between the lower bounds provided by the linear relaxations of the considered formulations. In addition, we also give some families of tight valid inequalities which are proven to be facet of the polytope defined by assignment constraints (see Section 3.2.) 3.1. Comparison of formulations We denote by zl(·) the value of the objective function of DOMPlevaluated at the point (·), by Plthe polytope defining the feasible set of the linear relaxation of formulation DOMPl, and by PI lthe convex hull of the integer solutions within that polytope. Consider the following mapping f: [0,1]n3×[0,1]n−→ [0,1]n2×[0,1]n×[0,1]nG (xk ij, yj)7−→ (xij, yj, ukh) 9 (b) PI 3⊂g(PI 2) To prove the converse, we observe by that every integer solution in PI 3provides a projected integer solution and we have seen in Theorem 2 that this point satisfies the inequalities defining P2. 2. PI 1=g(PI 2) (a) g(PI 2)⊂PI 1 Every integer solution in P2is projected into an integer solution and, from Theorem 1 this point satisfies the inequalities defining P1. (b) PI 1⊂g(PI 2) Conversely, let (xk ij, yj)∈PI 1and its projection (xij , yj, ukh) by means of (35) and (36). It is easy to see that (3), (5) and (6) are equivalent to (12), (13) and (11), respectively. Since all the variables are positive the point (xk ij, yj) satisfies n X Cij ≥c(h) xk ij ≥ n X Cij ≥c(h+1) xk ij, and (14) holds. Furthermore, (16) holds by the change of variable defined in (35) and (36). In inequalities (7), we can order the cost without loss of generality. Next by (4), there will be an unique variable with value one in each position. Therefore we obtain the following constraint for each different cost, n X Cij ≥c(h) xk+1 ij ≥ n X Cij ≥c(h) xk ij h= 1, . . . , G, k = 1, . . . , n −1 and point (xij, yj, ukh) satisfies (15). Furthermore, this point is integer. 3. PI 3=PI 4 (a) PI 4⊂PI 3 Let (xk ij, yj)∈PI 4. In particular, the weak order constraint is satisfied. By (4) we have n X i0=1 :i0j0 n X j0=1 ij xk i0j0≤1 and n X i0=1 :i0j0 n X j0=1 ij xk−1 i0j0≤1, where both inequalities are a sum of binary variables. So there are two possibilities: n X i0=1 n X j0=1: i0j0ij xk i0j0+ n X i0=1 n X j0=1: i0j0ij xk−1 i0j0≤1,or (43) n X i0=1 n X j0=1: i0j0ij xk i0j0+ n X i0=1 n X j0=1: i0j0ij xk−1 i0j0= 2.(44) In case (44), there will be two different costs ˆıˆ ≺˜ i˜ j(note that ˆıˆ = ˜ i˜ jis not possible by (6)) such that xk ˆıˆ = 1 and xk−1 ˜ i˜ j= 1. This fact contradicts the ordering relationship. Thus, the only possibility is that case (43) holds and, in consequence, constraint (21) is satisfied. 16 (b) PI 3⊂PI 4 Let (xk ij, yj) be an integer point satisfying PI 3. This integer point belongs to PI 4because if for all kwe add constraints (21) for all iand jwe obtain (32). Corollary 3. PI 1=PI 3=PI 4=g(PI 2). 3.2. On the polytope defined by the assignment constraints The goal of this section is to provide some results about the facial structure of the polytope corresponding to the assignment constraints of formulation DOMP3. We restrict ourselves to the analysis of formulation DOMP3since, according to Corollary 2, it gives the tightest formulation to DOMP. Similar studies of (simpler) polytopes related to location problem have been carried out by e.g. Arbib et al. [1], Guignard [8], Cornu´ejols and Thizy [4], de Farias Jr. [5] and Vasilyev et al. [18]. For a given set J of popen facilities, we define the assignment polytope P3(J) of DOMP3as follows: X j∈J n X k=1 xk ij ≤1i= 1, . . . , n (45) n X i=1 X j∈J xk ij ≤1k= 1, . . . , n (46) n X k=1 xk ij ≤1i= 1, . . . , n, j ∈J(47) X i0j0ij xk−1 i0j0+X i0j0ij xk i0j0≤1i= 1, . . . , n, j ∈J, k = 2, . . . , n (48) xk ij ≥0i, k = 1, . . . , n, j ∈J(49) We show which of the constraints that describe this polytope are facet inducing. Clearly, this contributes to the good quality of the LP relaxation of DOMP3, since this assignment polytope represents the underlying structure of the problem once the set of open facilities is determined. We summarize the polyhedral properties of, PI 3(J), the convex hull of P3(J)∩{0,1}n2×pin the following result. Proposition 2. 1. dim(P3(J)) = n2p. 2. Constraints (45) and (49) induce facets of PI 3(J). Proof. Constraints (45)-(49) define a particular packing polytope which has been studied by Padberg [15] among others. The results are consequences of this observation and the fact that variables appearing in (45) define maximal cliques in the conflict graph associated to the problem, see Padberg [15]. Observe that constraints (47) do not induce facets of P3(J) since they are dominated by constraints (45). Let us denote by (ij)s, s = 1, . . . , n2the couple client ifacility jwhere Cij is the s-th lower cost over the cost matrix C. For instance, (ij)1and (ij)n2are, respectively, the most and the least preferred couples in the sorted list of costs of the cost matrix C. Definition 1. We call a pair ij generic if for some feasible set Jsuch that j∈J, it satisfies 1. Let ˆıˆ be the couple client ˆı facility ˆ ∈Jsuch that it does not exist any other couple client i0facility j0∈Jsatisfying ˆıˆ ≺i0j0≺ij, i.e. ˆıˆ is the couple immediately before ij. Then, ˆı 6=i. 17 2. Let ˜ı˜ be the couple client ˜ı facility ˜ ∈Jsuch that it does not exist any other couple client i0facility j0∈Jsatisfying ij ≺i0j0≺˜ı˜, i.e. ˆıˆ is the couple immediately after ij. Then, ˜ı 6=i. Intuitively, a pair ij is generic with respect to a feasible solution set Jif the remaining feasible allocation costs are well distributed around it. That is, there are costs of different clients surrounding the ij cost (Cij) in the sorted list of costs. Proposition 3. If ij is generic for the feasible set Jthen X i0j0ij xk−1 i0j0+X i0j0ij xk i0j0+ n X l=k+1 xl (ij)1+ k−2 X l=1 xl (ij)n2≤1k= 2,· · · , n (50) is facet defining for the assignment polytope PI 3(J). ijîĵ ĩǰ 3 4 2 1 (ij) (ij) 1 n2 Figure 5: Constraints (50)’s scheme. Proof. According to Definition 1, let us denote by ˆıˆ and ˜ı˜ the couples immediately before and after the couple ij in the sorted list of costs (see Figure 5). Let us prove that this family of constraints are facet defining inequalities showing that they are maximal cliques, i.e. it does not exist any variable which does not belong to the clique and that it is incompatible with all of those that appear in it. Next, we show that for each variable that does not belong to the considered clique, there is one in the clique compatible with it. 1. For all i0j0≺ij,l≤k−1 if i06=i, variables xl i0j0and xk ij are compatible. Otherwise, variables xl i0j0 and xk ˆıˆ are compatible. (Case 1 in Figure 5.) 2. For all i0j0ij,l < k −1 and i0j06= (ij)n2if i06=ivariables xl i0j0and xk ij are compatible. Otherwise, variables xl i0j0and xk ˆıˆ are compatible. (Case 2 in Figure 5.) 3. For all i0j0ij,l=k+ 1, .., .n and i0j06= (ij)1if i06=i, variables xl i0j0and xk−1 ij are compatible. Otherwise, variables xl i0j0and xk−1 ˜ı˜ are compatible. (Case 3 in Figure 5.) 4. For all i0j0ij,l=k, ..., n if i06=ivariables xl i0j0and xk−1 ij are compatible. Otherwise, variables xl i0j0 and xk−1 ˜ı˜ are compatible. (Case 4 in Figure 5.) Remark that in our computational experiments, we did not reinforce constraints (48) because our preprocessing procedure (see Section 4) sets to zero all the variables appearing in these reinforcements associated with the costs of couples (ij)1and (ij)n2. Thus, in most cases in our formulation DOMP3, after preprocessing those constraints are already facet inducing. 18 Notation λ-vector Name T1 (1,1,...,1,1) p-median T2 (0,0,...,0,1) p-center T3 (0,0,...,0,0,1,1,...,1,1 | {z } k )k-centrum T4 (0,0,...,0,0 | {z } k1 ,1,1,...,1,1,0,0,...,0,0 | {z } k2 ) (k1+k2)-trimmed mean T5 (0,1,0,1,0,1,0,1, . . . ) – T6 (. . . , 0,0,1,0,0,1) – T7 λrandom Random T8 (α, α, . . . , α, α, 1) Centdian Table 3: Types of λ-vectors used in experiments 4. Computational study In order to test the performance of our new formulations for DOMP, we have performed intensive computational tests comparing results with respect to previous available formulations of DOMP (see Boland et al. [3], Dom´ınguez-Mar´ın [6], Nickel [13], Nickel and Puerto [14]) and Mar´ın et al. [11]). 4.1. Description of the test instances We use two different types of instances. First, we consider random instances in which the elements of the cost matrix are integer numbers randomly generated between 10000 and 100000. The second set of instances consists in p-median instances from OR Lib, Beasley [2]. Regarding the random instances; we vary the number of clients nin {10,20,30,40,50}and for each n, we consider three possible values for the number of facilities to be open: p=n 4,n 3,n 2. As for the p−median instances from Beasley’s library, we have selected graphs corresponding to p− med1, . . . , p−med20 with up to 400 nodes from the original data. Each set of nodes is divided in two disjoint subsets containing, respectively, the set of clients and the set of facilities. (The reader may observe that in these instances we force these two sets to be disjoint). Next, each cost Cij is computed as the shortest path between client iand facility jin the resulting complete graph induced by the above described set of nodes. Eight different types of λ-vectors are tested. Their description is provided in Table 3. We consider, among others, p-median, p-center, p-k-centrum, p-trimmed mean, random and p-α-centdian problems. In the k-centrum case, k=n 2. In the (k1+k2)-trimmed mean, k1 = k2 = n 10 . When λis taken as a random vector we generate 5 instances which contain values randomly drawn between 1 and 100. Finally, we use α= 0.5 in the centdian case. For each type of data and each possible values of parameters n,pand vector λ, the results presented consist in average value over five instances. This results in an overall number of 900 tested instances. 4.2. Preprocessing We use two different preprocessings to set some assignment variables to zero, the first one is based on feasibility and the second on optimality. Claim 1. (Feasibility based preprocessing) 1. Let l(h) = |{i:Cij ≤c(h), i, j = 1, . . . , n}| and u(h) = |{i:Cij ≥c(h), i, j = 1, . . . , n}|. Then, (a) vkh = 0 and ukh = 1 k=l(h)+1, . . . , n, h = 1 . . . , G.; and (b) vkh = 0 and ukh = 0 k= 1, . . . , n −u(h), h = 1 . . . , G. 2. Let l(ij) = |{i0: minˆ i0ˆ ≺ij, i06=i}| and u(ij) = |{i0: maxˆ i0ˆ ij, i06=i}|. Then, (a) xk ij = 0 for all i, j = 1, . . . , n, k =l(ij)+2, . . . , n; and 19 (b) xk ij = 0 for all i, j = 1, . . . , n, k = 1, . . . , n −u(ij)−1. This claim formalizes the fact that a given cost Cij can appear in position kfeasible solution only if there are at least k−1 allocations costs lower than or equal to and n−kgreater than or equal to Cij. Claim 2. Let ij be a couple client facility. If |{j0:Cij > Cij0}| > n −p, then xij = 0 and xk ij = 0 for k= 1, . . . , n. This second claim removes some feasible solutions of the problems which cannot be optimal because they are dominated for other feasible solutions with smaller objective value. Table 4 shows the percentage of variables vkh,ukh,xij and xk ij fixed by our preprocessing based on Claim 1 and 2. Notice that the percentage of two or three index variables fixed to zero in random instances is almost equal because assignment costs in random instances are almost all different. Furthermore, we observe that the percentage of fixed variables slightly decreases with n. Observe that Beasley instances have a larger number of ties in the distances (allocation costs). For this reason, these instances are solved using our compact formulations DOMP3Cand DOMP4C, whereas those with random data are solved with the formulations that do not take adventage of ties, namely DOMP3and DOMP4. Hence, Beasley instances are preprocessed with Claim 1.1 (for vvariables) and Random instances with Claim 1.2 (for xk ij variables). n (Random instances) 10 20 30 40 50 Average Claim 1.2 16.40% 9.17% 6.21% 4.67% 3.93% 8.07% Claim 2 20.00% 25.00% 30.00% 29.99% 29.99% 27.00% Total 29.68% 29.61% 30.43% 32.28% 31.97% 31.40% n (Beasley instances) 50 100 150 200 Average Claim 1.1 27.03% 31.86% 33.42% 31.01% 30.83% Claim 2 28.91% 26.85% 25.87% 25.26% 26.72% Total 27.38% 30.33% 29.67% 27.12% 28.62% Table 4: Number of variables fixed by preprocessing 4.3. Computational results All our experiments have been carried out on a PC with two Intel Xeon processors with 3.46 GHz and 48 GB of RAM. The models were written in Mosel and solved using Xpress IVE 7.3, To have a clean comparison of our solution approaches, all automatic cuts from Xpress have been disabled. We now report a summary of our computational experiments. Detailed information can be found in the Supplementary material included in the Appendix of this paper. In particular, we report results for the different types of lambda vectors from Table 3. 4.3.1. Random allocation costs data sets. Table 5 provides a comparison of the LP-relaxations of models DOMP1,DOMP2,DOMP3,DOMP4 and DOMP4∩1, averaging for n= 20 and n= 40 and all possible values of pand λ. The last model DOMP4∩1consists in DOMP4to which constraints (7) of DOMP1have been appended. Specifically, we report the integrality gap defined as GAP =z∗−zLP l z∗, where z∗and zLP lrepresent the optimal value of DOMPland its LP-relaxation, respectively. We point out that the reported GAP in all formulations is computed after the application of the preprocessing results developed in this paper. We observe that the values of the integrality gap vary between 2.29% and 7.34%, obtained in formulations DOMP3the best and DOMP2or DOMP1the worst. In all cases, the LP gaps are good but specially in formulation DOMP3 which seems to be rather tight. We remark also the small difference that is obtained adding constraints (7) to the formulation DOMP4in terms of integrality gap. Moreover, we point out that we have obtained 20 the same LP gap with formulations DOMP1and DOMP2for all the tested instances. (It is still an open question whether this is also theoretically true.) Thus, from Table 5 one could conclude that the best formulation is DOMP3. In spite of that, the large number of inequalities (O(n3)) used in the model makes it rather slow whenever the number of clients nis of moderate size (n > 50). In order to define a solution approach which presents the best performance, we have conducted a preliminary computational test with instance sizes n= 10,20,30. Our first strategy consists in solving DOMP3with a pure branch-and-bound. Our second strategy, DOMP4∩1(B&B), solves DOMP4∩1with a pure branchand-bound. The third approach, DOMP4∩1(B&C−3), starts by solving the LP-relaxation of DOMP4 and then adds inequalities (7) at the root node and order constraints (21) as long as they are violated by the current solution of the LP. The reader may observe that both families of inequalities are cliques in the conflict graph induced by the three index variables of our formulation. Therefore, order constraints could in principle be added by standard clique cuts generation techniques implemented in Xpress. Nevertheless, our own implementation is more efficient since the separation of the entire family of valid inequalities (21) can be performed in O(n3) by sequentially updating the L.H.S. value of the order constraints when switching from a couple ij to the adjacent one in the same position k. The fourth and last strategy DOMP4(B&C−3) is a branch and cut algorithm based upon DOMP4adding only valid inequalities from (21). Formulation GAP DOMP17,34% DOMP27,34% DOMP32,29% DOMP46,27% DOMP4∩16,06% Table 5: Average integrality gaps Solution approach Time (s) #nodes DOMP3(B&B) 463.24(8) 46.20 DOMP4∩1(B&B) 18.64 1389.51 DOMP4∩1(B&C−3) 52.25 24.94 DOMP4(B&C−3) 39.39 29.37 Table 6: CPU-Time and Number of nodes of the different formulations for n= 10,20,30. Our results are reported in Table 6. There we have included the CPU-times and the number of nodes in the B&B tree for solving instances to optimality within 2 hours of CPU-time. The numbers between parentheses indicate the number of unsolved instances within the time limit. We observe that on average DOMP4∩1(B&B) is the strategy that solves problems faster even though it has to visit the largest number of nodes in the B&B tree. This is explained by the fact that it is the most compact formulation (with the smallest number of inequalities) and therefore, it can be easily solved at each node. On the other hand, DOMP3(B&B) is the heaviest one (in terms of LP representation) giving rise to worse CPU-times although it visits few nodes in the searching phase. In between, we found the two branch-and-cut procedures that we have tested DOMP4∩1(B&C−3) and DOMP4(B&C−3). In the implementation of this two B&C approaches we have tested to separate maximal clique inequalities over the conflict graph as an alternative to our own separation procedure. Nevertheless, our separation algorithm applied on inequalities (7) and (21) gives better results. From the above two tables, we can conclude that the best strategies to be tested in the intensive computational tests are DOMP4∩1(B&B) and, at times, DOMP4(B&C−3). To finish, Table 7 allows us to compare our best strategy, i.e. DOMP4∩1(B&B), with a branch-andbound approach based on DOMP2, as well as to determine the size of instances that can be solved within a reasonable time limit of two hours. This table is organized in four columns. The first two columns show the size nand pof the instances. The last two columns show the average CPU-time in seconds necessary for solving those instances applying formulation DOMP2(B&B) and DOMP4∩1(B&B). The numbers between parentheses indicate the number of unsolved instances within the time limit of 2 hours. From Table 7 we remark that DOMP4∩1(B&B) performs similarly as DOMP2(B&B) (it improves the behavior of DOMP2(B&B) only for some instance sizes). It is better for data instances with n < 40 in all combinations of p. In addition, for larger n, i.e. n= 40,50, it is also better except if pis relatively small as compared with n. Furthermore, we also observe that for n= 50 both models fail to solve some instances for the smallest tested value of p= 12. Finally, Table 8 allows us to conclude that three index with aggregated 21 scheduling constraints performs the best whenever the number of facilities to be located is not too small compared to the number of possible locations. n p Time (# unsolved) DOMP2(B&B)DOMP4∩1(B&B) 10 2 1.12 0.42 10 3 0.54 0.25 10 5 0.18 0.09 20 5 22.21 5.80 20 6 9.05 4.09 20 10 3.33 1.54 30 7 161.32 121.04 30 10 66.85 22.70 30 15 33.15 11.87 40 10 449.11 625.68 40 13 232.79 174.62 40 20 116.99 49.13 50 12 1870.80(2) 3184.54(6) 50 16 988.08 804.36 50 25 771.51(2) 157.44 Average 315.14 344.24 Table 7: Summary of results with random matrices pDOMP2(B&B)DOMP4∩1(B&B) n 4500.91(2) 787.49(6) n 3259.46 201.20 n 2185.03(2) 44.01 Table 8: CPU-Time of the different formulations for different values of p. 4.3.2. Beasley’s data set. The large number of ties within the cost matrices of this data set suggests that on this second part of the study DOMP4Cis the appropriate formulation to solve the problems. For each cost matrix we solve each instance with DOMP2(B&B) and DOMP4C(B&C−3), i.e. formulation DOMP4Cwithin a branch and cut scheme separating inequalities (30). For each value of nwe solve those problems for the number of open facilities p, suggested in the original data from Beasley’s library, and for all the considered vectors of λshown in Table 3. Table 9 shows the average results for these instances. Detailed information for each λcan be found in the electronic appendix. This table is organized in five columns. The first three columns show the name of the instance problem and its size nand p. The last two columns show the CPU-time for solving those instances applying strategies DOMP2(B&B) and DOMP4C(B&C−3). The numbers between parentheses indicate the number of instances that could not be solved to optimality within the time limit of 2 hours. We can see that in 11 out of 20 instances DOMP4C(B&C−3) is faster than DOMP2(B&B). This behavior confirms that both formulations have a rather similar performance. In spite of that, we observe that the use of our new formulation outperforms DOMP2provided that nis of moderate size n < 50 or whenever the size of prelative to nis not too small, namely p/n ≥0.2. This behavior allows us to conclude that DOMP4C(B&C−3) is advisable to be used, at least, in those cases. 22 Problem n p Time (#unsolved) DOMP2(B&B)DOMP4C(B&C−3) pmed1 50 5 94.65 94.66 pmed2 50 10 79.24 40.58 pmed3 50 10 103.84 37.58 pmed4 50 20 34.94 8.57 pmed5 50 33 20.44 3.90 pmed6 100 5 805.02 4757.27(5) pmed7 100 10 617.60 2708.85(2) pmed8 100 20 848.95(1) 543.26 pmed9 100 40 130.28 45.87 pmed10 100 67 45.84 22.01 pmed11 150 5 1353.85(1) 3702.40(2) pmed12 150 10 1667.06(1) 4235.25(4) pmed13 150 30 2030.82(1) 3808.45(3) pmed14 150 60 509.60 302.93 pmed15 150 100 88.48 52.17 pmed16 200 5 2760.99(1) 5940.66(6) pmed17 200 10 2940.22(1) 5670.05(6) pmed18 200 40 1830.27(1) 4129.70(3) pmed19 200 80 358.11 264.35 pmed20 200 133 256.21 150.76 Table 9: Summary of results using Beasley’s data set 5. Concluding remarks This paper presents new formulations for the Discrete Ordered Median Problem based on order constraints (21) that are valid for the general non free self-service case. Furthermore, we prove theoretical relationships, in terms of their LP-gap, for different formulations of DOMP. According to the theoretical and computational results obtained in this paper, the main quality of the new formulations is that they provide substantial improvement of the integrality gap with respect to previously known ones. We have observed that the LP-gap of DOMP1and DOMP2is always equal. It is currently an open question whether this property holds in general. This question will be a subject of our future research. This paper has also opened another interesting line of research that consists in finding extensions of some of the existing formulations to exploit special structures of the lambda coefficients, as for instance the one in Mar´ın et al. [11]. Extensions based on the results in [11] seem to require additional variables to handle the non free self-service case. A similar rationale can be also applied to the new formulations in this paper. Further theoretical and computational comparisons of the above mentioned new approach will be the subject of a follow up paper. Acknowledgement The research of the first author was partially supported by the Interuniversity Attraction Poles Programme initiated by the Belgian Science Policy Office and Spanish MTM2012-36163-C06-04 . The research of the second and third authors was partially supported by the project FQM-5849 (Junta de Andaluc´ıa\FEDER) and MTM2010-19576-C02-01 (MICINN, Spain). [1] C. Arbib, M. Labb´e, and M. Servilio. Scheduling two chains of unit jobs on one machine: A polyhedral study. Networks, 58-2:103–113, 2011. [2] J.E. Beasley. OR-Library, 2012. people.brunel.ac.uk/~mastjjb/jeb/info.html. [3] N. Boland, P. Dom´ınguez-Mar´ın, S. Nickel, and J. Puerto. Exact procedures for solving the discrete ordered median problem. Computers & Operations Research, 33(11):3270–3300, 2006. ISSN 0305-0548. [4] G. Cornu´ejols and J.M. Thizy. Some facets of the simple plant location polytope. Mathematical Programming, (23):50–74, 1982. 23 [5] I. R. de Farias Jr. A family of facets for the uncapacitated p-median polytope. Operations Research Letters, 28(4):161–167, 2001. [6] P. Dom´ınguez-Mar´ın. The Discrete Ordered Median Problem: Models and Solution Methods. Kluwer, 2003. [7] E. Fern´andez, J. Puerto, A.M. Rodr´ıguez-Ch´ıa. On Discrete Optimization with Ordering. Annals of Operations Research, 207(1):83–96, 2013. [8] M. Guignard. Fractional vertices cuts and facets of the simple plant location problem. Mathematical Programming Study, (12):152–160, 1980. [9] J. Kalcsics, S. Nickel, J. Puerto, A. M. Rodr´ıguez-Ch´ıa, Distribution Systems Design With Role Dependent Objectives EJOR, 202 (2010) 491-501. [10] A. Mar´ın, S. Nickel, J. Puerto, and S. Velten. A flexible model and efficient solution strategies for discrete location problems. Discrete Applied Mathematics, 157(5):1128–1145, 2009. ISSN 0166-218X. [11] A. Mar´ın, S. Nickel, and S. Velten. An extended covering model for flexible discrete and equity location problems. Mathematical Methods of Operations Research, 71 (1): 125–163, 2010. [12] G.L. Nemhauser and L.E. Trotter. Vertex packings: Structural properties and algorithms. Mathematical Programming, 8:232–248, 1975. [13] S. Nickel. Discrete ordered weber problems. In Operations Research Proceedings 2000, pages 71–76. Springer Verlag, 2001. [14] S. Nickel and J. Puerto. Location Theory: A Unified Approach. Springer Verlag, 2005. [15] M. W. Padberg. On the Facial Structure of Set Packing Polyhedra. Mathematical Programming, 5:199–215, 1973. [16] J. Puerto. A new formulation of the capacitated discrete ordered median problem with {0,1}-assignment. In Operations Research Proceedings 2007, volume 1, pages 165–170. Springer, 2008. ISBN 978-3-540-77902-5. [17] J. Puerto, A. B. Ramos, and A. M. Rodr´ıguez-Ch´ıa. Single-allocation ordered median hub location problems. Comput. Oper. Res., 38(2):559–570, February 2011. ISSN 0305-0548. [18] I. Vasilyev, X. Klimentova, and M. Boccia. Polyhedral study of simple plant location problem with order. Operations Research Letters, (2):153–158, 2013. 24 AppendixA. Supplementary material for the paper “A Comparative Study of Formulations and Solution Methods for the Discrete Ordered p-Median Problem” This electronic supplement reports detailed integrality GAP and CPU times for each of the eight type of λ-vector that have been considered in the computational results of the paper. We follow the same notation so that T1, ..., T8, refer respectively to p-median, p-center, p-k-centrum, (k1+k2)-trimmed mean, the sequence of 0,1,0,1..., the sequence of . . . , 0,0,1,0,0,1, random elements and p-α-centdian (see Table 3 in the paper). AppendixA.1. GAP In this section, we detail the result of Table 5 in the paper comparing the average integrality gap of the formulations in the paper for the different types of λ-vectors T1, . . . ,T8. Table A.10: Average integrality GAP Formulation λ-vector DOMP1DOMP2DOMP3DOMP4DOMP4∩1 T1 0,14% 0,14% 0,14% 0,14% 0,14% T2 36,04% 36,04% 18,81% 31,59% 31,57% T3 16,38% 16,38% 4,53% 13,07% 13,07% T4 5,53% 5,53% 0,14% 3,29% 3,29% T5 2,16% 2,16% 0,91% 3,42% 2,16% T6 4,01% 4,01% 1,78% 5,38% 3,99% T7 3,32% 3,32% 0,77% 3,26% 3,00% T8 3,69% 3,69% 1,02% 3,30% 3,30% 25