The Facility Dispersion Problem: Mixed-Integer Programming Models
A comprehensive Python tutorial on the p-dispersion and maxisum models
In some facility location problems, facilities need to be positioned so that the influence of one doesn’t overshadow or adversely affect the others. Whether driven by motives of risk mitigation, avoidance of competition, or other considerations, mastering the solutions to these challenges is a crucial skill in the operations research toolkit.
Kuby (1987) proposed two distinct mixed-integer programming (MIP) formulations for these problems: the p-dispersion problem and the maxisum problem. In this article, both will be implemented using the Python library Pyomo and the HiGHS solver.
Besides the two models, some useful modeling resources will be presented. First, a strategy for linearizing the product of binary variables in the context of MIP, although considering the maximization objectives, doesn’t need to be explicit in this problem. Second, a max-min MIP formulation, which intends to maximize something smaller than any parameter of a set of items if the item is selected. Finally, a strategy for solving multiple objectives with a well-defined hierarchy of priorities that combines elements of the two models.
Those interested not yet familiar with numerical optimization might find it helpful to take a look at my previous stories about Linear Programming and the Branch & Bound method before proceeding with this one.
As usual, you can find the complete code in this git repository.
Which of these locations would you choose to place 5 facilities?
The product of binary variables
When defining the basic elements of this problem, one can use binary variables that assume the value of 1 if a location is selected and 0 otherwise. Let us denote these variables xᵢ. Suppose the distances (or some other dispersion metric) between two locations are computed in advance and let us denote these as dᵢⱼ. How could we calculate the dispersion of any pair of facilities selected?
It is somehow intuitive in such a situation to formulate it using the product of the binary variables xᵢ and xⱼ to compute the dissimilarity obtained when they are both included in the solution. This approach is equivalent to a logical AND statement. However, this would lead to a quadratic formulation, therefore one wouldn’t be able to apply linear solvers. Fortunately, there is a way to formulate the product of binary variables in MIP using a few constraints.
Consider a directed graph G(V, A) and zᵢⱼ a new binary variable that indicates both nodes i and j are selected. If either i or j is not selected, zᵢⱼ must be 0. This leads to the first group of constraints:
So far, even if both i and j are selected, zᵢⱼ can assume the value of 0. Therefore we must include an additional constraint such that when i and j are selected zᵢⱼ becomes 1.
When maximizing something proportional to zᵢⱼ, as in the maxisum problem, the second constraint group presented is unnecessary, as it would make no sense not to compute a gain proportional to zᵢⱼ if it is feasible. However, this can be useful when formulating other MIP models.
In the following section, let us apply these concepts to the maxisum problem.
The maxisum model
The discrete maxisum problem must define the location of p facilities among N discrete nodes to maximize the sum of distances (or the average of distances) computed between all pairs of selected facilities. Therefore, consider the facilities are nodes displaced in a directed graph G(V, A). The weight of each arc from i to j is the distance (dispersion) metric dᵢⱼ known in advance. Also, consider binary variables xᵢ to indicate if a location i is selected and zᵢⱼ to indicate if both i and j are selected.
The objective can be stated as the following:
Besides the constraints to linearize the product of binary variables presented in the previous section, it is necessary to include a constraint to guarantee that p locations are selected.
Thus, we can state the constraints of the problem as:
Let us put that into code using Python. To do so, we must start by importing the libraries which will be used. The library numpy will be used for linear algebra calculations and to create random coordinate points; the functions squareform and pdist from scipy will be used to compute distances given a matrix of coordinates; matplotlib will be our visualization tool; and pyomo the algebraic modeling language (AML) (with solver HiGHS).
import numpy as np
from scipy.spatial.distance import squareform, pdist
import matplotlib.pyplot as plt
import pyomo.environ as pyo
from pyomo.contrib.appsi.solvers.highs import HighsWe must create N points and define how many of those must be selected as facility locations. By fixing the random seed (12) in every code execution the points represented in the introduction should be obtained.
# Fix random seeds
np.random.seed(12)
# Create random points
N = 25
p = 5
coordinates = np.random.random((N, 2))
# Calculate pairwise distances
distances = squareform(pdist(coordinates))We now have the necessary elements to start our pyomo model.
There are two approaches for modeling a problem in pyomo: Abstract and Concrete models. In the first approach, the algebraic expressions of the problem are defined before some data values are supplied, whereas, in the second, the model instance is created immediately as its elements are defined. You can find more about these approaches in the library documentation or in the book by Bynum et al. (2021). Throughout this article, we will adopt the Concrete model formulation.
model = pyo.ConcreteModel()In this model, we have two sets: nodes and arcs. As we are considering a complete directed graph, there exist arcs between every pair of two nodes except from the node to itself.
# Sets of nodes and arcs
model.V = pyo.Set(initialize=range(N))
model.A = pyo.Set(initialize=[(i, j) for i in model.V for j in model.V if i != j])The parameters provided beforehand are the number of nodes that must be selected and the distances between pairs of nodes.
# Parameters
model.d = pyo.Param(model.A, initialize={(i, j): distances[i, j] for (i, j) in model.A})
model.p = pyo.Param(initialize=p)Then, let us include the decision variables.
# Decision variables
model.x = pyo.Var(model.V, within=pyo.Binary)
model.z = pyo.Var(model.A, within=pyo.Binary)And the constraints…
# p nodes are selected
def p_selection(model):
return sum(model.x[:]) == model.p
# If starting node is not selected, the arc is 0
def dispersion_c1(model, i, j):
return model.z[i, j] <= model.x[i]
# If ending node is not selected, the arc is 0
def dispersion_c2(model, i, j):
return model.z[i, j] <= model.x[j]
# Include constraints in model
model.p_selection = pyo.Constraint(rule=p_selection)
model.dispersion_c1 = pyo.Constraint(model.A, rule=dispersion_c1)
model.dispersion_c2 = pyo.Constraint(model.A, rule=dispersion_c2)Finally, we must create the objective function.
def disp_obj(model):
return sum(model.z[i, j] * model.d[i, j] for (i, j) in model.A)
model.obj = pyo.Objective(rule=disp_obj, sense=pyo.maximize)Now the maxisum model is ready to be solved. We must then instantiate a solver compatible with pyomo to handle it. The Highs solver is available in Pyomo (check the imports) since version 6.4.3 if the highspy package is installed. So make sure to run a pip install highspy .
solver = Highs()
solver.solve(model)And, after approximately 120s of computing time, we have the following results:
Notice that, even though the total dispersion is maximized, two facilities in the upper left corner are still quite close to each other, which can be undesirable. Thus, alternatively to the maxisum formulation we have the p-dispersion model in which we maximize the minimum distance between any pair of selected nodes.
The p-dispersion model
In the p-dispersion model, we need an additional decision variable to compute the minimum distance between all pairs of selected nodes, which is our goal to maximize. Let us denote this variable D. We must create a big M constraint that ensures that if both i and j are selected, D is lesser than or equal to dᵢⱼ; otherwise, we must ensure that D is unlimited. If we keep the formulation with zᵢⱼ as the product of binary variables, we can formulate this additional constraint as the following.
In this formulation, M is a fixed parameter arbitrarily large used to formulate a disjunctive rule. A good choice of M should be large enough to guarantee the feasibility of the constraint (i, j) for any value of D when zᵢⱼ is zero, but as small as possible so that the linear relaxation of the model is similar to the integer version, which makes convergence easier. In this problem, the largest value of dᵢⱼ can be an interesting alternative.
Although this formulation would work well for this model, a more concise approach can be used, in which variables zᵢⱼ are not even included in the model.
In this formulation, either xᵢ or xⱼ equal to zero is a sufficient condition to guarantee that the inequality is valid for any value of D. And the objective becomes simply maximizing D.
In Python code next, consider we have a new model with the same sets and parameters of the previous one as well as the group of decision variables x.
# Max-min constraint
def maxmin_rule(model, i, j):
return model.D <= model.d[i, j] + model.M * (1 - model.x[i]) + model.M * (1 - model.x[j])
# New parameter big M
model.M = max(model.d[:, :])
# New variable
model.D = pyo.Var(within=pyo.NonNegativeReals)
# New constraint
model.maxmin_rule = pyo.Constraint(model.A, rule=maxmin_rule)
# Objective
model.obj = pyo.Objective(expr=model.D, sense=pyo.maximize)And after calling the solver it took less than 1.2s to obtain the following results.
Which seems great as the locations are well distributed in space.
Is there a way to improve this distribution?
A multicriteria approach
Remember the objective function of the p-dispersion model depends only on the minimum distance between any pair of nodes selected. Therefore several solutions can be obtained using combinations of the two points defining this distance and others with distances to each other greater than or equal to the objective itself. Can we refine our solution by selecting the best of these alternatives? That leads to the two-step “multicriteria approach” proposed by Kuby (1987).
In the first step, the p-dispersion model is solved to optimality, and the current objective is stored in a parameter d_opt. Then, a maxisum model with an additional constraint ensuring that D ≥ d_opt is solved to obtain, among the optimal solutions to the p-dispersion model, the one with the maximum total dispersion.
To do that in Python, consider you have a p-dispersion model instantiated with also decision variables and constraints of the maxisum model.
# D must be optimal
def composed_constr(model):
return model.D >= model.d_opt
# Solve p-dispersion
solver.solve(model)
# New parameter
model.d_opt = pyo.Param(initialize=model.obj())
# Deactivate old objective
model.obj.deactivate()
# Solution will not make the current D worse
model.composed_constr = pyo.Constraint(rule=composed_constr)
# New objective
model.obj_disp = pyo.Objective(rule=disp_obj, sense=pyo.maximize)
solver.solve(model)And this is surprisingly fast because when the solver enters the second objective the space of feasible alternatives is substantially small. In less than one (additional) second, it leads to the following results.
Further reading
When customers aren’t uniformly distributed, facilities have a limited capacity, or perhaps the adequate number of facilities is unknown beforehand, probably you are facing a different Facility Location Problem. You can find the formulation of Balinski (1965) implemented in Python using PuLP in this amazing story by Nicolo Cosimo Albanese.
Conclusions
In this article, two mixed-integer programming models for the facility dispersion problem were implemented in Python using Pyomo. As previously verified in the literature, the maxisum model can lead to an uneven distribution of elements in space whereas the p-dispersion model produced solutions with locations well spread in space and evenly distributed. A maxisum objective can be used to refine solutions of the p-dispersion model by selecting the one with the greatest total dispersion from the set of optimal solutions. The complete code used in these examples is available for further use.
References
Balinski, M. L. 1965. Integer programming: methods, uses, computations. Management science, 12(3), 253–313.
Bynum, M. L., Hackebeil, G. A., Hart, W. E., Laird, C. D., Nicholson, B. L., Siirola, J. D., … & Woodruff, D. L. 2021. Pyomo-optimization modeling in python (Vol. 67). Berlin/Heidelberg, Germany: Springer.
Kuby, M. J. 1987. Programming models for facility dispersion: The p‐dispersion and maxisum dispersion problems. Geographical Analysis, 19(4), 315–329.
