read_solve_mclp_gurobi.py

# Adapted from https://github.com/cyang-kth/maximum-coverage-location

import numpy as np
from gurobipy import *
from numpy import random
from scipy.spatial import distance_matrix
import argparse

def mclp(points, sites, demands, B, radius):
    print('----- Configurations -----')
    print('  Number of points %g' % points.shape[0])
    print('  Number of sites %g' % sites.shape[0])
    print('  B %g' % B)
    print('  Radius %g' % radius)
    import time
    start = time.time()
    I = sites.shape[0]
    J = points.shape[0]
    D = distance_matrix(points, sites)
    mask1 = D<=radius
    D[mask1]=1
    D[~mask1]=0
    m = Model()
    y = {}
    z = {}
    for i in range(I):
      y[i] = m.addVar(vtype=GRB.BINARY, name="y%d" % i)
    for j in range(J):
      z[j] = m.addVar(vtype=GRB.BINARY, name="x%d" % j)

    m.update()
    m.addConstr(quicksum(y[i] for i in range(I)) <= B)

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

    m.setObjective(quicksum(z[j]*demands[j] for j in range(J)),GRB.MAXIMIZE)
    m.setParam('OutputFlag', 1)
    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]=="z":
               solution.append(int(v.varName[1:]))
    opt_sites = sites[solution]
    return opt_sites,m.objVal

parser = argparse.ArgumentParser(description='Solve MCLP instance by gurobi.')
parser.add_argument('file_path', type=str, help='Path to the file to be processed.')
args = parser.parse_args()

print(f'The file path is: {args.file_path}')

filepath = args.file_path
budget = 10
coverage_radius = 5.5
dformat = {'names': ('type', 'index', 'x', 'y', 'demand'), 'formats': ('S1', 'i4', 'f4', 'f4', 'i4')}
pt_type, pt_index, pt_x, pt_y, pt_demand = np.loadtxt(filepath, delimiter=None, skiprows = 1, dtype = dformat, unpack = True)
customer_mask = pt_type == b'C'
facility_mask = pt_type == b'F'

customer_pts = np.stack((pt_x[customer_mask], pt_y[customer_mask]), 1)
facility_pts = np.stack((pt_x[facility_mask], pt_y[facility_mask]), 1)

mclp(customer_pts, facility_pts, pt_demand, budget, coverage_radius)