Backup Covering Location Problem (BCLP)

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_BCLP_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.bclp([coverage1, coverage2], max_supply={coverage1: 5, coverage2: 19})

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_BCLP_14_0.png

Find the total demand that is covered

[8]:
covered_demand = d.query(f"GEOID10 in ({[f'{i}' for i in problem.selected_demand(coverage1)]})")
print(f'{(covered_demand["Population"].sum() / d["Population"].sum()*100):0.2f}%')
100.00%

Find the demand that is provided backup coverage

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

Find total number of facilities

[10]:
print(len(selected_locations) + len(selected_locations2))
24