Location Set Covering Problem (LSCP)

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")
coverage2 = Coverage.from_geodataframes(d, s2, "GEOID10", "ORIG_ID")

Create the problem

[5]:

problem = Problem.lscp([coverage1, coverage2])

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 number of locations required

[8]:

print(problem.pulp_problem.objective.value())