# Symmetric undirected $p$-median instance with fractional LP solution?

The $$p$$-median problem is NP-hard, so its LP relaxation does not naturally have all-integer solutions. However, it very often does; in fact, it can be hard to find an instance for which the LP relaxation solution is not integer.

The only examples that I know of for which the LP solution is non-integer use asymmetric distances (or, equivalently, directed graphs)—for example, the examples by ReVelle and Swain (1970) and Baïou and Barahona (2011).

Is there a symmetric, undirected $$p$$-median instance for which the optimal solution to the LP relaxation is not integer? I'm especially looking for small examples (say, 5 nodes).

• Which formulation do you have in mind? Aug 9 '19 at 2:25
• The standard one (at least I’d consider it standard). For example see (2.1)-(2.6) of Daskin and Maass 2015. Aug 9 '19 at 2:30

I think I've found an instance with four nodes and $$p = 2$$ via brute force (a lot of randomized instances). I've attached my Python script as well. I relaxed the Daskin and Maass (2015) formulation and assumed $$I = J$$.

Nodes: $$I = J = \{1, 2, 3, 4\}$$

Demands: $$d = (75, 34, 40, 40)$$

Costs (or distances):

$$c = \begin{bmatrix} 0 & 39 & 95 & 7 \\ 39 & 0 & 83 & 18 \\ 95 & 83 & 0 & 16 \\ 7 & 18 & 16 & 0\\ \end{bmatrix}$$

Gurobi gives the following solution:

$$y = (0.5, 0.5, 0.5, 0.5 )$$

$$x = \begin{bmatrix} 0.5 & 0 & 0 & 0.5 \\ 0 & 0.5 & 0 & 0 \\ 0 & 0 & 0.5 & 0 \\ 0.5 & 0.5 & 0.5 & 0.5\\ \end{bmatrix}$$

The optimal demand-weighted cost is $$1028.5$$.

Code:

# Edited to add code highlighting
from gurobipy import *
import numpy as np
from random import *

## DEFINE MODEL ###

def p_median_relaxed(I, J, c, d, p):
model = Model("relaxed")

model.setParam('OutputFlag', 0)

######### VARIABLES AND BOUNDS #########
x = {}
y = {} # the (binary) facility location variable being relaxed

for i in I:
for j in J:
x[i, j] = model.addVar(vtype = GRB.CONTINUOUS, lb = 0)

for i in I:
y[i] = model.addVar(vtype = GRB.CONTINUOUS, lb = 0, ub = 1)

######### CONSTRAINTS #########

# (2.2)
for j in J:
model.addConstr(quicksum(x[i, j] for i in I) == 1, name = "c1_{}".format(j))

# (2.3)
model.addConstr(quicksum(y[i] for i in I) == p, name = "c2")

# (2.4)
for i in I:
for j in J:
model.addConstr(x[i, j] - y[i] <= 0, name = "c3_{}_{}".format(i, j))

######### OBJECTIVE ##########

model.setObjective(quicksum((d[j] * c[i, j] * x[i, j]) for i in I for j in J), GRB.MINIMIZE)

############ SOLVE ###########

model.optimize()

### STORE AND RETURN SOLUTION ###

x_val = {}
y_val = {}

for i in I:
y_val[i] = y[i].x

for j in J:
x_val[i, j] = x[i, j].x

return(x_val, y_val, model.objVal)

## CREATE RANDOMIZED INSTANCES AND SOLVE ##

siz = 4 # size of network

I = range(1, 5)
J = range(1, 5)
c = {}
d = {}

found_flag = 0

for ii in range(1000):
p = randrange(2, siz)

for j in J:
d[j] = randrange(0, 100, 1)

for i in I:

if i == j:
c[i, j] = 0
elif i < j:
c[i, j] = randrange(0, 100, 1)

for i in I:
for j in J:
if i > j:
c[i, j] = c[j, i]

(x_val, y_val, cost) = p_median_relaxed(I, J, c, d, p)

for i in I:
if y_val[i] <= 0.999:
if y_val[i] >= 0.001:
print(x_val)
print(y_val)
print(c)
print(d)
print(p)
print(cost)
print(' ')

found_flag = 1
break

if found_flag == 1:
break


I put this together quickly, so please let me know if I've made any errors or incorrect assumptions.

• Yep, this is perfect. I should have tried your random approach. Thanks! Aug 10 '19 at 2:29