# Capacitated Maximum Coverage Location Problem, Python and Gurobi

I am building a variant of the maximum coverage location model and want to limit the amount of points that each "facility" can cover. I am using Gurobi optimization . I have tried using the AddConstr() function but have failed. Each K represents a facility and r a radius. My goal is to add a constraint that would limit the amount of points inside that radius. points is a 2D array of x-y coordinates. K is the inputted number of facilities, and M is the total number of points. Would I be able to accomplish this inside this function or would I need to add another variable?

import numpy as np
from scipy.spatial import distance_matrix
from gurobipy import *
from scipy.spatial import ConvexHull
from shapely.geometry import Polygon, Point
from numpy import random

def generate_candidate_sites(points,M=100):
hull = ConvexHull(points)
polygon_points = points[hull.vertices]
poly = Polygon(polygon_points)
min_x, min_y, max_x, max_y = poly.bounds
sites = []
while len(sites) < M:
random_point = Point([random.uniform(min_x, max_x),
random.uniform(min_y, max_y)])
if (random_point.within(poly)):
sites.append(random_point)
return np.array([(p.x,p.y) for p in sites])

print('----- Configurations -----')
print('  Number of points %g' % points.shape[0])
print('  K %g' % K)
print('  M %g' % M)
import time
start = time.time()
sites = generate_candidate_sites(points,M)
J = sites.shape[0]
I = points.shape[0]
D = distance_matrix(points,sites)
# Build model
m = Model()
x = {}
y = {}
for i in range(I):
y[i] = m.addVar(vtype=GRB.BINARY, name="y%d" % i)
for j in range(J):
x[j] = m.addVar(vtype=GRB.BINARY, name="x%d" % j)

m.update()
m.addConstr(quicksum(x[j] for j in range(J)) == K)

for i in range(I):
m.addConstr(quicksum(x[j] for j in np.where(D[i]==1)[0]) >= y[i])

m.setObjective(quicksum(y[i]for i in range(I)),GRB.MAXIMIZE)
m.setParam('OutputFlag', 0)
m.optimize()
end = time.time()
print('----- Output -----')
print('  Running time : %s seconds' % float(end-start))
print('  Optimal coverage points: %g' % m.objVal)

solution = []
if m.status == GRB.Status.OPTIMAL:
for v in m.getVars():
# print v.varName,v.x
if v.x==1 and v.varName[0]=="x":
solution.append(int(v.varName[1:]))
opt_sites = sites[solution]
return opt_sites,m.objVal

def plot_input(points):
from matplotlib import pyplot as plt
fig = plt.figure(figsize=(8,8))
plt.scatter(points[:,0],points[:,1],c='C0')
ax = plt.gca()
ax.axis('equal')
ax.tick_params(axis='both',left=False, top=False, right=False,
bottom=False, labelleft=False, labeltop=False,
labelright=False, labelbottom=False)

from matplotlib import pyplot as plt
fig = plt.figure(figsize=(8,8))
plt.scatter(points[:,0],points[:,1],c='C0')
ax = plt.gca()
plt.scatter(opt_sites[:,0],opt_sites[:,1],c='C1',marker='+')
for site in opt_sites:
ax.axis('equal')
ax.tick_params(axis='both',left=False, top=False, right=False,
bottom=False, labelleft=False, labeltop=False,
labelright=False, labelbottom=False)

Code to run a small sample program

import numpy as np
Npoints = 300
from sklearn.datasets import make_moons
points,_ = make_moons(Npoints,noise=0.15)
K = 20
M = 100

• This is rather a modelling question than a programming question, so providing the mathematical model (written in LaTeX) instead of the code makes it easier to help. IMHO, the code is only necessary if you already know how to formulate a constraint and struggle to implement it.
– joni
Jun 18, 2021 at 5:48
• Could you please clarify the model? I am guessing x[j] represents whether or not to locate facility. What does y[j] represent? I think you might need a variable z[i, j] that assigns point i to facility j. Jun 18, 2021 at 14:21
• @joni Just added the mathematical model I am using. Jun 18, 2021 at 14:59
• @PrameshKumar Yes! I just added clarification. Jun 18, 2021 at 15:00

Similar to the last constraint, you can define a neighborhood of points around a facility, and specify the following constraint:

$$\sum_{i \in N(j)} y_{i} \le \max{\rm Points}, \forall j \in J$$

where $$N(j) = \{i\in I: d_{i, j} \le r\}$$.

The following code snippet should work:

for j in range(J):
m.addConstr(quicksum(y[i] for i in np.where(D[:,j]==1)[0]) <= 10)
• based on the constraint that you mentioned it might assign the same demand point to the different facility that causes some strange results. For example, demand point one is being assigned to facilities one and two simultaneously. have you faced that? Jun 21, 2021 at 9:07
• Yes, this is possible but @izc2300 seems to accept that solution when I tried to ran their previous script. Jun 25, 2021 at 1:25
• that's right. In the real situation it causes something like the symmetry. It would be better if you could update the answer to cover that. 👍 Jun 25, 2021 at 10:53