2
$\begingroup$

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])

def mclp(points,K,radius,M):
    print('----- Configurations -----')
    print('  Number of points %g' % points.shape[0])
    print('  K %g' % K)
    print('  Radius %g' % radius)
    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)
    mask1 = D<=radius
    D[mask1]=1
    D[~mask1]=0
    # Build model
    m = Model()
    # Add variables
    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()
    # Add constraints
    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)

def plot_result(points,opt_sites,radius):
    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:
        circle = plt.Circle(site, radius, color='C1',fill=False,lw=2)
        ax.add_artist(circle)
    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
radius = 0.2
M = 100
opt_sites,f = mclp(points,K,radius,M)
plot_result(points,opt_sites,radius)

Mathematical Formulation

$\endgroup$
4
  • 1
    $\begingroup$ 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. $\endgroup$
    – joni
    Jun 18 at 5:48
  • $\begingroup$ 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. $\endgroup$ Jun 18 at 14:21
  • $\begingroup$ @joni Just added the mathematical model I am using. $\endgroup$
    – izc2300
    Jun 18 at 14:59
  • $\begingroup$ @PrameshKumar Yes! I just added clarification. $\endgroup$
    – izc2300
    Jun 18 at 15:00
2
$\begingroup$

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)
$\endgroup$
3
  • $\begingroup$ 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? $\endgroup$
    – A.Omidi
    Jun 21 at 9:07
  • $\begingroup$ Yes, this is possible but @izc2300 seems to accept that solution when I tried to ran their previous script. $\endgroup$ Jun 25 at 1:25
  • $\begingroup$ 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. 👍 $\endgroup$
    – A.Omidi
    Jun 25 at 10:53

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.