capacitated p-median with predefined facilities

spopt/locate/base.py

@@ -16,6 +16,10 @@
 }
 
 
+class SpecificationError(pulp.PulpError):
+    pass
+
+
 class LocateSolver(BaseSpOptExactSolver):
     """Base class for the ``locate`` package."""
 
@@ -449,37 +453,86 @@ def add_facility_constraint(obj: T_FacModel, p_facilities: int) -> None:
 
     @staticmethod
     def add_predefined_facility_constraint(
-        obj: T_FacModel, predefined_fac: np.array
+        obj: T_FacModel,
+        predefined_fac: np.array,
+        demand: np.array = None,
+        facility_capacity: np.array = None,
     ) -> None:
         """
-        Create predefined demand constraints.
+        Create predefined supply constraints.
 
         Parameters
         ----------
 
         obj : T_FacModel
             A bounded type of the ``LocateSolver`` class.
-        facility_indexes : numpy.array
+        predefined_fac : numpy.array
             Indexes of facilities that are already located (zero-indexed).
+        demand : numpy.array (default None)
+            A 1D array of service load or population demand.
+        facility_capacity : numpy.array (default None)
+            The capacity of each facility.
 
         Returns
         -------
 
         None
 
         """
+        if predefined_fac.ndim == 2:
+            n, k = predefined_fac.shape
+            if k != 1:
+                raise ValueError(
+                    "predefined facilties array must only be "
+                    "of shape (n_supply, 1) or (n_supply,)"
+                )
+            predefined_fac = predefined_fac.squeeze()
+
+        n_predefined = len(predefined_fac)
+
         if hasattr(obj, "fac_vars"):
             fac_vars = getattr(obj, "fac_vars")
-            for ind in range(len(predefined_fac)):
-                if predefined_fac[ind]:
-                    fac_vars[ind].setInitialValue(1)
-                    fac_vars[ind].fixValue()
+            n_facilities = len(fac_vars)
+
+            if n_facilities > n_predefined:  # treat as indices
+                dummies = np.zeros_like(fac_vars)
+                dummies[predefined_fac] = 1
+            elif n_facilities == n_predefined:  # treat as dummies
+                dummies = predefined_fac.copy()
+            else:
+                raise ValueError(
+                    "More preselected facilities were provided than supply sites. "
+                    "Expected fewer preselected facilities than supply sites. Check"
+                    " the shape of the predefined faciltiies & supply sites provided."
+                )
+
+            for i, dummy in enumerate(dummies):
+                if dummy:
+                    fac_vars[i].setInitialValue(1)
+                    fac_vars[i].fixValue()
+
         else:
             raise AttributeError(
                 "Before setting predefined facility constraints "
                 "facility variables must be set."
             )
 
+        # To add the capacity fulfill constraint
+        if (
+            (facility_capacity is not None)
+            and hasattr(obj, "cli_assgn_vars")
+            and hasattr(obj, "fac_vars")
+        ):
+            fac_vars = getattr(obj, "fac_vars")
+            cli_vars = getattr(obj, "cli_assgn_vars")
+            model = getattr(obj, "problem")
+
+            for j in predefined_fac:
+                model += (
+                    pulp.lpSum(demand[i] * cli_vars[i, j] for i in range(len(cli_vars)))
+                    == fac_vars[j] * facility_capacity[j]
+                )
+
     @staticmethod
     def add_facility_capacity_constraint(
         obj: T_FacModel,

spopt/locate/p_median.py

@@ -11,6 +11,7 @@
     FacilityModelBuilder,
     LocateSolver,
     MeanDistanceMixin,
+    SpecificationError,
 )
 
 
@@ -129,6 +130,8 @@ def from_cost_matrix(
         weights: np.array,
         p_facilities: int,
         predefined_facilities_arr: np.array = None,
+        facility_capacities: np.array = None,
+        fulfill_predefined_fac: bool = False,
         name: str = "p-median",
     ):
         """
@@ -145,6 +148,10 @@ def from_cost_matrix(
             The number of facilities to be located.
         predefined_facilities_arr : numpy.array (default None)
             Predefined facilities that must appear in the solution.
+        facility_capacity : numpy.array (default None)
+            The capacity of each facility.
+        fulfill_predefined_fac : bool (default False)
+            If the predefined facilities need to be fulfilled.
         name : str (default 'p-median')
             The problem name.
 
@@ -246,10 +253,50 @@ def from_cost_matrix(
         )
 
         if predefined_facilities_arr is not None:
-            FacilityModelBuilder.add_predefined_facility_constraint(
-                p_median, predefined_facilities_arr
+            if fulfill_predefined_fac and facility_capacities is not None:
+                sum_predefined_fac_cap = np.sum(
+                    facility_capacities[predefined_facilities_arr]
+                )
+                if sum_predefined_fac_cap <= weights.sum():
+                    FacilityModelBuilder.add_predefined_facility_constraint(
+                        p_median,
+                        predefined_facilities_arr,
+                        weights,
+                        facility_capacities,
+                    )
+                else:
+                    raise SpecificationError(
+                        "Problem is infeasible. The predefined facilities can't be "
+                        "fulfilled, because their capacity is larger than the total "
+                        f"demand {weights.sum()}."
+                    )
+            elif fulfill_predefined_fac and facility_capacities is None:
+                raise SpecificationError(
+                    "Data on the capacity of the facility is missing, "
+                    "so the model cannot be calculated."
+                )
+            else:
+                FacilityModelBuilder.add_predefined_facility_constraint(
+                    p_median, predefined_facilities_arr
+                )
+
+        if facility_capacities is not None:
+            sorted_capacities = np.sort(facility_capacities)
+            highest_possible_capacity = sorted_capacities[-p_facilities:].sum()
+            if highest_possible_capacity < weights.sum():
+                raise SpecificationError(
+                    "Problem is infeasible. The highest possible capacity "
+                    f"{highest_possible_capacity}, coming from the {p_facilities} "
+                    "sites with the highest capacity, is smaller than "
+                    f"the total demand {weights.sum()}."
+                )
+            FacilityModelBuilder.add_facility_capacity_constraint(
+                p_median,
+                weights,
+                facility_capacities,
+                range(len(weights)),
+                range(len(facility_capacities)),
             )
-
         p_median.__add_obj(r_cli, r_fac)
 
         FacilityModelBuilder.add_facility_constraint(p_median, p_facilities)
@@ -267,7 +314,9 @@ def from_geodataframe(
         facility_col: str,
         weights_cols: str,
         p_facilities: int,
+        facility_capacity_col: str = None,
         predefined_facility_col: str = None,
+        fulfill_predefined_fac: bool = False,
         distance_metric: str = "euclidean",
         name: str = "p-median",
     ):
@@ -294,6 +343,10 @@ def from_geodataframe(
            The number of facilities to be located.
         predefined_facility_col : str (default None)
             Column name representing facilities are already defined.
+        facility_capacities_col: str (default None)
+            Column name representing the capacities of each facility.
+        fulfill_predefined_fac : bool (default False)
+            If the predefined facilities need to be fulfilled.
         distance_metric : str (default 'euclidean')
             A metric used for the distance calculations supported by
             `scipy.spatial.distance.cdist <https://docs.scipy.org/doc/scipy/reference/generated/scipy.spatial.distance.cdist.html>`_.
@@ -379,6 +432,10 @@ def from_geodataframe(
         if predefined_facility_col is not None:
             predefined_facilities_arr = gdf_fac[predefined_facility_col].to_numpy()
 
+        facility_capacities = None
+        if facility_capacity_col is not None:
+            facility_capacities = gdf_fac[facility_capacity_col].to_numpy()
+
         service_load = gdf_demand[weights_cols].to_numpy()
         dem = gdf_demand[demand_col]
         fac = gdf_fac[facility_col]
@@ -410,7 +467,13 @@ def from_geodataframe(
         distances = cdist(dem_data, fac_data, distance_metric)
 
         return cls.from_cost_matrix(
-            distances, service_load, p_facilities, predefined_facilities_arr, name
+            cost_matrix=distances,
+            weights=service_load,
+            p_facilities=p_facilities,
+            predefined_facilities_arr=predefined_facilities_arr,
+            facility_capacities=facility_capacities,
+            fulfill_predefined_fac=fulfill_predefined_fac,
+            name=("capacitated" + name if facility_capacities is not None else name),
         )
 
     def facility_client_array(self) -> None:

spopt/tests/data/example_subject_schools.csv

@@ -0,0 +1,75 @@
+,SE2 PP: Code,SE2 PP: PC,PL: Subject,Count,priority
+0,IOE00043,NW51RL,Mathematics,1,2.0
+1,IOE00044,NW11RX,Mathematics,1,3.0
+2,IOE00045,NW23RT,Mathematics,1,3.0
+3,IOE00046,NW51UJ,Mathematics,1,3.0
+4,IOE00081,SE181QF,Mathematics,1,2.0
+5,IOE00128,E28LS,Mathematics,1,2.0
+6,IOE00131,E96NR,Mathematics,1,2.0
+7,IOE00166,W68RB,Mathematics,1,2.0
+8,IOE00172,SW64UN,Mathematics,1,3.0
+9,IOE00218,N65LY,Mathematics,4,3.0
+10,IOE00223,N43LS,Mathematics,1,3.0
+11,IOE00276,SW165UQ,Mathematics,1,2.0
+12,IOE00348,SE264RD,Mathematics,1,2.0
+13,IOE00352,SE146TJ,Mathematics,1,2.0
+14,IOE00403,SE166AT,Mathematics,1,2.0
+15,IOE00406,SE220AT,Mathematics,1,1.0
+16,IOE00452,E20PX,Mathematics,1,3.0
+17,IOE00454,E14SD,Mathematics,2,3.0
+18,IOE00460,E10LB,Mathematics,1,2.0
+19,IOE00471,E26NW,Mathematics,1,2.0
+20,IOE00506,SW128JZ,Mathematics,1,3.0
+21,IOE00509,SW170AQ,Mathematics,3,3.0
+22,IOE00538,NW65SN,Mathematics,1,2.0
+23,IOE00553,SW1E5HJ,Mathematics,1,1.0
+24,IOE00572,IG119AG,Mathematics,2,2.0
+25,IOE00579,RM94UN,Mathematics,3,2.0
+26,IOE00582,RM66SB,Mathematics,1,2.0
+27,IOE00641,NW117HY,Mathematics,1,2.0
+28,IOE00664,N20SE,Mathematics,1,2.0
+29,IOE00674,N20SQ,Mathematics,1,1.0
+30,IOE00692,DA67QJ,Mathematics,1,2.0
+31,IOE00694,DA159NU,Mathematics,1,1.0
+32,IOE00739,HA30UH,Mathematics,1,2.0
+33,IOE00740,NW99JR,Mathematics,1,2.0
+34,IOE00745,HA98NA,Mathematics,1,2.0
+35,IOE00753,SE94QF,Mathematics,1,2.0
+36,IOE00812,SE256AE,Mathematics,1,2.0
+37,IOE00813,SE192JH,Mathematics,1,2.0
+38,IOE00818,CR78BT,Mathematics,1,3.0
+39,IOE00822,SE193HL,Mathematics,1,2.0
+40,IOE00831,CR07NJ,Mathematics,1,2.0
+41,IOE00863,UB13HZ,Mathematics,1,3.0
+42,IOE00867,W71JJ,Mathematics,2,1.0
+43,IOE00869,W30HW,Mathematics,1,2.0
+44,IOE00953,N170PG,Mathematics,1,1.0
+45,IOE00954,N101NE,Mathematics,1,3.0
+46,IOE00963,N153QR,Mathematics,2,2.0
+47,IOE00976,N225HN,Mathematics,1,2.0
+48,IOE00978,N227ST,Mathematics,1,2.0
+49,IOE01009,HA55RP,Mathematics,1,2.0
+50,IOE01013,HA29AH,Mathematics,1,3.0
+51,IOE01030,RM79NX,Mathematics,1,2.0
+52,IOE01062,UB82PR,Mathematics,1,2.0
+53,IOE01153,E126JB,Mathematics,1,2.0
+54,IOE01157,E138SJ,Mathematics,2,2.0
+55,IOE01165,E65JG,Mathematics,1,2.0
+56,IOE01169,E63SQ,Mathematics,1,2.0
+57,IOE01190,IG89LA,Mathematics,1,3.0
+58,IOE01191,IG87DQ,Mathematics,2,3.0
+59,IOE01223,TW13BB,Mathematics,1,2.0
+60,IOE01263,E173PY,Mathematics,1,2.0
+61,IOE01271,E49PJ,Mathematics,2,3.0
+62,IOE01275,E47LT,Mathematics,1,2.0
+63,IOE01406,SL12PU,Mathematics,2,2.0
+64,IOE01431,IG103JA,Mathematics,1,2.0
+65,IOE01462,SS69BZ,Mathematics,1,2.0
+66,IOE01586,WD232TY,Mathematics,2,3.0
+67,IOE01595,AL36DB,Mathematics,1,3.0
+68,IOE01599,WD234PA,Mathematics,1,2.0
+69,IOE01821,N111BF,Mathematics,2,2.0
+70,IOE01839,W21QZ,Mathematics,1,2.0
+71,IOE01856,N28GA,Mathematics,1,1.0
+72,IOE02395,N166PA,Mathematics,1,2.0
+73,IOE02497,SE193UG,Mathematics,1,2.0