Code with Sample data:\
import gurobipy as gp
from gurobipy import *
import numpy as np
from scipy.spatial.distance import pdist
from collections import namedtuple
from scipy.spatial import distance_matrix as dm
N =15 #Sites" Towers
R = 86
#Generate Data
Points = namedtuple('Points', 'x, y')
def gen_data(Points,N):
locations = []
for p in range(N):
x,y = np.random.uniform(1,10,size = 2)
points = Points(x,y)
if points.x**2 + points.y**2 <= R**2:
locations.append(points)
return np.array([[p.x,p.y] for p in locations],dtype='float16')
#return sites
#print(sites)
#Generate sites
sites = gen_data(Points,N)
towers = gen_data(Points,7)
#print(towers)
N,_ = sites.shape
I = [*range(N)]
T,_ = towers.shape
J = [*range(T)]
K = 5
maxpoints = 13 #per tower circle
minpoints = 2 #per tower circle
#Generate distance D(i,j)
D = dm(sites,towers).astype('float16')
print(D)
#Masking
mask = (D < 7)
D[mask] = 1; D[~mask]=0
model = Model('fac')
x = model.addVars(J,vtype='b',name='x')
y = model.addVars(I,vtype='b',name='y')
C1 = model.addConstr(x.sum() == K,'C1')
C2 = model.addConstrs((minpoints*x[j] <= quicksum(y[i] for i in np.where(D[:,j]==1)[0])\
for j in J),'C2')
#C3 = model.addConstrs((quicksum(x[j] for j in np.where(D[i,:]==1)[0])<= K*(1-y[i])+y[i]\
# for i in I),'NonOverlapping')
for j in J),'C2')
C5 = model.addConstrs((maxpoints >= quicksum(y[i] for i in np.where(D[:,j]==1)[0])\
for j in J),'C2')
C4 = model.addConstrs((y[i] == 1 for i in np.where(np.sum(D,axis=1)>=1)[0]),'C4')
obj = y.sum()
model.setObjective(obj,GRB.MAXIMIZE)
model.update()
model.optimize()
model.printAttr('x')