Empirical examples of facility location problems¶

This tutorial aims to show a facility location application. To achieve this goal the tuorial will solve mulitple real world facility location problems using a dataset describing an area of census tract 205, San Francisco. The general scenarios can be stated as: store sites should supply the demand in this census tract considering the distance between the two sets of sites: demand points and candidate supply sites. Four fundamental facility location models are utilized to highlight varying outcomes dependent on objectives: LSCP (Location Set Covering Problem), MCLP (Maximal Coverage Location Problem), P-Median and P-Center. For further information on these models, it’s recommended to see the notebooks that explain more deeply about each one: LSCP , MCLP , P-Median , P-Center .

Also, this tutorial demonstrates the use of different solvers that PULP supports.

%config InlineBackend.figure_format = "retina" %load_ext watermark %watermark

Last updated: 2023-12-10T13:35:40.591294-05:00 Python implementation: CPython Python version : 3.12.0 IPython version : 8.18.0 Compiler : Clang 15.0.7 OS : Darwin Release : 23.1.0 Machine : x86_64 Processor : i386 CPU cores : 8 Architecture: 64bit

import geopandas import matplotlib.pyplot as plt from matplotlib.patches import Patch import matplotlib.lines as mlines import matplotlib_scalebar from matplotlib_scalebar.scalebar import ScaleBar import numpy import pandas import pulp import shapely from shapely.geometry import Point import spopt import time import warnings %watermark -w %watermark -iv

Watermark: 2.4.3 pulp : 2.7.0 matplotlib : 3.8.2 geopandas : 0.14.1 numpy : 1.26.2 pandas : 2.1.3 shapely : 2.0.2 spopt : 0.5.1.dev53+g5cadae7 matplotlib_scalebar: 0.8.1

We use 4 data files as input: - network_distance is the distance between facility candidate sites calculated by ArcGIS Network Analyst Extension - demand_points represents the demand points with some important features for the facility location problem like population - facility_points represents the stores that are candidate facility sites - tract is the polygon of census tract 205.

All datasets are online on this repository.

DIRPATH = "../spopt/tests/data/"

network_distance = pandas.read_csv( DIRPATH + "SF_network_distance_candidateStore_16_censusTract_205_new.csv" ) network_distance

3280 rows × 4 columns

demand_points = pandas.read_csv( DIRPATH + "SF_demand_205_centroid_uniform_weight.csv", index_col=0 ) demand_points = demand_points.reset_index(drop=True) demand_points

205 rows × 11 columns

facility_points = pandas.read_csv(DIRPATH + "SF_store_site_16_longlat.csv", index_col=0) facility_points = facility_points.reset_index(drop=True) facility_points

study_area = geopandas.read_file(DIRPATH + "ServiceAreas_4.shp").dissolve() study_area

base = study_area.plot() base.axis("off");

To start modeling the problem assuming the arguments expected by spopt.locate , we should pass a numpy 2D array as a cost matrix. So, first we pivot the network_distance dataframe.

Note that the columns and rows are sorted in ascending order.

ntw_dist_piv = network_distance.pivot_table( values="distance", index="DestinationName", columns="name" ) ntw_dist_piv

205 rows × 16 columns

Here the pivot table is transformed to a numpy 2D array such as spopt.locate expected. The matrix has a shape of 205x16.

cost_matrix = ntw_dist_piv.to_numpy() cost_matrix.shape

(205, 16)

cost_matrix[:3, :3]

array([[11495.19045438, 20022.66650296, 10654.59373325], [10436.16991032, 19392.09477041, 10024.0220007 ], [10746.29681106, 19404.67285964, 10036.60008993]])

Now, as the rows and columns of our cost matrix are sorted, we have to sort our facility points and demand points geodataframes, too.

n_dem_pnts = demand_points.shape[0] n_fac_pnts = facility_points.shape[0]

process = lambda df: as_gdf(df).sort_values(by=["NAME"]).reset_index(drop=True) as_gdf = lambda df: geopandas.GeoDataFrame(df, geometry=pnts(df)) pnts = lambda df: geopandas.points_from_xy(df.long, df.lat)

facility_points = process(facility_points) demand_points = process(demand_points)

Reproject the input spatial data.

for _df in [facility_points, demand_points, study_area]: _df.set_crs("EPSG:4326", inplace=True) _df.to_crs("EPSG:7131", inplace=True)

Here the rest of parameters are set.

ai is the demand weight and in this case we model as population in year 2000 of this tract.

ai = demand_points["POP2000"].to_numpy() # maximum service radius (in meters) SERVICE_RADIUS = 5000 # number of candidate facilities in optimal solution P_FACILITIES = 4

Below, the method is used to plot the results of the four models that we will prepare to solve the problem.

dv_colors_arr = [ "darkcyan", "mediumseagreen", "saddlebrown", "darkslategray", "lightskyblue", "thistle", "lavender", "darkgoldenrod", "peachpuff", "coral", "mediumvioletred", "blueviolet", "fuchsia", "cyan", "limegreen", "mediumorchid", ] dv_colors = ": dv_colors_arr[i] for i in range(len(dv_colors_arr))> dv_colors

def plot_results(model, p, facs, clis=None, ax=None): """Visualize optimal solution sets and context.""" if not ax: multi_plot = False fig, ax = plt.subplots(figsize=(6, 9)) markersize, markersize_factor = 4, 4 else: ax.axis("off") multi_plot = True markersize, markersize_factor = 2, 2 ax.set_title(model.name, fontsize=15) # extract facility-client relationships for plotting (except for p-dispersion) plot_clis = isinstance(clis, geopandas.GeoDataFrame) if plot_clis: cli_points = <> fac_sites = <> for i, dv in enumerate(model.fac_vars): if dv.varValue: dv, predef = facs.loc[i, ["dv", "predefined_loc"]] fac_sites[dv] = [i, predef] if plot_clis: geom = clis.iloc[model.fac2cli[i]]["geometry"] cli_points[dv] = geom # study area and legend entries initialization study_area.plot(ax=ax, alpha=0.5, fc="tan", ec="k", zorder=1) _patch = Patch(alpha=0.5, fc="tan", ec="k", label="Dissolved Service Areas") legend_elements = [_patch] if plot_clis: # any clients that not asscociated with a facility if model.name.startswith("mclp"): c = "k" if model.n_cli_uncov: idx = [i for i, v in enumerate(model.cli2fac) if len(v) == 0] pnt_kws = dict(ax=ax, fc=c, ec=c, marker="s", markersize=7, zorder=2) clis.iloc[idx].plot(**pnt_kws) _label = f"Demand sites not covered ($n$=)" _mkws = dict(marker="s", markerfacecolor=c, markeredgecolor=c, linewidth=0) legend_elements.append(mlines.Line2D([], [], ms=3, label=_label, **_mkws)) # all candidate facilities facs.plot(ax=ax, fc="brown", marker="*", markersize=80, zorder=8) _label = f"Facility sites ($n$=)" _mkws = dict(marker="*", markerfacecolor="brown", markeredgecolor="brown") legend_elements.append(mlines.Line2D([], [], ms=7, lw=0, label=_label, **_mkws)) # facility-(client) symbology and legend entries zorder = 4 for fname, (fac, predef) in fac_sites.items(): cset = dv_colors[fname] if plot_clis: # clients geoms = cli_points[fname] gdf = geopandas.GeoDataFrame(geoms) gdf.plot(ax=ax, zorder=zorder, ec="k", fc=cset, markersize=100 * markersize) _label = f"Demand sites covered by " _mkws = dict(markerfacecolor=cset, markeredgecolor="k", ms=markersize + 7) legend_elements.append( mlines.Line2D([], [], marker="o", lw=0, label=_label, **_mkws) ) # facilities ec = "k" lw = 2 predef_label = "predefined" if model.name.endswith(predef_label) and predef: ec = "r" lw = 3 fname += f" ()" facs.iloc[[fac]].plot( ax=ax, marker="*", markersize=1000, zorder=9, fc=cset, ec=ec, lw=lw ) _mkws = dict(markerfacecolor=cset, markeredgecolor=ec, markeredgewidth=lw) legend_elements.append( mlines.Line2D([], [], marker="*", ms=20, lw=0, label=fname, **_mkws) ) # increment zorder up and markersize down for stacked client symbology zorder += 1 if plot_clis: markersize -= markersize_factor / p if not multi_plot: # legend kws = dict(loc="upper left", bbox_to_anchor=(1.05, 0.7)) plt.legend(handles=legend_elements, **kws)

LSCP¶

from spopt.locate import LSCP lscp = LSCP.from_cost_matrix(cost_matrix, SERVICE_RADIUS) lscp = lscp.solve(pulp.GLPK(msg=False)) lscp_objval = lscp.problem.objective.value() lscp_objval