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

Question

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).

Answer 1

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:

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.