Note
Go to the end to download the full example code.
P-median problem#
We demonstrate how to implement the model for the p-median problem with DOcplex, gurobipy, Xpress, and highspy along with the functionality from opti-extensions.
Implementation reference: ampl/colab.ampl.com — Copyright (c) 2022-2022, AMPL Optimization inc. (licensed under the MIT License)
Mathematical formulation#
Index-sets:
- Customers to be served:\(i \in CUST\)
- Facilities to consider:\(j \in FAC\)
Parameters:
- Cost of supplying customer \(i\) from facility \(j\):\(cost_{i, j} \in \mathbb{R}_{0}^{+} \quad \forall \; i \in CUST \; \& \; \forall \; j \in FAC\)
- Number of facilities to be opened:\(p \in \mathbb{R}^{+}\)
Variables:
- Whether to serve customer \(i\) from facility \(j\):\(x_{i, j} \in \mathbb{B} \quad \forall \; i \in CUST \; \& \; j \in FAC\)
- Whether to open facility \(j\) or not:\(y_{j} \in \mathbb{B} \quad \forall \; j \in FAC\)
Objective:
- Minimize the total cost of serving all customers:\(\min \; \sum_{i \in CUST} \sum_{j \in FAC} cost_{i, j} \times x_{i, j}\)
Constraints:
- Each customer must be assigned to exactly one facility:\(\sum_{j \in FAC} x_{i, j} = 1, \; \forall \; i \in CUST\)
- Ensure that facility \(j\) is open if it serves customer \(i\):\(x_{i, j} \leq y_{j}, \; \forall \; i \in CUST \; \& \; \forall \; j \in FAC\)
- \(p\) facilities must be opened:\(\sum_{j \in FAC} y_{j} = p, \; \forall \; j \in FAC\)
Generate input data#
import math
import random
from typing import NamedTuple
random.seed(27)
class Instance(NamedTuple):
customers: list[int]
facilities: list[int]
cost: dict[tuple[int, int], float]
p: int
def generate_instance(num_customers: int, num_facilities: int, p: int):
customers = list(range(1, num_customers + 1))
customer_coord = {i: (random.uniform(0, 100), random.uniform(0, 100)) for i in customers}
facilities = list(range(num_customers, num_facilities + num_customers + 1))
facility_coord = {j: (random.uniform(0, 100), random.uniform(0, 100)) for j in facilities}
cost = {
# Eucledian distance between customer `i` and facility `j`
(i, j): round(math.sqrt((i_coord[0] - j_coord[0]) ** 2 + (i_coord[1] - j_coord[1]) ** 2), 2)
for i, i_coord in customer_coord.items()
for j, j_coord in facility_coord.items()
}
return Instance(customers, facilities, cost, p)
instance = generate_instance(5, 5, 2)
Set up problem data#
Index-sets#
from opti_extensions import IndexSet1D, IndexSetND
CUST = IndexSet1D(instance.customers, name='CUST')
FAC = IndexSet1D(instance.facilities, name='FAC')
CUST_x_FAC = IndexSetND(CUST, FAC, names=('CUST', 'FAC'))
Parameters#
from opti_extensions import ParamDictND
cost = ParamDictND(instance.cost, key_names=('CUST', 'FAC'), value_name='COST')
p = instance.p
# fmt: off
Implement with DOcplex#
from docplex.mp.model import Model
from opti_extensions.docplex import add_variables, solve
# Instantiate model
model = Model(name='p-median')
# Add variables
x = add_variables(model, indexset=CUST_x_FAC, vartype='binary', name='x')
# Instead of:
# x = model.binary_var_dict(CUST_x_FAC, name='x')
y = add_variables(model, indexset=FAC, vartype='binary', name='y')
# Instead of:
# y = model.binary_var_dict(FAC, name='y')
# Set objective
model.minimize(
cost @ x
# Instead of:
# model.sum(cost[i, j] * x[i, j] for i in CUST for j in FAC)
)
# Add constraints
model.add_constraints_(
x.sum(i, '*') == 1
# Instead of:
# model.sum(x[i, j] for j in FAC) == 1
for i in CUST
)
model.add_constraints_(
x[i, j] <= y[j]
for i in CUST for j in FAC
)
model.add_constraint_(
y.sum() == p
# Instead of:
# model.sum(y[j] for j in FAC) == p
)
# Solve with additional logging output for problem and solution statistics
sol = solve(model, log_output=True)
# Instead of:
# sol = model.solve(log_output=True)
----------------------------- MILP problem statistics -----------------------------
Problem name : p-median
Objective sense : Minimize
Variables : 36 [Binary: 36]
Objective nonzeros : 30
Linear constraints : 36 [Less: 30, Equal: 6]
Nonzeros : 96
RHS nonzeros : 6
Variables : Min LB: 0.000000 Max UB: 1.000000
Objective nonzeros : Min : 9.370000 Max : 123.2700
Linear constraints :
Nonzeros : Min : 1.000000 Max : 1.000000
RHS nonzeros : Min : 1.000000 Max : 2.000000
------------------------------- CPLEX optimizer log -------------------------------
Version identifier: 22.1.2.0 | 2026-03-02 | af0ce9b93
CPXPARAM_Read_DataCheck 1
Found incumbent of value 245.220000 after 0.00 sec. (0.00 ticks)
Tried aggregator 1 time.
Reduced MIP has 36 rows, 36 columns, and 96 nonzeros.
Reduced MIP has 36 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.00 sec. (0.06 ticks)
Probing time = 0.00 sec. (0.02 ticks)
Tried aggregator 1 time.
Detecting symmetries...
Reduced MIP has 36 rows, 36 columns, and 96 nonzeros.
Reduced MIP has 36 binaries, 0 generals, 0 SOSs, and 0 indicators.
Presolve time = 0.00 sec. (0.08 ticks)
Probing time = 0.00 sec. (0.02 ticks)
Clique table members: 35.
MIP emphasis: balance optimality and feasibility.
MIP search method: dynamic search.
Parallel mode: deterministic, using up to 2 threads.
Root relaxation solution time = 0.00 sec. (0.06 ticks)
Nodes Cuts/
Node Left Objective IInf Best Integer Best Bound ItCnt Gap
* 0+ 0 245.2200 0.0000 100.00%
* 0+ 0 163.8300 0.0000 100.00%
0 0 cutoff 163.8300 163.8300 13 0.00%
0 0 cutoff 163.8300 163.8300 13 0.00%
Elapsed time = 0.00 sec. (0.33 ticks, tree = 0.01 MB, solutions = 2)
Root node processing (before b&c):
Real time = 0.00 sec. (0.34 ticks)
Parallel b&c, 2 threads:
Real time = 0.00 sec. (0.00 ticks)
Sync time (average) = 0.00 sec.
Wait time (average) = 0.00 sec.
------------
Total (root+branch&cut) = 0.00 sec. (0.34 ticks)
--------------------------- Solution quality statistics ---------------------------
Incumbent solution:
MILP objective 1.6383000000e+02
MILP solution norm |x| (Total, Max) 7.00000e+00 1.00000e+00
MILP solution error (Ax=b) (Total, Max) 0.00000e+00 0.00000e+00
MILP x bound error (Total, Max) 0.00000e+00 0.00000e+00
MILP x integrality error (Total, Max) 0.00000e+00 0.00000e+00
MILP slack bound error (Total, Max) 0.00000e+00 0.00000e+00
-------------------------------------------------------------------------------------
sol.display(print_zeros=False)
solution for: p-median
objective: 163.830
status: OPTIMAL_SOLUTION(2)
x_1_5 = 1
x_2_5 = 1
x_3_6 = 1
x_4_5 = 1
x_5_5 = 1
y_5 = 1
y_6 = 1
Implement with gurobipy#
from gurobipy import GRB, Model
from opti_extensions.gurobipy import addVars
# Instantiate model
model = Model(name='p-median')
# Add variables
x = addVars(model, indexset=CUST_x_FAC, vtype=GRB.BINARY, name='x')
# Instead of:
# x = model.addVars(CUST_x_FAC, vtype=GRB.BINARY, name='x')
y = addVars(model, indexset=FAC, vtype=GRB.BINARY, name='y')
# Instead of:
# y = model.addVars(FAC, vtype=GRB.BINARY, name='y')
# Set objective
model.setObjective(
cost @ x,
# Instead of:
# quicksum(cost[i, j] * x[i, j] for i in CUST for j in FAC)
sense=GRB.MINIMIZE,
)
# Add constraints
model.addConstrs(
x.sum(i, '*') == 1
# Instead of:
# quicksum(x[i, j] for j in FAC) == 1
for i in CUST
)
model.addConstrs(
x[i, j] <= y[j]
for i, j in CUST_x_FAC
)
model.addConstr(
y.sum() == p
# Instead of:
# quicksum(y[j] for j in FAC) == p
)
# Solve
model.optimize()
Gurobi Optimizer version 13.0.1 build v13.0.1rc0 (linux64 - "Ubuntu 22.04.3 LTS")
CPU model: AMD EPYC 7R13 Processor, instruction set [SSE2|AVX|AVX2]
Thread count: 1 physical cores, 2 logical processors, using up to 2 threads
Optimize a model with 36 rows, 36 columns and 96 nonzeros (Min)
Model fingerprint: 0x1488130c
Model has 30 linear objective coefficients
Variable types: 0 continuous, 36 integer (36 binary)
Coefficient statistics:
Matrix range [1e+00, 1e+00]
Objective range [9e+00, 1e+02]
Bounds range [1e+00, 1e+00]
RHS range [1e+00, 2e+00]
Found heuristic solution: objective 270.1900000
Presolve time: 0.00s
Presolved: 36 rows, 36 columns, 96 nonzeros
Variable types: 0 continuous, 36 integer (36 binary)
Root relaxation: objective 1.638300e+02, 21 iterations, 0.00 seconds (0.00 work units)
Nodes | Current Node | Objective Bounds | Work
Expl Unexpl | Obj Depth IntInf | Incumbent BestBd Gap | It/Node Time
* 0 0 0 163.8300000 163.83000 0.00% - 0s
Explored 1 nodes (21 simplex iterations) in 0.00 seconds (0.00 work units)
Thread count was 2 (of 2 available processors)
Solution count 2: 163.83 270.19
Optimal solution found (tolerance 1.00e-04)
Best objective 1.638300000000e+02, best bound 1.638300000000e+02, gap 0.0000%
model.printAttr('X')
Variable X
-------------------------
x[1,5] 1
x[2,5] 1
x[3,6] 1
x[4,5] 1
x[5,5] 1
y[5] 1
y[6] 1
Implement with Xpress#
import xpress as xp
from opti_extensions.xpress import addVariables
# Instantiate problem
prob = xp.problem(name='p-median')
# Add variables
x = addVariables(prob, indexset=CUST_x_FAC, vartype=xp.binary, name='x')
# Instead of:
# x = prob.addVariables(CUST_x_FAC, vartype=xp.binary, name='x')
y = addVariables(prob, indexset=FAC, vartype=xp.binary, name='y')
# Instead of:
# y = prob.addVariables(FAC, vartype=xp.binary, name='y')
# Set objective
prob.setObjective(
cost @ x,
# Instead of:
# xp.Sum(cost[i, j] * x[i, j] for i in CUST for j in FAC)
sense=xp.minimize,
)
# Add constraints
prob.addConstraint(
x.sum(i, '*') == 1
# Instead of:
# xp.Sum(x[i, j] for j in FAC) == 1
for i in CUST
)
prob.addConstraint(
x[i, j] <= y[j]
for i, j in CUST_x_FAC
)
prob.addConstraint(
y.sum() == p
# Instead of:
# xp.Sum(y[j] for j in FAC) == p
)
# Solve
prob.optimize()
FICO Xpress v9.8.1, Community, solve started 19:57:17, Apr 11, 2026
Heap usage: 455KB (peak 455KB, 84KB system)
Minimizing MILP p-median using up to 2 threads and up to 7734MB memory, with these control settings:
OUTPUTLOG = 1
NLPPOSTSOLVE = 1
XSLP_DELETIONCONTROL = 0
XSLP_OBJSENSE = 1
Original problem has:
36 rows 36 cols 96 elements 36 entities
Presolved problem has:
36 rows 36 cols 96 elements 36 entities
Presolve finished in 0 seconds
Heap usage: 1599KB (peak 1599KB, 84KB system)
Coefficient range original solved
Coefficients [min,max] : [ 1.00e+00, 1.00e+00] / [ 1.00e+00, 1.00e+00]
RHS and bounds [min,max] : [ 1.00e+00, 2.00e+00] / [ 1.00e+00, 2.00e+00]
Objective [min,max] : [ 9.37e+00, 1.23e+02] / [ 9.37e+00, 1.23e+02]
Autoscaling applied standard scaling
Will try to keep branch and bound tree memory usage below 6.7GB
*** Solution found: 245.220000 Time: 0.00 Heuristic: e ***
*** Solution found: 163.830000 Time: 0.00 Heuristic: k ***
Starting concurrent solve with dual (1 thread)
Concurrent-Solve, 0s
Dual
objective dual inf
D 163.83000 .0000000
-------- cutoff --------
Concurrent statistics:
Dual: 16 simplex iterations, 0.00s
Problem is cut off
*** Search completed ***
Uncrunching matrix
Final MIP objective : 1.638300000000000e+02
Final MIP bound : 1.638300000000000e+02
Solution time / primaldual integral : 0.01s/ 100.000000%
Work / work units per second : 0.00 / 0.05
Number of solutions found / nodes : 2 / 0
Max primal violation (abs/rel) : 0.0 / 0.0
Max integer violation (abs ) : 0.0
(<SolveStatus.COMPLETED: 3>, <SolStatus.OPTIMAL: 1>)
print(f'var: {"value" :>20}')
print('-' * 25)
for vd in (x, y):
for k, v in prob.getSolution(vd).items():
if abs(v) > 1e-6:
name = f'{vd.value_name}[{k}]:'
print(f'{name:<15} {v :>8.4f}')
var: value
-------------------------
x[(1, 5)]: 1.0000
x[(2, 5)]: 1.0000
x[(3, 6)]: 1.0000
x[(4, 5)]: 1.0000
x[(5, 5)]: 1.0000
y[5]: 1.0000
y[6]: 1.0000
Implement with highspy#
from highspy import Highs, HighsVarType
from opti_extensions.highspy import addVariables as addVariablesHighs
# Instantiate model
model = Highs()
model.silent()
# Add variables
x = addVariablesHighs(
model,
indexset=CUST_x_FAC,
type=HighsVarType.kInteger,
lb=0,
ub=1,
name_prefix='x',
)
# Instead of:
# x = model.addVariables(
# CUST_x_FAC, type=HighsVarType.kInteger, lb=0, ub=1, name_prefix='x'
#)
y = addVariablesHighs(
model,
indexset=FAC,
type=HighsVarType.kInteger,
lb=0,
ub=1,
name_prefix='y',
)
# Instead of:
# y = model.addVariables(FAC, type=HighsVarType.kInteger, lb=0, ub=1, name_prefix='y')
# Set objective
model.minimize(
cost @ x
# Instead of:
# Highs.qsum(cost[i, j] * x[i, j] for i in CUST for j in FAC)
)
# Add constraints
model.addConstrs(
x.sum(i, '*') == 1
# Instead of:
# Highs.qsum(x[i, j] for j in FAC) == 1
for i in CUST
)
model.addConstrs(
x[i, j] <= y[j]
for i, j in CUST_x_FAC
)
model.addConstr(
y.sum() == p
# Instead of:
# Highs.qsum(y[j] for j in FAC) == p
)
# Solve
model.solve()
<HighsStatus.kOk: 0>
print(f'{"var:":<15} {"value":>8}')
print('-' * 25)
for vd in (x, y):
for idx, var in vd.items():
val = model.val(var)
if abs(val) > 1e-6:
name = f'{vd.value_name}[{idx}]:'
print(f'{name:<15} {val:>8.4f}')
var: value
-------------------------
x[(1, 5)]: 1.0000
x[(2, 5)]: 1.0000
x[(3, 6)]: 1.0000
x[(4, 5)]: 1.0000
x[(5, 5)]: 1.0000
y[5]: 1.0000
y[6]: 1.0000
Total running time of the script: (0 minutes 0.031 seconds)