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.
The $p$-median problem can be formulated mathematically as an integer programming problem using the following model.
$M$ = set of candidate locations
$N$ = set of customer demand nodes
$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$
$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$
Minimize the demand-weighted total cost
$\min \sum_{i \in M} \sum_{j \in N} d_j c_{ij} x_{ij}$
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$
The following is an abstract Pyomo model for this problem:
!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.
This model is parameterized by three values: the number of facility locations, the number of customers, and the number of facilities. For example:
!cat p-median.dat
param m := 10; param n := 6; param p := 3;
Pyomo includes a pyomo
command that automates the construction and optimization of models. The GLPK solver can be used in this simple example:
!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
:
!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.