diff --git a/spopt/locate/base.py b/spopt/locate/base.py
index f55bfa42..b200518a 100644
--- a/spopt/locate/base.py
+++ b/spopt/locate/base.py
@@ -178,11 +178,11 @@ def add_facility_integer_variable(
         Parameters
         ----------
 
-        obj: T_FacModel
+        obj : T_FacModel
             A bounded type of the ``LocateSolver`` class.
-        range_facility: range
+        range_facility : range
             The range of facility points.
-        var_name: str
+        var_name : str
             A formatted string for the facility variable name.
 
         Returns
@@ -209,11 +209,11 @@ def add_client_integer_variable(
         Parameters
         ----------
 
-        obj: T_FacModel
+        obj : T_FacModel
             A bounded type of the ``LocateSolver`` class.
-        range_client: range
+        range_client : range
             The range of demand points.
-        var_name: str
+        var_name : str
             A formatted string for the demand variable name.
 
         Returns
@@ -244,15 +244,15 @@ def add_client_assign_integer_variable(
         Parameters
         ----------
 
-        obj: T_FacModel
+        obj : T_FacModel
             A bounded type of the ``LocateSolver`` class.
-        range_client: range
+        range_client : range
             The range of demand points.
-        range_facility: range
+        range_facility : range
             The range of facility points.
-        var_name: str
+        var_name : str
             A formatted string for the  client assignment variable name.
-        lp_category: pulp.LpVariable parameter
+        lp_category : pulp.LpVariable parameter
             The category this variable is in,
             ``pulp.LpInteger`` or ``pulp.LpContinuous``.
 
@@ -263,15 +263,20 @@ def add_client_assign_integer_variable(
 
         """
 
-        cli_assgn_vars = [
+        cli_assgn_vars = np.array(
             [
-                pulp.LpVariable(
-                    var_name.format(i=i, j=j), lowBound=0, upBound=1, cat=lp_category
-                )
-                for j in range_facility
+                [
+                    pulp.LpVariable(
+                        var_name.format(i=i, j=j),
+                        lowBound=0,
+                        upBound=1,
+                        cat=lp_category,
+                    )
+                    for j in range_facility
+                ]
+                for i in range_client
             ]
-            for i in range_client
-        ]
+        )
 
         setattr(obj, "cli_assgn_vars", cli_assgn_vars)
 
@@ -281,7 +286,7 @@ def add_weight_continuous_variable(obj: T_FacModel) -> None:
 
         Parameters
         ----------
-        obj: T_FacModel
+        obj : T_FacModel
             A bounded type of the ``LocateSolver`` class.
 
         Returns
@@ -301,7 +306,7 @@ def add_maximized_min_variable(obj: T_FacModel) -> None:
         Parameters
         ----------
 
-        obj: T_FacModel
+        obj : T_FacModel
             A bounded type of the ``LocateSolver`` class.
 
         Returns
@@ -316,10 +321,8 @@ def add_maximized_min_variable(obj: T_FacModel) -> None:
     @staticmethod
     def add_set_covering_constraint(
         obj: T_FacModel,
-        model: pulp.LpProblem,
-        ni: np.array,
-        range_facility: range,
         range_client: range,
+        range_facility: range,
     ) -> None:
         """
         Create set covering constraints.
@@ -329,17 +332,12 @@ def add_set_covering_constraint(
         Parameters
         ----------
 
-        obj: T_FacModel
+        obj : T_FacModel
             A bounded type of the ``LocateSolver`` class.
-        model: pulp.LpProblem
-            A ``pulp`` instance of an optimization model.
-        ni: numpy.array
-            A 2D array that defines candidate sites between facility
-            points within a distance to supply demand points.
-        range_facility: range
-            The range of facility points.
-        range_client: range
+        range_client : range
             The range of demand points.
+        range_facility : range
+            The range of facility points.
 
         Returns
         -------
@@ -349,9 +347,11 @@ def add_set_covering_constraint(
         """
         if hasattr(obj, "fac_vars"):
             fac_vars = getattr(obj, "fac_vars")
+            model = getattr(obj, "problem")
+            ni = getattr(obj, "aij")
             for i in range_client:
                 model += (
-                    pulp.lpSum([ni[i][j] * fac_vars[j] for j in range_facility]) >= 1
+                    pulp.lpSum([ni[i, j] * fac_vars[j] for j in range_facility]) >= 1
                 )
         else:
             raise AttributeError(
@@ -361,8 +361,6 @@ def add_set_covering_constraint(
     @staticmethod
     def add_backup_covering_constraint(
         obj: T_FacModel,
-        model: pulp.LpProblem,
-        ni: np.array,
         range_facility: range,
         range_client: range,
     ) -> None:
@@ -374,16 +372,11 @@ def add_backup_covering_constraint(
         Parameters
         ----------
 
-        obj: T_FacModel
+        obj : T_FacModel
             A bounded type of the ``LocateSolver`` class.
-        model: pulp.LpProblem
-            A ``pulp`` instance of an optimization model.
-        ni: numpy.array
-            A 2D array that defines candidate sites between facility
-            points within a distance to supply demand points.
-        range_facility: range
+        range_facility : range
             The range of facility points.
-        range_client: range
+        range_client : range
             The range of demand points.
 
         Returns
@@ -395,18 +388,20 @@ def add_backup_covering_constraint(
         if hasattr(obj, "fac_vars"):
             fac_vars = getattr(obj, "fac_vars")
             cli_vars = getattr(obj, "cli_vars")
+            model = getattr(obj, "problem")
+            ni = getattr(obj, "aij")
             for i in range_client:
                 if sum(ni[i]) >= 2:
                     model += (
                         pulp.lpSum(
-                            [int(ni[i][j]) * fac_vars[j] for j in range_facility]
+                            [int(ni[i, j]) * fac_vars[j] for j in range_facility]
                         )
                         >= 1 + 1 * cli_vars[i]
                     )
                 else:
                     model += (
                         pulp.lpSum(
-                            [int(ni[i][j]) * fac_vars[j] for j in range_facility]
+                            [int(ni[i, j]) * fac_vars[j] for j in range_facility]
                         )
                         >= 1 + 0 * cli_vars[i]
                     )
@@ -417,9 +412,7 @@ def add_backup_covering_constraint(
             )
 
     @staticmethod
-    def add_facility_constraint(
-        obj: T_FacModel, model: pulp.LpProblem, p_facilities: int
-    ) -> None:
+    def add_facility_constraint(obj: T_FacModel, p_facilities: int) -> None:
         """
         Create the facility constraint.
 
@@ -428,11 +421,9 @@ def add_facility_constraint(
         Parameters
         ----------
 
-        obj: T_FacModel
+        obj : T_FacModel
             A bounded type of the ``LocateSolver`` class.
-        model: pulp.LpProblem
-            A ``pulp`` instance of an optimization model.
-        p_facilities: int
+        p_facilities : int
             The number of facilities to be sited.
 
         Returns
@@ -443,6 +434,7 @@ def add_facility_constraint(
         """
         if hasattr(obj, "fac_vars"):
             fac_vars = getattr(obj, "fac_vars")
+            model = getattr(obj, "problem")
             model += pulp.lpSum(fac_vars) == p_facilities
         else:
             raise AttributeError(
@@ -451,7 +443,7 @@ def add_facility_constraint(
 
     @staticmethod
     def add_predefined_facility_constraint(
-        obj: T_FacModel, model: pulp.LpProblem, predefined_fac: np.array
+        obj: T_FacModel, predefined_fac: np.array
     ) -> None:
         """
         Create predefined demand constraints.
@@ -459,11 +451,9 @@ def add_predefined_facility_constraint(
         Parameters
         ----------
 
-        obj: T_FacModel
+        obj : T_FacModel
             A bounded type of the ``LocateSolver`` class.
-        model: pulp.LpProblem
-            A ``pulp`` instance of an optimization model.
-        facility_indexes: numpy.array
+        facility_indexes : numpy.array
             Indexes of facilities that are already located (zero-indexed).
 
         Returns
@@ -474,6 +464,7 @@ def add_predefined_facility_constraint(
         """
         if hasattr(obj, "fac_vars"):
             fac_vars = getattr(obj, "fac_vars")
+            model = getattr(obj, "problem")
             for ind in range(len(predefined_fac)):
                 if predefined_fac[ind]:
                     fac_vars[ind].setInitialValue(1)
@@ -486,7 +477,7 @@ def add_predefined_facility_constraint(
 
     @staticmethod
     def add_facility_capacity_constraint(
-        obj: T_FacModel, model, cl_ni, dq_ni, range_facility, range_client
+        obj: T_FacModel, dq_ni, cl_ni, range_client, range_facility
     ) -> None:
         """
         set facility capacity constraint:
@@ -496,30 +487,31 @@ def add_facility_capacity_constraint(
 
         Parameters
         ----------
-        obj: T_FacModel
+        obj : T_FacModel
             bounded type of LocateSolver class
-        model: pulp.LpProblem
-            optimization model problem
-        cl_ni: np.array
-            one-dimensional array that defines capacity limits of facility points
-        dq_ni: np.array
+        dq_ni : np.array
             one-dimensional array that defines demand quantities for demand points
-        range_facility: range
-            range of facility points quantity
-        range_client: range
+        cl_ni : np.array
+            one-dimensional array that defines capacity limits of facility points
+        range_client : range
             range of demand points quantity
+        range_facility : range
+            range of facility points quantity
 
         Returns
         -------
+
         None
+
         """
         if hasattr(obj, "fac_vars") and hasattr(obj, "cli_assgn_vars"):
             fac_vars = getattr(obj, "fac_vars")
             cli_assn_vars = getattr(obj, "cli_assgn_vars")
+            model = getattr(obj, "problem")
 
             for j in range_facility:
                 model += (
-                    pulp.lpSum([dq_ni[i] * cli_assn_vars[i][j] for i in range_client])
+                    pulp.lpSum([dq_ni[i] * cli_assn_vars[i, j] for i in range_client])
                     <= cl_ni[j] * fac_vars[j]
                 )
         else:
@@ -529,34 +521,34 @@ def add_facility_capacity_constraint(
 
     @staticmethod
     def add_client_demand_satisfaction_constraint(
-        obj: T_FacModel, model, ni, range_client, range_facility
+        obj: T_FacModel, range_client, range_facility
     ) -> None:
         """
 
         Parameters
         ----------
-        obj: T_FacModel
+        obj : T_FacModel
             bounded type of LocateSolver class
-        model: pulp.LpProblem
-            optimization model problem
-        ni: np.array
-            two-dimensional array that defines candidate sites between facility points within a distance to supply {i} demand point
-        range_client: range
+        range_client : range
             range of demand points quantity
-        range_facility: range
+        range_facility : range
             range of facility points quantity
 
         Returns
         -------
+
         None
+
         """
         if hasattr(obj, "fac_vars") and hasattr(obj, "cli_assgn_vars"):
             cli_assn_vars = getattr(obj, "cli_assgn_vars")
+            model = getattr(obj, "problem")
+            aij = getattr(obj, "aij")
 
             for i in range_client:
                 model += (
                     pulp.lpSum(
-                        [int(ni[i][j]) * cli_assn_vars[i][j] for j in range_facility]
+                        [int(aij[i, j]) * cli_assn_vars[i, j] for j in range_facility]
                     )
                     == 1
                 )
@@ -568,8 +560,6 @@ def add_client_demand_satisfaction_constraint(
     @staticmethod
     def add_maximal_coverage_constraint(
         obj: T_FacModel,
-        model: pulp.LpProblem,
-        ni: np.array,
         range_facility: range,
         range_client: range,
     ) -> None:
@@ -581,16 +571,11 @@ def add_maximal_coverage_constraint(
         Parameters
         ----------
 
-        obj: T_FacModel
+        obj : T_FacModel
             A bounded type of the ``LocateSolver`` class.
-        model: pulp.LpProblem
-            A ``pulp`` instance of an optimization model.
-        ni: numpy.array
-            A 2D array that defines candidate sites between facility
-            points within a distance to supply demand points.
-        range_facility: range
+        range_facility : range
             The range of facility points.
-        range_client: range
+        range_client : range
             The range of demand points.
 
         Returns
@@ -602,6 +587,8 @@ def add_maximal_coverage_constraint(
         if hasattr(obj, "fac_vars") and hasattr(obj, "cli_vars"):
             fac_vars = getattr(obj, "fac_vars")
             dem_vars = getattr(obj, "cli_vars")
+            model = getattr(obj, "problem")
+            ni = getattr(obj, "aij")
 
             for i in range_client:
                 model += (
@@ -617,7 +604,6 @@ def add_maximal_coverage_constraint(
     @staticmethod
     def add_assignment_constraint(
         obj: T_FacModel,
-        model: pulp.LpProblem,
         range_facility: range,
         range_client: range,
     ) -> None:
@@ -629,13 +615,11 @@ def add_assignment_constraint(
         Parameters
         ----------
 
-        obj: T_FacModel
+        obj : T_FacModel
             A bounded type of the ``LocateSolver`` class.
-        model: pulp.LpProblem
-            A ``pulp`` instance of an optimization model.
-        range_facility: range
+        range_facility : range
             The range of facility points.
-        range_client: range
+        range_client : range
             The range of demand points.
 
         Returns
@@ -646,9 +630,10 @@ def add_assignment_constraint(
         """
         if hasattr(obj, "cli_assgn_vars"):
             cli_assgn_vars = getattr(obj, "cli_assgn_vars")
+            model = getattr(obj, "problem")
 
             for i in range_client:
-                model += pulp.lpSum([cli_assgn_vars[i][j] for j in range_facility]) == 1
+                model += pulp.lpSum([cli_assgn_vars[i, j] for j in range_facility]) == 1
         else:
             raise AttributeError(
                 "Before setting assignment constraints "
@@ -658,7 +643,6 @@ def add_assignment_constraint(
     @staticmethod
     def add_opening_constraint(
         obj: T_FacModel,
-        model: pulp.LpProblem,
         range_facility: range,
         range_client: range,
     ) -> None:
@@ -670,13 +654,11 @@ def add_opening_constraint(
         Parameters
         ----------
 
-        obj: T_FacModel
+        obj : T_FacModel
             A bounded type of the ``LocateSolver`` class.
-        model: pulp.LpProblem
-            A ``pulp`` instance of an optimization model.
-        range_facility: range
+        range_facility : range
             The range of facility points.
-        range_client: range
+        range_client : range
             The range of demand points.
 
         Returns
@@ -688,10 +670,11 @@ def add_opening_constraint(
         if hasattr(obj, "cli_assgn_vars"):
             cli_assgn_vars = getattr(obj, "cli_assgn_vars")
             fac_vars = getattr(obj, "fac_vars")
+            model = getattr(obj, "problem")
 
             for i in range_client:
                 for j in range_facility:
-                    model += fac_vars[j] - cli_assgn_vars[i][j] >= 0
+                    model += fac_vars[j] - cli_assgn_vars[i, j] >= 0
         else:
             raise AttributeError(
                 "Before setting opening constraints "
@@ -701,7 +684,6 @@ def add_opening_constraint(
     @staticmethod
     def add_minimized_maximum_constraint(
         obj: T_FacModel,
-        model: pulp.LpProblem,
         cost_matrix: np.array,
         range_facility: range,
         range_client: range,
@@ -714,15 +696,13 @@ def add_minimized_maximum_constraint(
         Parameters
         ----------
 
-        obj: T_FacModel
+        obj : T_FacModel
             A bounded type of the ``LocateSolver`` class.
-        model: pulp.LpProblem
-            A ``pulp`` instance of an optimization model.
         cost_matrix : numpy.array
             A cost matrix in the form of a 2D array between origins and destinations.
-        range_facility: range
+        range_facility : range
             The range of facility points.
-        range_client: range
+        range_client : range
             The range of demand points.
 
         Returns
@@ -740,12 +720,13 @@ def add_minimized_maximum_constraint(
         if hasattr(obj, "cli_assgn_vars") and hasattr(obj, "weight_var"):
             cli_assgn_vars = getattr(obj, "cli_assgn_vars")
             weight_var = getattr(obj, "weight_var")
+            model = getattr(obj, "problem")
 
             for i in range_client:
                 model += (
                     pulp.lpSum(
                         [
-                            cli_assgn_vars[i][j] * cost_matrix[i][j]
+                            cli_assgn_vars[i, j] * cost_matrix[i, j]
                             for j in range_facility
                         ]
                     )
@@ -760,7 +741,6 @@ def add_minimized_maximum_constraint(
     @staticmethod
     def add_p_dispersion_interfacility_constraint(
         obj: T_FacModel,
-        model: pulp.LpProblem,
         cost_matrix: np.array,
         range_facility: range,
     ) -> None:
@@ -772,13 +752,11 @@ def add_p_dispersion_interfacility_constraint(
         Parameters
         ----------
 
-        obj: T_FacModel
+        obj : T_FacModel
             A bounded type of the ``LocateSolver`` class.
-        model: pulp.LpProblem
-            A ``pulp`` instance of an optimization model.
         cost_matrix : numpy.array
             A cost matrix in the form of a 2D array between all facility points.
-        range_facility: range
+        range_facility : range
             The range of facility points.
 
         Returns
@@ -789,13 +767,14 @@ def add_p_dispersion_interfacility_constraint(
         """
         if hasattr(obj, "disperse_var") and hasattr(obj, "fac_vars"):
             M = cost_matrix.max()
+            model = getattr(obj, "problem")
 
             for i in range_facility:
                 for j in range_facility:
                     if j <= i:
                         continue
                     else:
-                        dij = cost_matrix[i][j]
+                        dij = cost_matrix[i, j]
                         model += (
                             pulp.lpSum(
                                 [(dij + M * (2 - obj.fac_vars[i] - obj.fac_vars[j]))]
diff --git a/spopt/locate/coverage.py b/spopt/locate/coverage.py
index 87007677..832fb51b 100644
--- a/spopt/locate/coverage.py
+++ b/spopt/locate/coverage.py
@@ -181,8 +181,9 @@ def from_cost_matrix(
 
         """
 
+        n_cli = cost_matrix.shape[0]
+        r_cli = range(n_cli)
         r_fac = range(cost_matrix.shape[1])
-        r_cli = range(cost_matrix.shape[0])
 
         model = pulp.LpProblem(name, pulp.LpMinimize)
         lscp = LSCP(name, model)
@@ -198,13 +199,13 @@ def from_cost_matrix(
         lscp.aij[cost_matrix <= service_radius] = 1
 
         if (demand_quantity_arr is None) and (facility_capacity_arr is not None):
-            demand_quantity_arr = np.ones(cost_matrix.shape[0])
+            demand_quantity_arr = np.ones(n_cli)
 
         FacilityModelBuilder.add_facility_integer_variable(lscp, r_fac, "y[{i}]")
 
         if predefined_facilities_arr is not None:
             FacilityModelBuilder.add_predefined_facility_constraint(
-                lscp, lscp.problem, predefined_facilities_arr
+                lscp, predefined_facilities_arr
             )
 
         if demand_quantity_arr is not None:
@@ -222,22 +223,15 @@ def from_cost_matrix(
             )
 
             FacilityModelBuilder.add_facility_capacity_constraint(
-                lscp,
-                lscp.problem,
-                facility_capacity_arr,
-                demand_quantity_arr,
-                r_fac,
-                r_cli,
+                lscp, demand_quantity_arr, facility_capacity_arr, r_cli, r_fac
             )
 
             FacilityModelBuilder.add_client_demand_satisfaction_constraint(
-                lscp, lscp.problem, lscp.aij, r_cli, r_fac
+                lscp, r_cli, r_fac
             )
 
         else:
-            FacilityModelBuilder.add_set_covering_constraint(
-                lscp, lscp.problem, lscp.aij, r_fac, r_cli
-            )
+            FacilityModelBuilder.add_set_covering_constraint(lscp, r_cli, r_fac)
 
         lscp.__add_obj()
 
@@ -432,7 +426,7 @@ def facility_client_array(self) -> None:
             array_cli = []
             if fac_vars[j].value() > 0:
                 for i in range(self.aij.shape[0]):
-                    if self.aij[i][j] > 0:
+                    if self.aij[i, j] > 0:
                         array_cli.append(i)
 
             self.fac2cli.append(array_cli)
@@ -666,8 +660,8 @@ def from_cost_matrix(
             lscp = LSCP.from_cost_matrix(cost_matrix, service_radius)
         lscp.solve(solver)
 
-        r_fac = range(cost_matrix.shape[1])
         r_cli = range(cost_matrix.shape[0])
+        r_fac = range(cost_matrix.shape[1])
 
         model = pulp.LpProblem(name, pulp.LpMaximize)
 
@@ -682,16 +676,12 @@ def from_cost_matrix(
 
         if predefined_facilities_arr is not None:
             FacilityModelBuilder.add_predefined_facility_constraint(
-                lscpb, lscpb.problem, predefined_facilities_arr
+                lscpb, predefined_facilities_arr
             )
 
         lscpb.__add_obj()
-        FacilityModelBuilder.add_facility_constraint(
-            lscpb, lscpb.problem, lscpb.lscp_obj_value
-        )
-        FacilityModelBuilder.add_backup_covering_constraint(
-            lscpb, lscpb.problem, lscpb.aij, r_fac, r_cli
-        )
+        FacilityModelBuilder.add_facility_constraint(lscpb, lscpb.lscp_obj_value)
+        FacilityModelBuilder.add_backup_covering_constraint(lscpb, r_fac, r_cli)
 
         return lscpb
 
@@ -874,7 +864,7 @@ def facility_client_array(self) -> None:
             array_cli = []
             if fac_vars[j].value() > 0:
                 for i in range(self.aij.shape[0]):
-                    if self.aij[i][j] > 0:
+                    if self.aij[i, j] > 0:
                         array_cli.append(i)
 
             self.fac2cli.append(array_cli)
@@ -1101,8 +1091,10 @@ def from_cost_matrix(
         99.0
 
         """
+
+        n_cli = cost_matrix.shape[0]
+        r_cli = range(n_cli)
         r_fac = range(cost_matrix.shape[1])
-        r_cli = range(cost_matrix.shape[0])
 
         model = pulp.LpProblem(name, pulp.LpMaximize)
         mclp = MCLP(name, model)
@@ -1112,20 +1104,18 @@ def from_cost_matrix(
 
         mclp.aij = np.zeros(cost_matrix.shape)
         mclp.aij[cost_matrix <= service_radius] = 1
-        weights = np.reshape(weights, (cost_matrix.shape[0], 1))
+        weights = np.reshape(weights, (n_cli, 1))
 
         mclp.__add_obj(weights, r_cli)
 
         if predefined_facilities_arr is not None:
             FacilityModelBuilder.add_predefined_facility_constraint(
-                mclp, mclp.problem, predefined_facilities_arr
+                mclp, predefined_facilities_arr
             )
 
-        FacilityModelBuilder.add_maximal_coverage_constraint(
-            mclp, mclp.problem, mclp.aij, r_fac, r_cli
-        )
+        FacilityModelBuilder.add_maximal_coverage_constraint(mclp, r_fac, r_cli)
 
-        FacilityModelBuilder.add_facility_constraint(mclp, mclp.problem, p_facilities)
+        FacilityModelBuilder.add_facility_constraint(mclp, p_facilities)
 
         return mclp
 
@@ -1327,7 +1317,7 @@ def facility_client_array(self) -> None:
             if fac_vars[j].value() > 0:
                 for i in range(self.aij.shape[0]):
                     if cli_vars[i].value() > 0:
-                        if self.aij[i][j] > 0:
+                        if self.aij[i, j] > 0:
                             array_cli.append(i)
 
             self.fac2cli.append(array_cli)
diff --git a/spopt/locate/p_center.py b/spopt/locate/p_center.py
index d1ab3fcc..6b657c9c 100644
--- a/spopt/locate/p_center.py
+++ b/spopt/locate/p_center.py
@@ -200,22 +200,16 @@ def from_cost_matrix(
 
         if predefined_facilities_arr is not None:
             FacilityModelBuilder.add_predefined_facility_constraint(
-                p_center, p_center.problem, predefined_facilities_arr
+                p_center, predefined_facilities_arr
             )
 
         p_center.__add_obj()
 
-        FacilityModelBuilder.add_facility_constraint(
-            p_center, p_center.problem, p_facilities
-        )
-        FacilityModelBuilder.add_assignment_constraint(
-            p_center, p_center.problem, r_fac, r_cli
-        )
-        FacilityModelBuilder.add_opening_constraint(
-            p_center, p_center.problem, r_fac, r_cli
-        )
+        FacilityModelBuilder.add_facility_constraint(p_center, p_facilities)
+        FacilityModelBuilder.add_assignment_constraint(p_center, r_fac, r_cli)
+        FacilityModelBuilder.add_opening_constraint(p_center, r_fac, r_cli)
         FacilityModelBuilder.add_minimized_maximum_constraint(
-            p_center, p_center.problem, cost_matrix, r_fac, r_cli
+            p_center, cost_matrix, r_fac, r_cli
         )
 
         return p_center
@@ -388,7 +382,7 @@ def facility_client_array(self) -> None:
             array_cli = []
             if fac_vars[j].value() > 0:
                 for i in range(len(cli_vars)):
-                    if cli_vars[i][j].value() > 0:
+                    if cli_vars[i, j].value() > 0:
                         array_cli.append(i)
 
             self.fac2cli.append(array_cli)
diff --git a/spopt/locate/p_dispersion.py b/spopt/locate/p_dispersion.py
index 2efc2b21..fc3c0c5d 100644
--- a/spopt/locate/p_dispersion.py
+++ b/spopt/locate/p_dispersion.py
@@ -88,13 +88,13 @@ def from_cost_matrix(
         Parameters
         ----------
 
-        cost_matrix: np.array
+        cost_matrix : np.array
             A cost matrix in the form of a 2D array between origins and destinations.
         p_facilities : int
             The number of facilities to be located.
         predefined_facilities_arr : numpy.array (default None)
             Predefined facilities that must appear in the solution.
-        name: str (default 'P-Dispersion')
+        name : str (default 'P-Dispersion')
             The name of the problem.
 
         Returns
@@ -173,16 +173,16 @@ def from_cost_matrix(
         )
 
         FacilityModelBuilder.add_facility_constraint(
-            p_dispersion, p_dispersion.problem, p_dispersion.p_facilities
+            p_dispersion, p_dispersion.p_facilities
         )
 
         if predefined_facilities_arr is not None:
             FacilityModelBuilder.add_predefined_facility_constraint(
-                p_dispersion, p_dispersion.problem, predefined_facilities_arr
+                p_dispersion, predefined_facilities_arr
             )
 
         FacilityModelBuilder.add_p_dispersion_interfacility_constraint(
-            p_dispersion, p_dispersion.problem, cost_matrix, r_fac
+            p_dispersion, cost_matrix, r_fac
         )
 
         return p_dispersion
@@ -208,14 +208,14 @@ def from_geodataframe(
             Facility locations.
         facility_col : str
             Facility candidate sites geometry column name.
-        p_facilities: int
+        p_facilities : int
            The number of facilities to be located.
         predefined_facility_col : str (default None)
             Column name representing facilities are already defined.
         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>`_.
-        name: str (default 'P-Dispersion')
+        name : str (default 'P-Dispersion')
             The name of the problem.
 
         Returns
diff --git a/spopt/locate/p_median.py b/spopt/locate/p_median.py
index 01ad49d9..d8a30a5c 100644
--- a/spopt/locate/p_median.py
+++ b/spopt/locate/p_median.py
@@ -99,9 +99,9 @@ def __add_obj(self, range_clients: range, range_facility: range) -> None:
         Parameters
         ----------
 
-        range_clients: range
+        range_clients : range
             The range of demand points.
-        range_facility: range
+        range_facility : range
             The range of facility point.
 
         Returns
@@ -115,7 +115,7 @@ def __add_obj(self, range_clients: range, range_facility: range) -> None:
         self.problem += (
             pulp.lpSum(
                 [
-                    self.aij[i][j] * cli_assgn_vars[i][j]
+                    self.aij[i, j] * cli_assgn_vars[i, j]
                     for i in range_clients
                     for j in range_facility
                 ]
@@ -138,7 +138,7 @@ def from_cost_matrix(
         Parameters
         ----------
 
-        cost_matrix: numpy.array
+        cost_matrix : numpy.array
             A cost matrix in the form of a 2D array between origins and destinations.
         weights : numpy.array
             A 1D array of service load or population demand.
@@ -229,13 +229,14 @@ def from_cost_matrix(
         3.027
 
         """
-        r_cli = range(cost_matrix.shape[0])
+        n_cli = cost_matrix.shape[0]
+        r_cli = range(n_cli)
         r_fac = range(cost_matrix.shape[1])
 
         model = pulp.LpProblem(name, pulp.LpMinimize)
 
         weights_sum = weights.sum()
-        weights = np.reshape(weights, (cost_matrix.shape[0], 1))
+        weights = np.reshape(weights, (n_cli, 1))
         aij = weights * cost_matrix
 
         p_median = PMedian(name, model, aij, weights_sum)
@@ -247,20 +248,14 @@ def from_cost_matrix(
 
         if predefined_facilities_arr is not None:
             FacilityModelBuilder.add_predefined_facility_constraint(
-                p_median, p_median.problem, predefined_facilities_arr
+                p_median, predefined_facilities_arr
             )
 
         p_median.__add_obj(r_cli, r_fac)
 
-        FacilityModelBuilder.add_facility_constraint(
-            p_median, p_median.problem, p_facilities
-        )
-        FacilityModelBuilder.add_assignment_constraint(
-            p_median, p_median.problem, r_fac, r_cli
-        )
-        FacilityModelBuilder.add_opening_constraint(
-            p_median, p_median.problem, r_fac, r_cli
-        )
+        FacilityModelBuilder.add_facility_constraint(p_median, p_facilities)
+        FacilityModelBuilder.add_assignment_constraint(p_median, r_fac, r_cli)
+        FacilityModelBuilder.add_opening_constraint(p_median, r_fac, r_cli)
 
         return p_median
 
@@ -443,7 +438,7 @@ def facility_client_array(self) -> None:
             array_cli = []
             if fac_vars[j].value() > 0:
                 for i in range(len(cli_vars)):
-                    if cli_vars[i][j].value() > 0:
+                    if cli_vars[i, j].value() > 0:
                         array_cli.append(i)
 
             self.fac2cli.append(array_cli)
diff --git a/spopt/tests/test_locate.py b/spopt/tests/test_locate.py
index 8938f85b..01ffe1b1 100644
--- a/spopt/tests/test_locate.py
+++ b/spopt/tests/test_locate.py
@@ -869,58 +869,56 @@ def setup_method(self) -> None:
     def test_attribute_error_add_set_covering_constraint(self):
         with pytest.raises(AttributeError, match="Before setting coverage constraints"):
             dummy_class = LSCP("dummy", pulp.LpProblem("name"))
-            dummy_matrix = numpy.array([])
+            dummy_class.aij = numpy.array([])
             dummy_range = range(1)
             FacilityModelBuilder.add_set_covering_constraint(
-                dummy_class, dummy_class.problem, dummy_matrix, dummy_range, dummy_range
+                dummy_class, dummy_range, dummy_range
             )
 
     def test_attribute_error_add_facility_constraint(self):
         with pytest.raises(AttributeError, match="Before setting facility constraint"):
-            dummy_class = LSCP("dummy", pulp.LpProblem("name"))
+            dummy_class = MCLP("dummy", pulp.LpProblem("name"))
             dummy_p_facility = 1
-            FacilityModelBuilder.add_facility_constraint(
-                dummy_class, dummy_class.problem, 1
-            )
+            FacilityModelBuilder.add_facility_constraint(dummy_class, 1)
 
     def test_attribute_error_add_maximal_coverage_constraint(self):
         with pytest.raises(
             AttributeError, match="Before setting maximal coverage constraints"
         ):
-            dummy_class = LSCP("dummy", pulp.LpProblem("name"))
-            dummy_matrix = numpy.array([])
+            dummy_class = MCLP("dummy", pulp.LpProblem("name"))
+            dummy_class.aij = numpy.array([])
             dummy_range = range(1)
             FacilityModelBuilder.add_maximal_coverage_constraint(
-                dummy_class, dummy_class.problem, dummy_matrix, dummy_range, dummy_range
+                dummy_class, dummy_range, dummy_range
             )
 
     def test_attribute_error_add_assignment_constraint(self):
         with pytest.raises(
             AttributeError, match="Before setting assignment constraints"
         ):
-            dummy_class = LSCP("dummy", pulp.LpProblem("name"))
+            dummy_class = PMedian("dummy", pulp.LpProblem("name"), numpy.array([]), 1)
             dummy_range = range(1)
             FacilityModelBuilder.add_assignment_constraint(
-                dummy_class, dummy_class.problem, dummy_range, dummy_range
+                dummy_class, dummy_range, dummy_range
             )
 
     def test_attribute_error_add_opening_constraint(self):
         with pytest.raises(AttributeError, match="Before setting opening constraints"):
-            dummy_class = LSCP("dummy", pulp.LpProblem("name"))
+            dummy_class = PMedian("dummy", pulp.LpProblem("name"), numpy.array([]), 1)
             dummy_range = range(1)
             FacilityModelBuilder.add_opening_constraint(
-                dummy_class, dummy_class.problem, dummy_range, dummy_range
+                dummy_class, dummy_range, dummy_range
             )
 
     def test_attribute_error_add_minimized_maximum_constraint(self):
         with pytest.raises(
             AttributeError, match="Before setting minimized maximum constraints"
         ):
-            dummy_class = LSCP("dummy", pulp.LpProblem("name"))
             dummy_matrix = numpy.array([])
+            dummy_class = PCenter("dummy", pulp.LpProblem("name"), dummy_matrix)
             dummy_range = range(1)
             FacilityModelBuilder.add_minimized_maximum_constraint(
-                dummy_class, dummy_class.problem, dummy_matrix, dummy_range, dummy_range
+                dummy_class, dummy_matrix, dummy_range, dummy_range
             )
 
     def test_error_lscp_different_crs(self):
diff --git a/spopt/tests/test_lscpb.py b/spopt/tests/test_lscpb.py
index 486e6064..93e70056 100644
--- a/spopt/tests/test_lscpb.py
+++ b/spopt/tests/test_lscpb.py
@@ -330,18 +330,15 @@ def test_warning_lscpb_demand_geodataframe(self):
                 pulp.PULP_CBC_CMD(msg=False),
             )
 
-    def test_attribute_error_add_facility_constraint(self):
+    def test_attribute_error_add_backup_covering_constraint(self):
         with pytest.raises(AttributeError, match="Before setting backup coverage"):
             dummy_class = LSCPB(
                 "dummy", pulp.LpProblem("name"), pulp.PULP_CBC_CMD(msg=False)
             )
-            dummy_p_facility = 1
             dummy_fac_r = 0
             dummy_cli_r = 0
             FacilityModelBuilder.add_backup_covering_constraint(
                 dummy_class,
-                dummy_class.problem,
-                dummy_p_facility,
                 dummy_fac_r,
                 dummy_cli_r,
             )
diff --git a/spopt/tests/test_p_dispersion.py b/spopt/tests/test_p_dispersion.py
index c4eefbf9..22e31327 100644
--- a/spopt/tests/test_p_dispersion.py
+++ b/spopt/tests/test_p_dispersion.py
@@ -192,11 +192,7 @@ def test_attribute_error_add_facility_constraint(self):
         with pytest.raises(AttributeError, match="Before setting facility constraint"):
             dummy_p_facility = 1
             dummy_class = PDispersion("dummy", pulp.LpProblem("name"), dummy_p_facility)
-            FacilityModelBuilder.add_facility_constraint(
-                dummy_class,
-                dummy_class.problem,
-                dummy_p_facility,
-            )
+            FacilityModelBuilder.add_facility_constraint(dummy_class, dummy_p_facility)
 
     def test_attribute_error_add_p_dispersion_interfacility_constraint(self):
         with pytest.raises(
@@ -208,7 +204,6 @@ def test_attribute_error_add_p_dispersion_interfacility_constraint(self):
             dummy_class = PDispersion("dummy", pulp.LpProblem("name"), dummy_p_facility)
             FacilityModelBuilder.add_p_dispersion_interfacility_constraint(
                 dummy_class,
-                dummy_class.problem,
                 dummy_matrix,
                 dummy_range,
             )
@@ -218,9 +213,7 @@ def test_attribute_error_add_predefined_facility_constraint(self):
             dummy_p_facility = 1
             dummy_matrix = numpy.array([])
             dummy_class = PDispersion("dummy", pulp.LpProblem("name"), dummy_p_facility)
-            FacilityModelBuilder.add_facility_constraint(
-                dummy_class, dummy_class.problem, dummy_matrix
-            )
+            FacilityModelBuilder.add_facility_constraint(dummy_class, dummy_matrix)
 
     def test_warning_facility_geodataframe(self):
         with pytest.warns(