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.