3.6 Facility location problem#
This notebook illustrates Example 20 from Chapter 3 of the textbook Hands-On Mathematical Optimization with Python, Cambridge University Press, 2024.
For details about these models and more information we refer to the book.
Also notice that at the end of Chapter 3 you may find a number of exercises that show that versions of this problem aiming at maximizing the importance of the customers served, when facing constraints that limit the possibility of serving all, leads to even more striking evidence for the importance of selecting the right model: it can offer large factors of speedup in solving the relevant instances!
Preamble: Install Pyomo and solvers#
The following cell sets and verifies a Pyomo installation and the solver CBC. If run on Google Colab, the cell installs Pyomo and the CBC solver, while, if run elsewhere, it assumes they all have been previously installed.
import sys
if 'google.colab' in sys.modules:
%pip install pyomo >/dev/null 2>/dev/null
!apt-get install -y -qq coinor-cbc
Problem description#
Consider the problem of a supplier facing the task of fulfilling specific customer demands with minimal costs while simultaneously deciding how many facilities to build and where. In terms of data, we are given a set \(I\) of customers, a set \(J\) of possible locations, the cost \(c_j\) of building facility \(j\), and the cost \(h_{ij}\) incurred to satisfy the demands of customer \(i\) at facility \(j\).
In this notebook we will formulate two equivalent models for this problem and compare their performance.
First MILO formulation#
Introduce two sets of binary variables,
and
The resulting MILO is
Second MILO formulation#
Let \(n=|J|\) be the number of possible facility locations and \(m=|I|\) the number of customers. Note that, since \(\sum_j y_{ij} =1\) for every \(i \in I\), we can replace the \(n \times m\) constraints \(y_{ij} \leq x_j\) by only \(n\) constraints, namely
The second MILO is
This approach leads to a more concise mathematical formulation of the model. However, it may not be a good idea if we want to solve the problem using its linear relaxation. Indeed, by reducing the number of constraints, we inadvertently made the feasible region of the relaxation larger and less tight around the feasible integer points. This fact becomes clearly evident in the increased run-time required to solve the optimization problem when working with these weaker constraints, see figures below.
Instead of defining two separate Pyomo models, we will define a single common model and then add either the strong or weaker formulation of the “facility built before use” constraint.
import pyomo.environ as pyo
import numpy as np
import pandas as pd
import itertools as it
import matplotlib.pyplot as plt
from time import perf_counter as pc
from tqdm.notebook import tqdm
import subprocess
Solvers and options#
def ListAvailableSolvers():
shell_command = "pyomo help --solvers"
output = subprocess.check_output(shell_command, shell=True).decode()
return [
line.strip()[1:]
for line in output.split()
if line.strip().startswith("+") and not line.strip().endswith(")")
]
def GetSolverName(solver):
try:
return solver.name
except:
return "appsi_highs"
def SwitchCutsOff(solver):
solver_name = GetSolverName(solver)
if "cbc" in solver_name:
solver.options["Cuts"] = "off"
elif "cplex" in solver_name:
solver.options["mip_cuts_bqp"] = -1
solver.options["mip_cuts_cliques"] = -1
solver.options["mip_cuts_covers"] = -1
solver.options["mip_cuts_disjunctive"] = -1
solver.options["mip_cuts_flowcovers"] = -1
solver.options["mip_cuts_pathcut"] = -1
solver.options["mip_cuts_gomory"] = -1
solver.options["mip_cuts_gubcovers"] = -1
solver.options["mip_cuts_implied"] = -1
solver.options["mip_cuts_localimplied"] = -1
solver.options["mip_cuts_liftproj"] = -1
solver.options["mip_cuts_mircut"] = -1
solver.options["mip_cuts_mcfcut"] = -1
solver.options["mip_cuts_rlt"] = -1
solver.options["mip_cuts_zerohalfcut"] = -1
elif "gurobi" in solver_name:
solver.options["Cuts"] = 0
elif "highs" in solver_name:
pass # there is no way to do this, see
# https://ergo-code.github.io/HiGHS/dev/options/definitions/
elif "xpress" in solver_name:
solver.options["CUTSTRATEGY"] = 0
else:
pass
return solver
def SwitchPresolveOff(solver):
solver_name = GetSolverName(solver)
if "cbc" in solver_name:
solver.options["Presolve"] = "off"
elif "cplex" in solver_name:
solver.options["preprocessing_presolve"] = 0
elif "gurobi" in solver_name:
solver.options["Presolve"] = 0
elif "highs" in solver_name:
solver.options["presolve"] = "off"
elif "xpress" in solver_name:
solver.options["PRESOLVE"] = 0
else:
pass
return solver
def LimitSolveTime(solver, max_in_seconds):
solver_name = GetSolverName(solver)
if "cbc" in solver_name:
solver.options["seconds"] = max_in_seconds
elif "cplex" in solver_name:
solver.options["timelimit"] = max_in_seconds
elif "gurobi" in solver_name:
solver.options["TimeLimit"] = max_in_seconds
elif "highs" in solver_name:
solver.options["time_limit"] = max_in_seconds
elif "xpress" in solver_name:
solver.options["MAXTIME"] = max_in_seconds
else:
pass
return solver
def ClosureForMaxTime(max_in_seconds):
return lambda solver: LimitSolveTime(solver, max_in_seconds)
def compose(*functions):
def composed(*args, **kwargs):
result = args[0]
for func in functions:
result = func(result)
return result
return composed
initial_solvers = ListAvailableSolvers()
initial_solvers
['appsi_gurobi',
'appsi_highs',
'cbc',
'gdpopt',
'gdpopt.gloa',
'gdpopt.lbb',
'gdpopt.loa',
'gdpopt.ric',
'gurobi',
'gurobi_direct',
'gurobi_persistent',
'ipopt',
'mindtpy',
'mindtpy.ecp',
'mindtpy.fp',
'mindtpy.goa',
'mindtpy.oa',
'mosek',
'mosek_direct',
'mosek_persistent',
'mpec_minlp',
'mpec_nlp',
'multistart',
'scipy.fsolve',
'scipy.newton',
'scipy.root',
'scipy.secant-newton',
'trustregion']
The models#
def FacilityLocationCommon(installation, service):
model = pyo.ConcreteModel("Facility location")
model.facilities = pyo.Set(initialize=range(len(installation)))
model.customers = pyo.Set(initialize=range(len(service)))
model.x = pyo.Var(model.facilities, within=pyo.Binary)
model.y = pyo.Var(model.customers, model.facilities, within=pyo.Binary)
@model.Objective(sense=pyo.minimize)
def obj(model):
return sum(installation[j] * model.x[j] for j in model.facilities) + sum(
service[i][j] * model.y[i, j]
for i in model.customers
for j in model.facilities
)
@model.Constraint(model.customers)
def ChooseOneFacility(model, i):
return sum(model.y[i, j] for j in model.facilities) == 1
return model
def FacilityLocationWeak(installation, service):
model = FacilityLocationCommon(installation, service)
model.name += " weak model"
@model.Constraint(model.facilities)
def ServeIfOpen(model, j):
return (
sum(model.y[i, j] for i in model.customers)
<= len(model.customers) * model.x[j]
)
return model
def FacilityLocationStrong(installation, service):
model = FacilityLocationCommon(installation, service)
model.name += " strong model"
@model.Constraint(model.customers, model.facilities)
def ServeIfOpen(model, i, j):
return model.y[i, j] <= model.x[j]
return model
def GetSolution(model):
X = [model.x[j]() > 0.5 for j in model.facilities]
Y = [[model.y[i, j]() > 0.5 for j in model.facilities] for i in model.customers]
return X, Y, model.obj()
We introduce functions to generate random instances for the problem and to then visualize the optimal solution.
def GenerateFacilityLocationInstance(nofFacilities, nofCustumers):
facilities = range(nofFacilities)
customers = range(nofCustumers)
xC = np.random.randint(0, 100, nofCustumers)
yC = np.random.randint(0, 100, nofCustumers)
xF = np.random.randint(0, 100, nofFacilities)
yF = np.random.randint(0, 100, nofFacilities)
installation = np.random.randint(1000, 2000, nofFacilities)
dist = lambda i, j: ((xC[i] - xF[j]) ** 2 + (yC[i] - yF[j]) ** 2)
service = [[dist(i, j) for j in facilities] for i in customers]
return installation, service, xC, yC, xF, yF
def ShowFacilityLocation(xC, yC, xF, yF, X=[], Y=[], value=None, title=None, ax=None):
if ax is None:
ax = plt.gca()
[
ax.plot([xC[i], xF[j]], [yC[i], yF[j]], "g-")
for j in range(len(X))
if X[j]
for i in range(len(Y))
if Y[i][j]
]
ax.plot(xC, yC, "o")
ax.plot(xF, yF, "s")
if title:
ax.set_title(title)
elif value:
ax.set_title(f"Optimal value: {value}")
We then generate an instance of the problem with 10 facilities and 100 customers and initialize both a strong and a weak Pyomo models.
np.random.seed(2023)
installation, service, xC, yC, xF, yF = GenerateFacilityLocationInstance(20, 200)
ShowFacilityLocation(xC, yC, xF, yF)
weak = FacilityLocationWeak(installation, service)
strong = FacilityLocationStrong(installation, service)
def Solve(solver, model):
t = pc()
_ = solver.solve(model)
t = pc() - t
X, Y, v = GetSolution(model)
title = f"{model.name} solved by {GetSolverName(solver)} in {t:.2f} seconds {v:.0f}"
return X, Y, v, title
solver = pyo.SolverFactory("cbc")
ShowFacilityLocation(xC, yC, xF, yF, *Solve(solver, weak))
ShowFacilityLocation(xC, yC, xF, yF, *Solve(solver, strong))
However, this is not a proper performance comparison between the two model formulations, since under the hood cbc is already solving the instance in a clever way. The commercial solver Gurobi and the open-source solver HiGHS also have a similar feature.
If run on Google Colab, the following cell installs three additional commercial solvers; if run elsewhere, it assumes they all have been previously installed.
if 'google.colab' in sys.modules:
%pip install gurobipy >/dev/null 2>/dev/null
%pip install cplex >/dev/null 2>/dev/null
%pip install xpress >/dev/null 2>/dev/null
available_solvers = ListAvailableSolvers()
set(available_solvers) - set(initial_solvers)
set()
installation, service, xC, yC, xF, yF = GenerateFacilityLocationInstance(8, 88)
weak = FacilityLocationWeak(installation, service)
strong = FacilityLocationStrong(installation, service)
solver_names = ["cplex_direct", "gurobi_direct", "xpress_direct", "appsi_highs"]
solver_names = sorted(set(solver_names) & set(ListAvailableSolvers()))
_, axs = plt.subplots(len(solver_names), 2, figsize=(20, 10))
plt.subplots_adjust(hspace=0.5)
for i, solver_name in enumerate(solver_names):
solver = SwitchCutsOff(pyo.SolverFactory(solver_name))
ShowFacilityLocation(xC, yC, xF, yF, *Solve(solver, weak), ax=axs[i, 0])
ShowFacilityLocation(xC, yC, xF, yF, *Solve(solver, strong), ax=axs[i, 1])
Running HiGHS 1.5.3 [date: 2023-05-16, git hash: 594fa5a9d]
Copyright (c) 2023 HiGHS under MIT licence terms
Running HiGHS 1.5.3 [date: 2023-05-16, git hash: 594fa5a9d]
Copyright (c) 2023 HiGHS under MIT licence terms
We now run a more extensive performance comparison considering increasingly larger instances of the problem. We will see that the strong formulation is consistently faster than the weak formulation.
def SolveInstances(dimensions, solvers, seed, modify=lambda s: s):
np.random.seed(seed)
model = {"weak": FacilityLocationWeak, "strong": FacilityLocationStrong}
df = pd.DataFrame(
index=dimensions,
columns=[f"{s}_{m}" for s, m in it.product(solvers, model.keys())],
)
for n, m in tqdm(df.index):
installation, service, *_ = GenerateFacilityLocationInstance(n, m)
instance = {
option: model[option](installation, service) for option in model.keys()
}
values = dict.fromkeys(solvers)
for column in df.columns:
solver = column[: column.rfind("_")]
option = column[column.rfind("_") + 1 :]
t = pc()
modify(pyo.SolverFactory(solver)).solve(instance[option])
df.at[(n, m), column] = pc() - t
*_, v = GetSolution(instance[option])
values[solver] = v
if len(set(values.values())) > 1:
print(values)
return df
We solve instances with 10 to 100 facilities and 100 to 1000 customers. We then plot the run-time of the strong and weak formulations as a function of the number of facilities.
if "google.colab" in sys.modules:
solvers = ["cbc"]
else:
solvers = ["gurobi_direct", "cplex_direct"]
dimensions = [(n, 10 * n) for n in range(10, 51, 5)]
SolveInstances(dimensions, solvers, 2023).plot(figsize=(8, 5))
plt.show()
modify = compose(ClosureForMaxTime(10), SwitchCutsOff)
SolveInstances(dimensions, solvers, 2023, modify).plot(figsize=(8, 5))
plt.show()
WARNING: Loading a SolverResults object with an 'aborted' status, but
containing a solution
WARNING: Loading a SolverResults object with an 'aborted' status, but
containing a solution
WARNING: Loading a SolverResults object with an 'aborted' status, but
containing a solution
WARNING: Loading a SolverResults object with an 'aborted' status, but
containing a solution
WARNING: Loading a SolverResults object with an 'aborted' status, but
containing a solution
WARNING: Loading a SolverResults object with an 'aborted' status, but
containing a solution
SolveInstances(dimensions, solvers, 2023, SwitchPresolveOff).plot(figsize=(8, 5))
plt.show()
