The $p$-Median Problem

Summary

The goal of the $p$-median problem is to locating $p$ facilities to minimize the demand weighted average distance between demand nodes and the nearest of the selected facilities. Hakimi (1964, 1965) first considered this problem for the design of network switch centers. However, this problem has been used to model a wide range of applications, such as warehouse location, depot location, school districting and sensor placement.

Problem Statement

The $p$-median problem can be formulated mathematically as an integer programming problem using the following model.

Sets

$M$ = set of candidate locations
$N$ = set of customer demand nodes

Parameters

$p$ = number of facilities to locate
$d_j$ = demand of customer $j$, $\forall j \in N$
$c_{ij}$ = unit cost of satisfying customer $j$ from facility $i$, $\forall i \in M, \forall j \in N$

Variables

$x_{ij}$ = fraction of the demand of customer $j$ that is supplied by facility $i$, $\forall i \in M, \forall j \in N$
$y_i$ = a binary value that is $1$ is a facility is located at location $i$, $\forall i \in M$

Objective

Minimize the demand-weighted total cost
$\min \sum_{i \in M} \sum_{j \in N} d_j c_{ij} x_{ij}$

Constraints

All of the demand for customer $j$ must be satisfied
$\sum_{i \in M} x_{ij} = 1$, $\forall j \in N$

Exactly $p$ facilities are located
$\sum_{i \in M} y_i = p$

Demand nodes can only be assigned to open facilities
$x_{ij} \leq y_i$, $\forall i \in M, \forall j \in N$

The assignment variables must be non-negative
$x_{ij} \geq 0$, $\forall i \in M, \forall j \in N$

Pyomo Formulation

The following is an abstract Pyomo model for this problem:



In [1]:

    
!cat p-median.py









    



from pyomo.environ import *
import random

random.seed(1000)

model = AbstractModel()

# Number of candidate locations
model.m = Param(within=PositiveIntegers)
# Number of customers
model.n = Param(within=PositiveIntegers)
# Set of candidate locations
model.M = RangeSet(1,model.m)
# Set of customer nodes
model.N = RangeSet(1,model.n)

# Number of facilities
model.p = Param(within=RangeSet(1,model.n))
# d[j] - demand of customer j
model.d = Param(model.N, default=1.0)
# c[i,j] - unit cost of satisfying customer j from facility i
model.c = Param(model.M, model.N, initialize=lambda i, j, model : random.uniform(1.0,2.0), within=Reals)

# x[i,j] - fraction of the demand of customer j that is supplied by facility i
model.x = Var(model.M, model.N, bounds=(0.0,1.0))
# y[i] - a binary value that is 1 is a facility is located at location i
model.y = Var(model.M, within=Binary)

# Minimize the demand-weighted total cost
def cost_(model):
    return sum(model.d[j]*model.c[i,j]*model.x[i,j] for i in model.M for j in model.N)
model.cost = Objective(rule=cost_)

# All of the demand for customer j must be satisfied
def demand_(model, j):
    return sum(model.x[i,j] for i in model.M) == 1.0
model.demand = Constraint(model.N, rule=demand_)

# Exactly p facilities are located
def facilities_(model):
    return sum(model.y[i] for i in model.M) == model.p
model.facilities = Constraint(rule=facilities_)

# Demand nodes can only be assigned to open facilities 
def openfac_(model, i, j):
    return model.x[i,j] <= model.y[i]
model.openfac = Constraint(model.M, model.N, rule=openfac_)

This model is simplified in several respects. First, the candidate locations and customer locations are treated as numeric ranges. Second, the demand values, $d_j$ are initialized with a default value of $1$. Finally, the cost values, $c_{ij}$ are randomly assigned.

Model Data

This model is parameterized by three values: the number of facility locations, the number of customers, and the number of facilities. For example:



In [2]:

    
!cat p-median.dat









    



param m := 10;
param n := 6;
param p := 3;

Solution

Pyomo includes a pyomo command that automates the construction and optimization of models. The GLPK solver can be used in this simple example:



In [3]:

    
!pyomo solve --solver=glpk p-median.py p-median.dat









    



[    0.00] Setting up Pyomo environment
[    0.00] Applying Pyomo preprocessing actions
[    0.00] Creating model
[    0.02] Applying solver
[    0.06] Processing results
    Number of solutions: 1
    Solution Information
      Gap: 0.0
      Status: optimal
      Function Value: 6.431184939357673
    Solver results file: results.json
[    0.07] Applying Pyomo postprocessing actions
[    0.07] Pyomo Finished

By default, the optimization results are stored in the file results.yml:



In [4]:

    
!cat results.yml









    



# ==========================================================
# = Solver Results                                         =
# ==========================================================
# ----------------------------------------------------------
#   Problem Information
# ----------------------------------------------------------
Problem: 
- Name: unknown
  Lower bound: 6.43118493936
  Upper bound: 6.43118493936
  Number of objectives: 1
  Number of constraints: 68
  Number of variables: 71
  Number of nonzeros: 191
  Sense: minimize
# ----------------------------------------------------------
#   Solver Information
# ----------------------------------------------------------
Solver: 
- Status: ok
  Termination condition: optimal
  Statistics: 
    Branch and bound: 
      Number of bounded subproblems: 1
      Number of created subproblems: 1
  Error rc: 0
  Time: 0.0117330551147
# ----------------------------------------------------------
#   Solution Information
# ----------------------------------------------------------
Solution: 
- number of solutions: 1
  number of solutions displayed: 1
- Gap: 0.0
  Status: optimal
  Message: None
  Objective:
    cost:
      Value: 6.43118493936
  Variable:
    x[6,5]:
      Value: 1
    y[3]:
      Value: 1
    x[6,2]:
      Value: 1
    x[9,6]:
      Value: 1
    y[9]:
      Value: 1
    x[3,4]:
      Value: 1
    y[6]:
      Value: 1
    x[6,3]:
      Value: 1
    x[6,1]:
      Value: 1
  Constraint: No values

This solution places facilities at locations 3, 6 and 9. Facility 3 meets the demand of customer 4, facility 6 meets the demand of customers 1, 2, 3 and 5, and facility 9 meets the demand of customer 6.

References

S.L. Hakimi (1964) Optimum location of switching centers and the absolute centers and medians of a graph. Oper Res 12:450–459
S.L. Hakimi (1965) Optimum distribution of switching centers in a communication network and some related graph theoretic problems. Oper Res 13:462–475