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);
![../_images/examples_MCLP_6_0.png](../_images/examples_MCLP_6_0.png)
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);
![../_images/examples_MCLP_14_0.png](../_images/examples_MCLP_14_0.png)
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))
15