The discrete maxisum problem must define the placement of p facilities amongst N discrete nodes to maximise the sum of distances (or the common of distances) computed between all pairs of chosen facilities. Subsequently, consider the facilities are nodes displaced in a directed graph G(V, A). The burden of every arc from i to j is the gap (dispersion) metric dᵢⱼ known prematurely. Also, consider binary variables xᵢ to point if a location i is chosen and zᵢⱼ to point if each i and j are chosen.

The target may be stated as the next:

Besides the constraints to linearize the product of binary variables presented within the previous section, it’s mandatory to incorporate a constraint to ensure that p locations are chosen.

Thus, we will state the constraints of the issue as:

Allow us to put that into code using Python. To accomplish that, we must start by importing the libraries which shall be used. The library numpy shall be used for linear algebra calculations and to create random coordinate points; the functions squareform and pdist from scipy shall be used to compute distances given a matrix of coordinates; matplotlib shall 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 Highs

We must create N points and define how a lot of those have to be chosen as facility locations. By fixing the random seed (12) in every code execution the points represented within the introduction ought to 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 mandatory elements to start out our pyomo model.

There are two approaches for modeling an issue in pyomo: Abstract and Concrete models. In the primary approach, the algebraic expressions of the issue are defined before some data values are supplied, whereas, within the second, the model instance is created immediately as its elements are defined. You could find more about these approaches within the library documentation or within the book by Bynum et al. (2021). Throughout this text, we are going to adopt the Concrete model formulation.

model = pyo.ConcreteModel()

On this model, we’ve got two sets: nodes and arcs. As we’re considering a whole 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 variety of nodes that have to be chosen 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, allow us to include the choice variables.

# Decision variables
model.x = pyo.Var(model.V, inside=pyo.Binary)
model.z = pyo.Var(model.A, inside=pyo.Binary)

And the constraints…

# p nodes are chosen
def p_selection(model):
return sum(model.x[:]) == model.p
# If starting node isn't chosen, the arc is 0
def dispersion_c1(model, i, j):
return model.z[i, j] <= model.x[i]
# If ending node isn't chosen, 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 target 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 able to be solved. We must then instantiate a solver compatible with pyomo to handle it. The Highs solver is obtainable in Pyomo (check the imports) since version 6.4.3 if the highspy package is installed. So be certain to run a pip install highspy .

solver = Highs()
solver.solve(model)

And, after roughly 120s of computing time, we’ve got the next results:

Notice that, although the entire dispersion is maximized, two facilities within the upper left corner are still quite close to one another, which may be undesirable. Thus, alternatively to the maxisum formulation we’ve got the p–dispersion model by which we maximize the minimum distance between any pair of chosen nodes.