Maximum Covering Location Problem (MCLP)

Import the required modules

[1]:

import warnings
warnings.filterwarnings('ignore')
import geopandas
import matplotlib.pyplot as plt
import pulp
from allagash import Coverage, Problem
%matplotlib inline

Load the data

[2]:

d = geopandas.read_file("sample_data/demand_point.shp")
s = geopandas.read_file("sample_data/facility_service_areas.shp")
s2 = geopandas.read_file("sample_data/facility2_service_areas.shp")

Plot the data

[3]:

fig, ax = plt.subplots(figsize=(15, 15))
ax.set_aspect('equal')
d.plot(ax=ax, color='blue')
s.plot(ax=ax, color='yellow', alpha=0.3)
s2.plot(ax=ax, color='green', alpha=0.3);

Create the coverage

[4]:

coverage1 = Coverage.from_geodataframes(d, s, "GEOID10", "ORIG_ID", demand_col="Population")
coverage2 = Coverage.from_geodataframes(d, s2, "GEOID10", "ORIG_ID", demand_col="Population")

Create the problem

[5]:

problem = Problem.mclp([coverage1, coverage2], max_supply={coverage1: 5, coverage2: 10})

Solve the problem

[6]:

problem.solve(pulp.GLPK(msg=False));

Plot the selected locations

[7]:

selected_locations = s.query(f"ORIG_ID in ({','.join(problem.selected_supply(coverage1))})")
selected_locations2 = s2.query(f"ORIG_ID in ({','.join(problem.selected_supply(coverage2))})")

fig, ax = plt.subplots(figsize=(15, 15))
ax.set_aspect('equal')
d.plot(ax=ax, color='blue')
selected_locations.plot(ax=ax, color='yellow', alpha=0.3)
selected_locations2.plot(ax=ax, color='green', alpha=0.3);

Find the total demand that is covered

[8]:

print(f'{(problem.pulp_problem.objective.value() / d["Population"].sum()*100):0.2f}%')

95.50%

Find total number of facilities

[9]:

print(len(selected_locations) + len(selected_locations2))