Authors: Germano Barcelos, James Gaboardi, Levi J. Wolf, Qunshan Zhao
Location Set Covering is a problem realized by Toregas, et al. (1971). He figured out that emergency services must have placed according to a response time, since, there is a allowable maximum service time when it's discussed how handle an emergency activity. Therefore he proprosed a model named LSCP that:
Minimize the number of facilities needed and locate them so that every demand area is covered within a predefined maximal service distance or time. Church L., Murray, A. (2018)
LSCP can be written as:
$\begin{array} \displaystyle \textbf{Minimize} & \sum_{j=1}^{n}{x_j} && (1) \\ \displaystyle \textbf{Subject to:} & \sum_{j\in N_i}{x_j} \geq 1 & \forall i & (2) \\ & x_j \in {0,1} & \forall j & (3) \\ \end{array}$
$\begin{array} \displaystyle \textbf{Where:}\\ & & \displaystyle i & \small = & \textrm{index referencing nodes of the network as demand} \\ & & j & \small = & \textrm{index referencing nodes of the network as potential facility sites} \\ & & S & \small = & \textrm{maximal acceptable service distance or time standard} \\ & & d_{ij} & \small = & \textrm{shortest distance or travel time between nodes} i \textrm{and} j \\ & & N_i & \small = & \{j | d_{ij} < S\} \\ & & x_j & \small = & \begin{cases} 1, \text{if a facility is located at node} j\\ 0, \text{otherwise} \\ \end{cases} \end{array}$
This excerpt above was quoted from Church L., Murray, A. (2018)
This tutorial solves LSCP using spopt.locate.coverage.LSCP
instance that depends on a array 2D representing the costs between facilities candidate sites and demand points. For that it uses a lattice 10x10 with simulated points to calculate the costs.
from spopt.locate.coverage import LSCP
from spopt.locate.util import simulated_geo_points
import numpy
import geopandas
import pulp
import spaghetti
from shapely.geometry import Point
import matplotlib.pyplot as plt
/home/gegen07/anaconda3/envs/spopt/lib/python3.9/site-packages/spaghetti/network.py:36: FutureWarning: The next major release of pysal/spaghetti (2.0.0) will drop support for all ``libpysal.cg`` geometries. This change is a first step in refactoring ``spaghetti`` that is expected to result in dramatically reduced runtimes for network instantiation and operations. Users currently requiring network and point pattern input as ``libpysal.cg`` geometries should prepare for this simply by converting to ``shapely`` geometries. warnings.warn(f"{dep_msg}", FutureWarning)
Since the model needs a distance cost matrix we should define some variables. In the comments, it's defined what these variables are for but solver. The solver, assigned below as pulp.PULP_CBC_CMD
, is an interface to optimization solver developed by COIN-OR. If you want to use another optimization interface as Gurobi or CPLEX see this guide that explains how to achieve this.
CLIENT_COUNT = 100 # quantity demand points
FACILITY_COUNT = 5 # quantity supply points
SERVICE_RADIUS = 8 # maximum service radius in meters
# Random seeds for reproducibility
CLIENT_SEED = 5
FACILITY_SEED = 6
solver = pulp.PULP_CBC_CMD(msg=False, warmStart=True)
Create lattice 10x10 with 9 vertical lines in interior.
lattice = spaghetti.regular_lattice((0, 0, 10, 10), 9, exterior=True)
ntw = spaghetti.Network(in_data=lattice)
Transform spaghetti instance into geodataframe.
street = spaghetti.element_as_gdf(ntw, arcs=True)
street_buffered = geopandas.GeoDataFrame(
geopandas.GeoSeries(street["geometry"].buffer(0.2).unary_union),
crs=street.crs,
columns=["geometry"],
)
Plotting the network created by spaghetti we can verify that it seems a district with quarters and streets.
street.plot()
<AxesSubplot:>
The function simulated_geo_points
simulates points inside a network. In this case, it uses a lattice network 10x10 created by using spaghetti package.
Below we use the function defined above and simulate the points inside lattice bounds.
client_points = simulated_geo_points(street_buffered, needed=CLIENT_COUNT, seed=CLIENT_SEED)
facility_points = simulated_geo_points(
street_buffered, needed=FACILITY_COUNT, seed=FACILITY_SEED
)
Plotting the 100 client and 5 facility points we can see that the function generates dummy points to an area of 10x10 which is the area created by our lattice created on previous cells.
fig, ax = plt.subplots(figsize=(6, 6))
street.plot(ax=ax, alpha=0.8, zorder=1, label='streets')
facility_points.plot(ax=ax, color='red', zorder=2, label='facility candidate sites ($n$=5)')
client_points.plot(ax=ax, color='black', label='clients points ($n$=100)')
plt.legend(loc='upper left', bbox_to_anchor=(1.05, 1))
<matplotlib.legend.Legend at 0x7fb35ddcffa0>
To use cost matrix or geodataframes we have to pay attention in some details. The client and facility points simulated don't belong to network, so if we calculate the distances now we are supposed to receive a wrong result. Before calculating distances we snap points to the networok and then calculate the distances.
Below we snap points that is not spatially belong to network and create new real points geodataframes
ntw.snapobservations(client_points, "clients", attribute=True)
clients_snapped = spaghetti.element_as_gdf(
ntw, pp_name="clients", snapped=True
)
ntw.snapobservations(facility_points, "facilities", attribute=True)
facilities_snapped = spaghetti.element_as_gdf(
ntw, pp_name="facilities", snapped=True
)
Now the plot seems more organized as the points belong to network. The network created is plotted below with facility points and clients points:
fig, ax = plt.subplots(figsize=(6, 6))
street.plot(ax=ax, alpha=0.8, zorder=1, label='streets')
facilities_snapped.plot(ax=ax, color='red', zorder=2, label='facility candidate sites ($n$=5)')
clients_snapped.plot(ax=ax, color='black', label='clients points ($n$=100)')
plt.legend(loc='upper left', bbox_to_anchor=(1.05, 1))
<matplotlib.legend.Legend at 0x7fb36771a400>
Calculate distance between clients and facilities.
cost_matrix = ntw.allneighbordistances(
sourcepattern=ntw.pointpatterns["clients"],
destpattern=ntw.pointpatterns["facilities"],
)
The expected result here is a Dijkstra distance between clients and facilities points, so we our case an array 2D 100x5.
cost_matrix
array([[12.60302601, 3.93598651, 8.16571655, 6.04319467, 5.65607701], [13.10096347, 4.43392397, 8.66365401, 6.54113213, 5.15813955], [ 6.9095462 , 4.2425067 , 2.47223674, 0.34971486, 5.34955682], [ 2.98196832, 7.84581224, 3.45534114, 3.57786302, 6.25374871], [ 7.5002892 , 6.32806975, 4.55779979, 6.43527791, 11.75939222], [ 0.60209077, 11.42987132, 5.03940023, 7.16192211, 9.8378078 ], [ 5.37335867, 6.20113923, 2.43086927, 4.30834738, 9.6324617 ], [ 5.40801577, 5.41976478, 3.02929369, 1.15181557, 4.85108725], [ 3.68807115, 8.51585171, 2.12538061, 4.24790249, 7.94717417], [14.22503627, 4.60274429, 9.78772681, 7.66520493, 4.98931924], [10.32521229, 4.99225179, 7.38272288, 9.260201 , 14.58431531], [ 6.65436171, 7.98732222, 5.59685112, 3.719373 , 2.58135531], [11.55510375, 1.11193575, 7.11779429, 5.37988496, 10.70399927], [10.90832519, 1.75871431, 6.47101573, 6.02666352, 11.35077783], [ 9.29354019, 9.53424036, 7.14376926, 5.26629115, 0.05782317], [11.25279502, 3.57498553, 6.81548556, 4.69296368, 6.01707799], [ 6.14400601, 11.47696651, 9.08649542, 7.2090173 , 3.09171102], [10.43008909, 2.23695041, 5.99277963, 6.50489962, 11.82901393], [ 1.79838406, 11.13134457, 4.74087347, 6.86339535, 9.53928104], [ 2.93052752, 7.89725303, 3.50678194, 3.62930382, 6.30518951], [11.55272282, 6.21976231, 8.61023341, 10.48771153, 15.81182584], [ 8.83964081, 3.66742137, 5.89715141, 7.77462952, 13.09874384], [ 4.11777697, 9.45073748, 7.06026638, 5.18278826, 7.11794005], [ 8.69768642, 8.63527408, 5.75519701, 7.63267513, 12.95678945], [ 8.2652832 , 6.56249735, 4.79222739, 2.66970551, 3.02956617], [ 1.71437731, 9.6185832 , 3.2281121 , 5.35063398, 8.02651967], [ 4.30308213, 6.52469842, 2.13422733, 2.25674921, 5.95602089], [ 9.31612329, 8.64908379, 6.25861269, 4.38113458, 0.94297974], [ 2.86540683, 13.69318738, 7.30271629, 9.42523817, 12.10112386], [ 8.95995574, 2.29291624, 4.52264628, 6.4001244 , 11.72423871], [10.54288208, 7.87584258, 6.10557262, 3.98305074, 1.71622094], [ 8.58885878, 8.74410173, 5.64636937, 7.52384749, 12.8479618 ], [ 2.51163835, 12.82132215, 6.43085106, 8.55337294, 11.22925863], [ 5.19213144, 5.63564912, 0.75482198, 1.74285727, 7.06697159], [ 4.1276352 , 13.2053253 , 6.81485421, 8.93737609, 11.61326178], [ 3.99217608, 6.83560448, 0.44513338, 2.94281263, 6.75645905], [ 5.88198594, 11.21494644, 8.82447535, 6.94699723, 5.35373109], [ 8.24225403, 4.58552653, 3.80494457, 1.68242269, 5.006537 ], [10.89255004, 6.22551054, 6.45524058, 4.3327187 , 3.36655299], [ 6.58504851, 11.91800902, 9.52753792, 7.6500598 , 2.65066851], [ 5.44204086, 8.77500136, 6.38453026, 4.50705215, 3.79367617], [ 5.56289993, 7.26488062, 4.87440953, 2.99693141, 3.6728171 ], [ 7.96716366, 10.86061689, 8.4701458 , 6.59266768, 1.26855337], [ 7.9603294 , 5.3726311 , 5.01783999, 6.89531811, 12.21943243], [ 8.68198919, 4.65097132, 5.73949978, 7.6169779 , 12.94109221], [ 9.06064716, 8.39360767, 6.00313657, 4.12565845, 1.19845586], [15.325265 , 4.65822551, 10.88795554, 8.76543366, 6.08954798], [ 3.51444772, 7.81851278, 1.95175718, 3.92572094, 7.77355074], [ 3.33469883, 14.16247938, 7.77200828, 9.89453017, 12.57041585], [ 4.46482284, 6.36295772, 1.40731225, 2.0950085 , 5.79428018], [11.20742649, 1.459613 , 6.77011704, 5.72756222, 11.05167653], [11.15442417, 5.67335639, 6.71711471, 4.59459283, 3.91870714], [ 5.17021584, 5.65756471, 0.73290638, 2.6103845 , 7.93449881], [ 5.54400588, 10.87696639, 8.48649529, 6.60901717, 5.28490286], [ 5.28695668, 8.04600382, 2.34446727, 4.22194539, 9.5460597 ], [ 7.33259845, 6.66555896, 4.27508786, 2.39760974, 2.92650457], [ 8.08642618, 10.74135437, 8.35088328, 6.47340516, 1.14929085], [ 7.97403829, 2.85374226, 3.53672884, 4.96095042, 10.28506473], [ 5.04455411, 6.2884064 , 2.1020647 , 3.97954282, 9.30365713], [ 8.05520721, 3.2777533 , 5.1127178 , 6.99019592, 12.31431023], [ 8.033197 , 3.2997635 , 5.09070759, 6.96818571, 12.29230002], [ 4.88391014, 5.94387041, 3.55339931, 1.6759212 , 4.62480712], [ 3.38092176, 9.44685879, 6.32341117, 5.17890958, 7.85479527], [ 5.83945489, 5.17241539, 2.78194429, 0.90446618, 4.41964814], [10.25764123, 4.57013932, 5.82033178, 3.69780989, 5.02192421], [ 3.16471551, 8.168245 , 1.7777739 , 3.90029578, 7.59956747], [ 8.83620663, 8.49675387, 5.89371722, 7.77119534, 13.09530965], [ 7.60754658, 6.94050708, 4.55003599, 2.67255787, 2.65155644], [ 4.14555919, 9.4785197 , 7.0880486 , 5.21057048, 5.09015784], [ 7.24126831, 4.57422881, 2.80395885, 0.68143697, 5.01783472], [ 5.70322513, 8.53100569, 2.76073572, 4.63821384, 9.96232815], [ 9.27617639, 9.55160416, 7.16113307, 5.28365495, 0.04045936], [ 2.5651854 , 11.39296595, 5.00249486, 7.12501674, 9.80090243], [14.22296519, 3.5559257 , 9.78565573, 7.66313385, 6.03613783], [ 8.33806089, 2.48971967, 3.90075143, 5.77822955, 11.10234386], [14.34079531, 3.51301476, 9.90348585, 7.78096397, 7.91830771], [ 7.55811406, 6.89107456, 4.50060346, 2.62312535, 2.70098897], [ 9.54667188, 8.87963238, 6.48916129, 4.61168317, 0.71243114], [ 6.99771477, 3.83006578, 2.56040532, 2.43788343, 7.76199775], [10.85478728, 4.18774778, 6.41747782, 4.29495594, 5.40431574], [ 6.89563349, 8.43732701, 3.95314408, 5.8306222 , 11.15473651], [12.29945454, 3.63241504, 7.86214508, 5.7396232 , 5.95964848], [ 6.57929244, 6.75366806, 3.63680304, 5.51428115, 10.83839547], [ 8.35675866, 8.47102189, 6.0805508 , 4.20307268, 1.90234436], [11.26183 , 1.40520949, 6.82452055, 5.67315871, 10.99727302], [ 6.92663397, 8.25959447, 5.86912337, 3.99164526, 2.30908306], [ 6.97410775, 3.8536728 , 2.53679829, 3.96088096, 9.28499527], [10.00715257, 8.82062799, 6.94964198, 4.92783614, 0.77143554], [ 8.83013405, 7.9976465 , 5.77262346, 3.89514534, 1.59441703], [ 6.69445759, 4.63850292, 2.86823295, 4.74571107, 10.06982539], [ 2.60649588, 11.43427644, 5.04380534, 7.16632722, 9.84221291], [ 9.01806225, 4.31489826, 6.07557284, 7.95305096, 13.27716527], [ 7.49191577, 3.84104474, 4.07077477, 5.94825289, 11.2723672 ], [ 7.80056437, 9.53239613, 4.85807497, 6.73555308, 12.0596674 ], [ 8.85156915, 2.48139135, 4.71112139, 6.58859951, 11.91271382], [10.04988811, 2.61715138, 5.61257866, 6.8851006 , 12.20921491], [ 3.68039673, 7.65256378, 1.26209268, 3.38461456, 7.08388625], [10.04984807, 7.28311243, 7.10735867, 8.98483678, 14.3089511 ], [ 8.34309643, 10.48468413, 8.09421303, 6.21673491, 0.8926206 ], [14.48203148, 3.65425093, 10.04472202, 7.92220014, 7.77707154]])
With LSCP.from_cost_matrix
we model LSC problem to cover all demand points with $p$ facility points within max_coverage
meters as service radius using cost matrix calculated previously.
lscp_from_cost_matrix = LSCP.from_cost_matrix(cost_matrix, SERVICE_RADIUS)
lscp_from_cost_matrix = lscp_from_cost_matrix.solve(solver)
Expected result is an instance of LSCP.
lscp_from_cost_matrix
<spopt.locate.coverage.LSCP at 0x7fb35db4d3a0>
Assigning predefined location using a geodataframe column
facilities_snapped['predefined_loc'] = numpy.array([0, 0, 0, 0, 1])
facilities_snapped
id | geometry | comp_label | predefined_loc | |
---|---|---|---|---|
0 | 0 | POINT (9.00000 3.25259) | 0 | 0 |
1 | 1 | POINT (0.91963 6.00000) | 0 | 0 |
2 | 2 | POINT (5.31010 4.00000) | 0 | 0 |
3 | 3 | POINT (5.18758 6.00000) | 0 | 0 |
4 | 4 | POINT (6.51169 10.00000) | 0 | 1 |
With LSCP.from_geodataframe
we model the LSC problem to cover all demand points with $p$ facility points within max_coverage
meters as service radius using geodataframes without calculating the cost matrix previously.
lscp_from_geodataframe = LSCP.from_geodataframe(
clients_snapped, facilities_snapped, "geometry", "geometry", SERVICE_RADIUS, distance_metric="euclidean"
)
lscp_from_geodataframe = lscp_from_geodataframe.solve(solver)
Expected result is an instance of LSCP.
lscp_from_geodataframe
<spopt.locate.coverage.LSCP at 0x7fb35db3a040>
Modelling LSCP with preselected facilities
lscp_preselected_from_geodataframe = LSCP.from_geodataframe(
clients_snapped, facilities_snapped, "geometry", "geometry", SERVICE_RADIUS, predefined_facility_col="predefined_loc", distance_metric="euclidean"
)
lscp_preselected_from_geodataframe = lscp_preselected_from_geodataframe.solve(solver)
The cell below describe the plotting of the results. For each method from LSCP class (from_cost_matrix, from_geodataframe) there is a plot displaying the facility site that was selected with a star colored and the points covered with the same color. Sometimes the demand points will be colored with not expected colors, it represents the coverage overlapping.
from matplotlib.patches import Patch
import matplotlib.lines as mlines
dv_colors_arr = [
"darkcyan",
"mediumseagreen",
"cyan",
"darkslategray",
"lightskyblue",
"limegreen",
"darkgoldenrod",
"peachpuff",
"coral",
"mediumvioletred",
"blueviolet",
"fuchsia",
"thistle",
"lavender",
"saddlebrown",
]
dv_colors = { f"y{i}":dv_colors_arr[i] for i in range(len(dv_colors_arr))}
def plot_results(model, facility_points):
arr_points = []
fac_sites = []
for i in range(FACILITY_COUNT):
if model.fac2cli[i]:
geom = client_points.iloc[model.fac2cli[i]]['geometry']
arr_points.append(geom)
fac_sites.append(i)
fig, ax = plt.subplots(figsize=(6, 6))
legend_elements = []
street.plot(ax=ax, alpha=1, color='black', zorder=1)
legend_elements.append(mlines.Line2D(
[],
[],
color='black',
label='streets',
))
facility_points.plot(ax=ax, color='brown', marker="*", markersize=80, zorder=2)
legend_elements.append(mlines.Line2D(
[],
[],
color='brown',
marker="*",
linewidth=0,
label=f'facility sites ($n$={FACILITY_COUNT})'
))
for i in range(len(arr_points)):
gdf = geopandas.GeoDataFrame(arr_points[i])
l = f"y{fac_sites[i]}"
label = f"coverage_points by y{fac_sites[i]}"
legend_elements.append(Patch(facecolor=dv_colors[l], edgecolor="k", label=label))
gdf.plot(ax=ax, zorder=3, alpha=0.7, edgecolor="k", color=dv_colors[l], label=label)
facility_points.iloc[[fac_sites[i]]].plot(ax=ax,
marker="*",
markersize=200 * 3.0,
alpha=0.8,
zorder=4,
edgecolor="k",
facecolor=dv_colors[l])
legend_elements.append(mlines.Line2D(
[],
[],
color=dv_colors[l],
marker="*",
ms=20 / 2,
markeredgecolor="k",
linewidth=0,
alpha=0.8,
label=f"y{fac_sites[i]} facility selected",
))
plt.title("LSCP", fontweight="bold")
plt.legend(handles = legend_elements, loc='upper left', bbox_to_anchor=(1.05, 1))
lscp_from_cost_matrix.facility_client_array()
plot_results(lscp_from_cost_matrix, facility_points)
lscp_from_geodataframe.facility_client_array()
plot_results(lscp_from_geodataframe, facility_points)
You may notice that the models are different. This result is expected as the distance between facility and demand points is calculated with different metrics. The cost matrix is calculated with dijkstra distance while the distance using geodataframe is calculated with euclidean distance.
But why it needs just one facility point to cover all of those demand points? It can be explained by the nature of the problem. The problem was configured in a synthetic manner, the street is created with 10x10 lattice and the max_coverage parameter is 8 meters, so this result is not weird at all. You can change the max_coverage parameter to 2 meters and you will obtain a different result but be aware with how many points will be covered.
lscp_preselected_from_geodataframe.facility_client_array()
plot_results(lscp_preselected_from_geodataframe, facility_points)