key: cord-0124602-ftglxojz authors: Tian, Qingyun; Lin, Yun Hui; He, Dongdong title: Decomposition approach for Stackelberg P-median problem with user preferences date: 2021-09-18 journal: nan DOI: nan sha: d4b1e9a1344b37ebe982814bf947d0aa08035bad doc_id: 124602 cord_uid: ftglxojz The P-median facility location problem with user preferences (PUP) studies an operator that locates P facilities to serve customers/users in a cost-efficient manner, upon anticipating customer preferences and choices. The problem can be visualized as a leader-follower game in which the operator is the leader that opens facilities, whereas the customer is the follower who observes the operator's location decision at first and then seeks services from the most preferred facility. Such a modeling perspective is of practical importance as we have witnessed its applications to various problems, such as the establishment of power plants in energy markets and the location of healthcare service centers for COVID-19 Vaccination. Despite that a considerable number of solution methodologies have been proposed, many of them are heuristic methods whose solution quality cannot be easily verified. Moreover, due to the hardness of the problems, existing exact approaches have limited performance. Motivated by these observations, we aim to develop an efficient exact algorithm for solving large-scale PUP models. We first propose a branch-and-cut decomposition algorithm and then design accelerated techniques to further enhance the performance. Using a broad testbed, we show that our algorithm outperforms various exact approaches by a large margin, and the advantage can go up to several orders of magnitude in terms of computational time in some datasets. Finally, we conduct sensitivity analysis to draw additional implications and to highlight the importance of considering user preferences when they exist. The classical facility location problem (FLP) studies the strategic planning of an operator who provides services to customers. The problem involves the determination of a set of open facilities from candidate sites, aiming to minimize the cost of satisfying customer's demand/or maximize some objectives of interest (e.g., to maximize the profit of the operator or to cover as many customers as possible). To date, a large number of models have been proposed to address the FLPs under various settings. For comprehensive discussions of models, algorithms, and applications, we refer readers to Melo et al. (2009) and Daskin (2011) . In the FLP literature, it is commonly assumed that customers (users) are allocated to their closest facilities or the operator acts as a centralized decision-maker that assigns customers to facilities and serves them in a cost-efficient fashion. However, this allocation/assignment rule is relatively simple and does not account for customer preferences. In many scenarios, customers may not prefer the assigned facility. They may instead seek services from some other facility, from which it could be costly for the operator to serve them. This implies that there could be a mismatch between the preferences of the customers and the desirable customer allocation/assignment of the operator. Therefore, when locating the facilities, the operator needs to consider this mismatch, giving rise to a variant of the FLP that is often referred to as the facility location problem with user preferences (FLPUP) or with order (Camacho-Vallejo et al., 2014b; Hansen et al., 2004; Marić et al., 2012; Vasil'ev et al., 2009; Vasilyev and Klimentova, 2010) . When the number of facilities to be set up is fixed, the problem is called the P-median problem with user preferences (PUP) Camacho-Vallejo et al., 2014a; Díaz et al., 2017) . In the context of FLPUP and PUP, customers are studied as individual decision-makers that select facilities to serve them based on their preferences. The resultant problem is visualized from a fully decentralized perspective and thus can be characterized as a Stackelberg leader-follower game, in which the operator is the leader and the customer is the follower. Figure 1 illustrates the decision process of both the operator and the customer. At the initial stage, the customer and the candidate facilities are known as discrete nodes. At Stage 1, the operator selects locations from the candidate sites. Customers, as followers, observe the operator's decisions at first and then request services from facilities at Stage 2. Finally, the operator proceeds to serve customers according to the requests. Operator's candidate facilities Operator's selected facilities Operator's services Customers' requests It is worth noting that when locating facilities at Stage 1, the operator needs to anticipate how customers choose facilities at Stage 2 because the decentralized customers want to be served by the most preferred facility that is not necessarily cost-efficient from the operator's side. Customers' requests could thus have a significant impact on the operators' service profile and the service provision cost at Stage 3. Since facilities are generally built for a long-term purpose and it will be costly and difficult to alter the decision, it is important for the operator to consider customer preferences and anticipate customers' choices. Such a modeling perspective is of practical importance. Recently, we have witnessed several applications of the FLPUP and PUP to real-world problems, such as the establishment of power plants in energy markets (Lotfi et al., 2021) and the location of healthcare service centers for COVID-19 Vaccination (Cabezas et al., 2021) . This paper investigates a PUP model. Due to the hidden bilevel structure, the resultant mathematical model is computationally challenging. Until now, we are unaware of any existing exact solution approach that can solve large-scale PUP instances within a reasonable computational time. This paper thus aims to develop an efficient exact algorithm for large-scale PUP instances. Below, we briefly review the related literature with a focus on the solution approaches to clearly position our contributions. Literature review. The first study that considered customer preferences is Hanjoul and Peeters (1987) where each customer has a preference ordering towards the set of potential locations, and the simple plant location problem was presented with a set of constraints on preference ordering. A greedy heuristic algorithm was proposed to conduct the numerical tests on multiple instances. Since then, there is a proliferation of research on the FLPUP and PUP. Although the problem has a Stackelbeg structure (as illustrated in Figure 1 ), single-level formulation indeed exists. A typical exact approach thus consists of formulating the FLPUP and PUP as a single-level mixed-integer linear programming (MILP) and then solving the resultant MILP using existing solvers (possibly strengthened by some valid inequalities). In general, the single-level formulation can be divided into two research streams. The first one is based on the primal-dual information of the lower-level problem and leverages the KKT conditions to recast the bilevel model into a mixed-integer program with complementarity constraints. Then the non-convex complementarity constraints are linearized, giving rise to a MILP that can be directly solved by off-the-shelf solvers (e.g., see Camacho-Vallejo et al. (2014a) ; Casas-Ramírez and Camacho-Vallejo (2017); Cao and Chen (2006) ; Lotfi et al. (2021) ). However, due to the use of the KKT conditions, such an approach needs to introduce a large number of variables (i.e., the dual variables) that dramatically increases the formulation size. As a consequence, its performance is typically limited. The second stream, essentially, utilizes the closest assignment constraints (CACs) to subtly rewrite the lowerlevel problem into a set of linear constraints without additional variables. For a comprehensive discussion on the CACs, we refer readers to Espejo et al. (2012) . To date, a considerable number of studies have adopted the CAC-based reformulation. For example, Cánovas et al. (2007) proposed a strengthened CAC to speed up solving the model of Hanjoul and Peeters (1987) . Vasil'ev et al. (2009) provided several valid inequalities when the bilevel problem was reformulated into a MILP using CACs and demonstrated the tightness of the inequalities in terms of their linear relaxations and integrality gaps. Then Vasilyev and Klimentova (2010) implemented a branch-and-cut method based on the proposed family of valid inequalities in Vasil'ev et al. (2009) . A recent application can be found in Casas-Ramírez and Camacho-Vallejo (2017) . We point out that the above approaches are exact and can yield the proven optimal solution, provided that the reformulated MILP model can be solved optimally. However, the MILP model itself is a challenging problem that cannot scale well on the problem size, rendering these approaches computational prohibitive for large-scale instances. Therefore, heuristic approaches are gaining popularity. Recent years have witnessed the use of Lagrangian-based heuristics to various problems involving customer preferences (e.g., see Cabezas and García (2018) ; Cabezas et al. (2021) ; Lee and Lee (2012) ). Moreover, there is an explosion of the related research on metaheuristics. Marić et al. (2012) proposed three metaheuristic methods for solving FLPUP, i.e., Particle Swarm Optimization (PSO), Simulated Annealing (SA), and a combination of Reduced and Basic Variable Neighborhood Search Method (RVNS-VNS). To solve a large-scale bilevel FLPUP, Camacho-Vallejo et al. (2014b) developed a population-based evolutionary algorithm (EA). A similar EA was applied to solve a bilevel facility location problem with cardinality constraints and preferences (Calvete et al., 2020) . Casas-Ramírez et al. (2018) obtained approximated solutions to the problem by using a cross-entropy method to search for improved location decisions. In Casas-Ramírez and Camacho-Vallejo (2017), the authors investigated the PUP model and designed a scatter search metaheuristic algorithm to efficiently handle large-scale problems. In the context of maximum coverage, Díaz et al. (2017) presented a bilevel model considering that customers will choose the facilities within a coverage radius. They proposed a greedy randomized adaptive search procedure (GRASP) heuristic and a hybrid GRASP-Tabu heuristic to find near-optimal solutions. Recently, Mrkela and Stanimirović (2021) investigated a similar problem as in Díaz et al. (2017) . Instead of locating P facilities, they considered a limited budget for locating facilities, and an efficient Variable Neighborhood Search (VNS) was proposed as the solution approach. Our contributions. Despite that there is a considerable number of research works dedicated to developing solution methodologies for the FLPUP and PUP, many of them are heuristic methods whose solution quality cannot be easily verified. Moreover, due to the hardness of the problems, existing exact approaches based on single-level reformulation techniques have limited performance. Such an observation motivates this paper to develop an efficient exact algorithm that can be applied to solve large-scale PUP models. In summary, our contributions are three-fold. (i)A branch-and-cut Benders decomposition algorithm. By exploring the problem structure, we prove that in a CAC-based MILP model (Casas-Ramírez and Camacho-Vallejo, 2017), the high-dimensional binary variables used to model customer preferences and facility choices can be relaxed to continuous variables. Using this condition, we propose a branch-and-cut algorithm, where we project out the preference-related variables and solve a master problem defined on a low-dimensional decision space. (ii) An enhanced Benders approach with analytical separation. The decomposition algorithm is expected to be more efficient than the MILP model since the branch-and-cut searching tree only involves the location variable of the original MILP and the problem size is thus substantially smaller. However, our numerical experiment reveals that such an algorithm cannot effectively expedite the solution process because generating the Benders cut (i.e., Benders separation) at nodes of the tree is not trivial, and significant computational time can be spent on the separation function. In light of this, we design an analytical approach for Benders separation to further speed up the algorithm. It turns out that our specialized separation is so efficient that we observe a significant performance improvement (with up to orders of magnitude in terms of computational time). Using a broad testbed, our computational experiments support that the enhanced Benders approach is able to handle large-scale instances satisfactorily. (iii) Additional managerial implications. We conduct sensitivity analysis using instances from a standard benchmark library. Our analysis reveals that when user preferences indeed exist (such as in the context of energy markets (Lotfi et al., 2021) and healthcare service (Cabezas et al., 2021) , the operator must consider these preferences and correctly anticipate user's choices; otherwise, the operator could suffer from a substantially higher cost, and opening more facilities could ironically result in additional service costs. The rest of the paper is organized as follows. We introduce the problem and the model in Section 2. We then develop the branch-and-cut Benders decomposition and the acceleration technique in Section 3. Extensive numerical studies are conducted in Section 4 to demonstrate the effectiveness of our algorithm and provide additional managerial implications. Finally, we conclude the paper in Section 5. This section describes the model formulation of the P-median problem with user preferences (PUP). Table 1 provides the main notations. : set of customers. J : set of candidate facilities. Parameters c ij : cost of serving customer i by facility j, ∀i ∈ I, j ∈ J g ij : disutility of facility j to customer i, ∀i ∈ I, j ∈ J π ij : normalized disutility of facility j to customer i, defined as g ij /(max j∈J g ij ), ∀i ∈ I, j ∈ J P : the number of open facilities Decision Variables x j : binary. If facility j is open; 0, otherwise, ∀j ∈ J y ij : binary. If customer i seeks the service from facility j; 0, otherwise, ∀i ∈ I, j ∈ J The virtual decision process of [PUP] is as explained in Figure 1 . We use J to denote the set of candidate facilities and I to denote the set of customers. At Stage 1, the operator will open P facilities, selecting from set J. We define a binary variable x j , ∀j ∈ J, which is 1 if facility j is open; 0, otherwise. After P facilities are built, customers then determine which facilities to seek services from at Stage 2. To reflect customers' choices, we define a binary variable y ij such that y ij = 1 if customer i requests services from facility j, ∀i ∈ I, j ∈ J. Following the standard assumption in the FLPUP and the PUP literature, we assume that facility j has a disutility to customer i, i.e., g ij , and customer i prefers facility j over facility k, if g ij < g ik . In other words, when both facility j and facility k are open, customer i will not select facility k if g ij < g ik . With these definitions, given a location decision x, customer preferences to the facilities can be modeled by which describes the decision problem at Stage 2. The objective function imposes that the customer will minimize the disutility. This is equivalent to state that he/she will request the service from the most preferred facility located by the operator. Given the solution y * from the above problem, the operator incurs a service provision cost of φ = i∈I j∈J c ij y * ij at Stage 3, where c ij is the cost of serving customer i from facility j (possibly weighted by the demand size or the relative "importance" of the customer). Note that the operator also wants to minimize the total service cost; therefore, when locating facilities at Stage 1, the operator needs to solve the following problem where Constraint (2d) states that y * is obtained by solving Problem (1). Therefore, [PUP] is indeed a mixed-integer bilevel program. Remark 1. Through this paper, we assume that customer preferences for facilities are different, i.e., elements in the vector g i· are distinct (g ij = g ik if j = k). Under this assumption, an operator's decision will only lead to a unique solution of y (i.e., Problem 1 has an unique solution), allowing us to avoid the discussion on optimistic and pessimistic strategies of the follower in the context of bilevel optimization. The bilevel structure of [PUP] renders the problem computationally challenging. Here, we present a single-level reformulation model of [PUP] to circumvent the difficulties of handling Problem (1) in [PUP] . In fact, minimizing the objective function (1a) in Problem (1) can be represented by the following inequality (3) imposes that, among all the possible allocations, customer i will be allocated to the facility with the minimum value of g i· . Essentially, this equation belongs to a type of closest assignment constraints (CACs). It was proposed by Berman et al. (2009) and has been used by Casas-Ramírez and Camacho-Vallejo (2017) in a PUP model. Now, define π ij = g ij /G i , i.e., the vector g i· is normalized by the large value in it to get π i· . We have the following single-level reformulation model for [PUP] min i∈I j∈J where Ω is the decision space for the operator This set is defined to simplify the notation for our later discussions. Now, [SRM] is a mixed-integer linear program (MILP) and is ready to be solved by modern MILP solvers. Remark 2. We point out that other types of CACs that can possibly be used to reformulate [PUP] into a single-level model have been discussed in detail in Espejo et al. (2012) . Among all CACs, we choose (4d) since it only involves O(|I| · |J|) constraints, whereas others may contain up to O(|I| · |J| · |J|) constraints. Remark 3. Besides the CAC-based reformulation, the primal-dual reformulation has also been discussed in the literature Casas-Ramírez et al., 2018; Lotfi et al., 2021) . The fundamental idea is to replace Problem (1) with its KKT conditions, thereby recasting [PUP] into a single-level primal-dual reformulation model [PDRM] . We provide the [PDRM] formulationl in Appendix A. We note that this approach introduces a large number of new variables into the formulation (i.e., there are |I| · (2|J| + 1) dual variables in the KKT conditions) and thus significantly increases the problem size. Nevertheless, the single-level reformulation approach leveraging the KKT conditions of the lower-level problem is probably the most widely used approach for solving bilevel programs; therefore, we will keep [PDRM] in this paper and use it as one of the benchmarks to our proposed algorithms. This section presents the solution methodology for [SRM] . In the literature, the single-level MILP model is typically solved using off-the-shelf solvers (e.g., CPLEX and Gurobi). As shown in Casas-Ramírez and Camacho-Vallejo (2017), such a straightforward approach cannot handle large-scale problems. Therefore, we develop a branch-and-cut Benders decomposition algorithm to expedite solving [SRM] , where Benders separation (i.e., the procedure of generating Benders cuts) at integer nodes of the searching tree is performed leveraging external solvers. However, we observe that such a standard Benders approach is still not efficient enough. Therefore, we further derive an acceleration technique to enhance the algorithm performance. Our decomposition approaches are inspired by the following observation. Lemma 1. In the [SRM] formulation, Constraint (4e) can be relaxed to y ij ≥ 0 without changing the solution and the objective. Proof. Suppose a location decision is made atx. Letτ be the set of open facilities, i.e.,τ = {∀j ∈ J |x j = 1}. Based on the condition y ij ≤x j , Constraint (4d) can be restated as j∈τ π ij y ij ≤ min j∈τ π ij , ∀i ∈ I. For notation simplicity, we drop subscript i in this proof. Let m denote the most preferred open facility for customer, i.e., m = arg min j∈τ π j . We have j∈τ π j y j ≤ π m . Note that π m < π j , ∀j ∈τ \ {m}, and j∈τ y j = 1. Therefore, j∈τ π j y j ≤ π m holds only if y m = 1 and y j = 0, ∀j ∈τ \ {m}, which also uses the condition that the elements in π i· (g i· ) are distinct. This indicates that relaxing y j ∈ {0, 1} to y j ≥ 0 will not change the solution since y j in [SRM] will automatically be an integer. This lemma implies that when the operator's decision is made, the remaining subproblem of [SRM] will become a linear program (LP) that prosesses strong duality. Therefore, it is possible to design a decomposition algorithm leveraging the dual information of the subproblem, which is exactly the idea of Benders decomposition. Modern Benders decomposition algorithms are typically implemented within the solver's branch-and-cut framework, and we have observed a considerable number of research works successfully employing such a technique to address facility location problems (Cordeau et al., 2019; Fischetti et al., 2016; Ljubić et al., 2012; Taherkhani et al., 2020) . In this paper, we design a branch-and-cut Benders decomposition to solve [SRM] . Our approach only maintains the location decision x in a master problem by projecting out the continuous variable y that is used to reflect the preferences and compute the service cost. Specifically, the master problem is defined as where Φ i (x) is obtained by solving the following subproblem Note that Φ i (x) is a convex function over x; therefore, [MP] is a mixed-integer convex optimization problem, and we can solve it by approximating Φ i (x) with its linear under-estimators, i.e., the Benders cuts. Given an integer solutionx from [MP], the Benders cut can be obtained by solving the dual problem of [SP i (x) ]. Let µ i and λ ij and v ij in (7) be the dual variables associated with the constraints. The dual subproblem is given by Let (λ,v,μ) be the optimal solution to the above LP. The following Benders cut defined atx arises. We can now solve [MP] using the following MILP formulation, i.e., the relaxed master problem, where Ξ is the set of dual variables. We can expand the set, each time we have an integer solution x and solve the corresponding [DSP i (x)]. These dual variables define Benders cuts that cut off non-optimal solutions and lead us to the proven optimal solution. We then illustrate our implementation of Benders decomposition. Modern advanced solvers, such as CPLEX and Gurobi, use branch-and-cut algorithms to solve (mixed-)integer programs. Apart from the general-purpose cuts that are embedded within the solvers, we can design a problemspecific separation function that uses a [rMP] solutionx as an input and generates violated cuts, which are inserted into the searching tree to improve the relaxation bound and/or cut off nonoptimal solutions. We perform Benders separation only at integer nodes because this is sufficient to guarantee the optimal convergence. More specifically, we use Gurobi 9.1.2 and rely on its default branching rules and embedded cutting planes. We initialize the master problem [rMP] without Benders cuts (i.e., Ξ = ∅). Now, when Gurobi visits an integer node, we retrieve the integer value ofx and generate the corresponding Benders cuts as in (9) by solving [DSP i (x)]. The Benders separation is performed within the callback function of Gurobi. Cuts are added as lazy cuts to the current node relaxation only if they violate the current solution (x,w) by more than 10 −5 of the relative violation, defined as the absolute violation of the cut divided by the current w i value. Due to the convexity of Φ i , these inserted cuts are globally valid and will eventually lead us to the optimal solution of [SRM] with a zero MIP gap. The above branch-and-cut Benders decomposition generates Benders cuts by solving the dual subproblem [DSP i (x)] using external solvers, each time an integer node is visited during the searching tree. However, [DSP i (x)] itself is a large-scale LP. For some problems, solving [DSP i (x)] is not trivial and may take up to several seconds. Noting that we typically need to visit a large number of integer nodes before the MIP gap vanishes, significant computational time could thus be consumed on Benders separation, causing a bottleneck in the branch-and-cut algorithm. In light of this, we design a specialized analytical approach for Benders separation, which is very fast and can, in general, expedite the algorithm by a large margin as demonstrated in our experiment. One should keep in mind that the procedure explained below relies on the condition that the solutionx is binary. Step 1: Solve [SP i (x)]. To start with, we point out that given a location decisionx, [SP i (x)] can be solved using a simple sorting algorithm since customers will seek the service from In other words, givenx, customer i will select facility m, and thus the optimal solution of y is Then, the leader incurs a cost of j∈J c ijȳij = c im to serve customer i, meaning that the optimal value of [SP i (x)] is Upon here, we have solved [SP i (x)] and obtained the objective c im and the solutionȳ. Step 2: Transform [DSP i (x)]. To solve [DSP i (x)], we first transform it into a problem of finding feasible solutions to a set of linear equations. Based on the value ofx j andȳ ij , there exist three cases. Case 1:x j = 1 andȳ ij = 0. In [SP i (x)], constraint (7c) is inactive; therefore, λ ij = 0. Moreover, constraint (7d) becomes π im ≤ π ij , which holds since m = arg min j∈τ π ij . This indicates that constraint (7d) is also inactive, leading to v ij = 0. Case 2: x j = 0 and y ij = 0. We have v ij = 0 since constraint (7d) is inactive. Case 3: x j = 1 and y ij = 1. We have j = m. Combining these cases, we conclude that j∈J v ij = v im since if j = m, then v ij = 0. Therefore, the objective function underx becomes The last equality holds since when j ∈τ , λ ij can be non-zero only if j = m according to Case 1. Note that, for customer i, Case 1 reflects the dual information of the set of open facilities excluding the most referred one, i.e.,τ \ {m}, whereas Case 2 reflects the dual information of the set of facilities are not open, i.e, J \τ . We can therefore rewritten [DSP i (x)] as Note that the optimal value is c im according to (14). We have c im = µ i − λ im − π im v im based on the strong duality of LP. Therefore, solving Problem (16) is equivalent to find out a feasible solution of the following linear equations where we replace µ i with c im + λ im + π im v im . Step 3: Obtain the dual variables. Now, we proceed to deriving an analytical feasible solution to Problem (17). Note that we can safely set λ im = 0 since to satisfy equations (17a) and (17b), λ im should be small enough. After some algebra, (17a) becomes Since where [z] + = max{z, 0}, and compute µ i by which utilizes the condition c im = µ i − λ im − π im v im and λ im is 0. Furthermore, based on (17b), we have Since λ ij ≥ 0, we set Note that λ ij = 0 whenx j = 1. We can express λ as Finally, we summarize the above three steps compactly into the following lemma. Lemma 2. Given an integer solutionx from [rMP], letτ = {∀j ∈ J |x j = 1} and m = arg min j∈τ π ij . An optimal dual variable (v, µ, λ) can be obtained analytically as where [z] + = max{z, 0}. To conclude the above discussion, given an integer solutionx, the optimal dual variables can be obtained by the procedure described in Lemma 2, which only involves a sorting algorithm (i.e., determining the most preferred facility) and some elemental matrix manipulation. Therefore, the analytical separation is very efficient. We emphasize that the integrality ofx is important for the analytical separation because Lemma 2 holds only ifx is an integer point. Therefore, similar to the standard branch-and-cut Benders decomposition, when we implement the accelerated Benders approach, analytical separation of Benders cuts is performed only at integer nodes. Moreover, following Lin and Tian (2021) , we conservatively set the integer feasibility tolerance and primal feasibility tolerance to the minimum values provided by Gurobi, i.e., the Gurobi parameters Int-FeasTol and FeasibilityTol are set to 10 −9 to enhance the numerical stability, despite doing so may potentially increase the overall computational time. Then the current "integer"x from the node is rounded to ensure the integrality. This section presents numerical experiments. We first conduct extensive computational experiments using a broad testbed to demonstrate the computational efficiency of the proposed Benders-based approaches. Then we discuss the importance of considering user preferences in facility location problems using instances from a standard library. We start with the computational experiment. Throughout this section, we refer to the standard Branch-and-cut Benders decomposition in Section 3.1 as Benders and the accelerated Benders approach with analytical separation in Section 3.2 as Benders-AS. Furthermore, we test the performance of two direct approaches where Gurobi is used to solve [SRM] presented in Section 2 and [PDRM] presented in Appendix A. They serve as benchmarks and are referred to as SRM and PDRM, respectively. All experiments are done on a 16 GB memory macOS computer with a 2.6 GHz Intel Core i7 processor, using Gurobi 9.1.2 as the solver and Python as the programming language. Moreover, the time limit is set to 7200 seconds. Our testbed consists of three datasets with different scales and structures. They are described below. PMPUP. A standard testbed for P-Median Problem with Users Preferences from the well-known facility location benchmark library Discrete Location Problems (see http://www.math.nsc.ru/AP/ benchmarks/Bilevel/bilevel-eng.html). In total, there are 30 instances (i.e., inst-333, inst-433, inst-533, ..., inst-3233 ) . RND. Medium-scale instances that are randomly generated. Candidate facilities and customers are generated from a two-dimensional uniform distribution [0, 100] 2 . The distance between facility j and customer i (i.e., l ij ) is represented by the Euclidean distance. Each customer has a demand d i , following a uniform distribution [0, 10]. The cost of serving customer i from facility j is computed as c ij = d i · l ij . The disutility g ij is randomly generated from [(1 − δ) · l ij , (1 + δ) · l ij ], where 0 ≤ δ < 1. That is, we allow g ij to deviate 100δ% from l ij , and thus, the nearest facility is not necessarily the most attractive facility to customers. Furthermore, we consider the following combi- ORLIB. Large-scale problems with 1000 customers from ORLib's facility location benchmark set, i.e., capa, capb, capc. They can be found at http://people.brunel.ac.uk/~mastjjb/jeb/orlib/ uncapinfo.html. The service cost c ij and the demand d i are given. We compute l ij = c ij /d i . Similar to RND instances, we generate g ij from [(1 − δ) · l ij , (1 + δ) · l ij ], where 0 ≤ δ < 1. Table 2 reports the computational time in seconds (CPU[s]) of the four solution approaches over the 30 instances in the PMPUP dataset. For each instance, CPU for the approach performing the best (i.e., the smallest value among PDRM, SRM, Benders, and Benders-AS) is highlighted in boldface. The last two rows show the average CPU (AVG) and the average relatively improvement (ARI) of Benders-AS over the others, computed as which indicates how much Benders-AS is faster than the benchmark approach in terms of the AVG. According to the table, PDRM is significantly slower than the others. As shown in (A.4), the PDRM model introduces a large number of additional (dual) variables and constraints; therefore, the size of the formulation increases dramatically, which adds to the difficulty of handling the problem. By contrast, the SRM model in (4) works in the original decision space and thus has substantially fewer variables and constraints. This explains why SRM is on average more than 2 times faster than PDRM in this dataset. Surprisingly, the performance of SRM is comparably good when benchmarking with Benders: The AVG of SRM is only slightly higher than that of Benders (273.9 seconds versus 263.4 seconds). Moreover, for these 30 instances, SRM has lower CPU in 16 instances. Therefore, employing the standard branch-and-cut Benders decomposition cannot effectively speed up the computation. In effect, the efficiency of the Benders-based approaches, to a large extent, depends on the speed of generating Benders cuts (i.e., the separation of Benders cuts) when a master solutionx is found. In Benders, such a process involves solving the dual subproblem [DSP i (x)] using external solvers, which may not be efficient since the solver requires compiling time and [DSP i (x)] itself is a large-scale LP that is not trivial to be solved by the LP algorithms. Considering that the separation of Benders cuts is typically performed a large number of times before the optimality is verified, it is within our expectation that Benders could be significantly slow down owing to the extensive efforts made on the separation. Fortunately, Section 3.2 has provided an analytical and efficient method for the separation. In general, when we apply the analytical separation, we observe substantial performance improvement, i.e., the AVG of Benders-AS is 67.23% and 60.87% shorter than SRM and Benders, respectively (see the last row of Table 2 ). Therefore, the branch-and-cut Benders decomposition should be supported by the well-designed analytical separation to be more practically powerful. We proceed to the result analysis on the RND instances. We drop PDRM since its performance is rather limited. The full computational results are provided in Appendix B. We summarize the main results in Figure 2 . Figure 2(a) shows the percentage of instances that are solved optimally (# instances solved) within any given CPU [s] . A point in the figure with coordinates (m, n) indicates that for n % of the instances, the required CPU is less than m seconds to solve them to optimality. According to Figure 2(a) , Benders-AS outperforms the others since the green solid-line is consistently above the others, meaning that given the same CPU on the x-axis, Benders-AS can successfully solves more instances. Moreover, using Benders-AS, all instances can be solved within 1000 seconds. By contrast, Benders requires roughly 2500 seconds to clear all instances, and SRM even fails in several instances within the 7200-second time limit. Therefore, in terms of the number of instances solved optimally, both Benders-based approaches outperform SRM in this dataset. Figure 2 (b) reports the boxplot of the computational time (in log-scale) for the above three approaches under two values of δ. The line in the box indicates the median of CPU under a specific value of δ. Note that a larger δ value generally indicates that customer preferences differ more from the operator's unit service costs. The difficulty of the instances certainly depends on δ, and intuitively, when δ = 0, the model reduces to the standard facility location problem without user preferences, which can be much more efficiently solved. Thus, we expect that under a larger δ value, the instances will be more challenging. Apparently, the results in Figure 2 agrees with our intuition as we observe that the instances requires longer CPU when δ = 0.5. Moreover, in terms of CPU, Benders-AS outperforms the others by roughly one order of magnitude, whereas the advantage of Benders over the SRM is still not significant. This observation once again highlights the effectiveness of the proposed analytical separation in expediting the branch-and-cut Benders decomposition. Our next experiment is to investigate the performance of the approaches on the large-scale ORLIB structured data (capa, capb, capc) . Table 3 reports the results. Here, the column rgap[%] is the relative exit gap in percentages when the solution process terminates. It is computed as |zbb − zopt|/|zbb| × 100, where zopt is the value of the current optimal feasible solution and zbb is the best bound at termination. If an instance is solved to optimality, then zbb = zopt or rgap < 0.01%. The column Ratio is computed by "the CPU of an approach divided by the CPU of Benders-AS" when the corresponding instance is solved optimally. For example, for capa under δ = 0.1 and P = 5, the two Ratio values under SRM and Benders are 89.8 and 62.7, stating that the CPUs by these two approaches are 89.8 and 62.7 times of the CPU by Benders-AS. Moreover, for each problem instance, the CPU for the approach performing the best is highlighted in boldface. For those instances where all approaches cannot terminate optimally within the time limit, the smallest rgap is highlighted in boldface. According to Table 3 , we have the following observations: (i) Consistent with the previous finding, the CPU increases with δ, i.e., when the user preferences differ more from the costs of the operator, the instances becomes harder. In particular, when δ = 0.5 and P = 10, the capa and capc instances are so changeling that none of the approaches can successfully solve them to optimality within the time limit; (ii) For all instances, Benders-AS is the best approach in this dataset. When the instance can be solved optimally, the CPU by Benders-AS is generally more than one order of magnitude shorter. This can be directly observed from results in the column Ratio. Meanwhile, when the instance cannot be solved optimally within the time limit, the rgap of Benders-AS is the smallest, and the value is less than 1%, which has been small enough to guarantee the solution quality. Based on the results of the above three experiments, we conclude that Benders-AS significantly improves the performance of the standard branch-and-cut Benders algorithm and outperforms the benchmark approaches by a large margin. In our final experiment, we investigate the importance of integrating user preferences into the facility location problem when the preference does exist. Specifically, we compare the total service costs of the operator (i.e., the objective function φ) under two scenarios: (i) the operator considers user preferences and anticipates user's choices when locating facilities. This scenario is exactly the [SRM] model presented in Section 2; (ii) the operator ignores user preferences and locates facilities based on the cost matrix c. In this scenario, the operator's decision is made based on the classical P-Medium problem (see Chapter 6 of Daskin (2011)). To facilitate our discussion, we introduce the following notations. Let x wt be the location decision when the operator locates facilities without considering user preferences. x wt can be obtained by setting g = c or simply deleting Constraint (4d) in the model. The operator's actual total service cost φ wt is evaluated by holding the location decision at x wt . Therefore, φ wt stands for the total service cost when the decision is made based on the assumption that user preferences do not exist (or are mistakenly assumed to coincide with the operator's service cost). Moreover, the optimal total service cost φ is the optimal objective function value when the operator indeed accounts for user preferences (i.e., the best objective of [SRM]). Then, we define the relative cost increase ∆[%] as which specifies the relative additional cost incurred when user preferences exist but is ignored by the operator when locating facilities. Our experiment is conducted on the PMPUP dataset since the original problem presented in the benchmark library has a similar structure to ours. Table 4 reports the operator's total service costs under three values of P . The first P value is set to 14 since the original dataset imposes P = 14 for all instances (except for Inst-533, marked with * in the table, where P = 13). We observe that for the original instances (i.e., the column "P = 14"), the average ∆ is 2.44%; therefore, ignoring user preferences will lead to an additional 2.44% of the average service cost. Furthermore, we observe that ∆ increases dramatically when P increases. In particular, the average ∆ under P = 30 blows up to 84.24%. This result is astonishing as it indicates that the operator will almost double its service cost; therefore, the preferences must be considered if they exist. Another interesting observation is that when user preferences are ignored, opening more facilities may lead to a higher service cost. For example, in the Inst-333, when the operator opens 20 facilities, φ wt is 139; however, when 10 more facilities are open, φ wt increases to 142. To give a direct view, we plot how φ and φ wt evolve with P in Figure 3 . The figure shows φ wt could increase when P increases. On the contrary, we can clearly observe the decreasing trend of φ, stating that when the preferences are considered in the operator's decision stage, opening more facilities will indeed reduce the service cost. To summarize, the above experiment demonstrates that the operator must take user preferences into consideration and correctly anticipate user's choices; otherwise, the operator could suffer from a substantially higher cost, and opening more facilities could unexpectedly result in additional service costs as well. These results further justify the usefulness of the P-median problem with user preferences. There are limitations and potential future research directions. Firstly, as shown in Fischetti et al. (2017) , one may further enhance the Benders decomposition by designing a proper cut loop strategy and stabilization at the root node. However, fractional location solutions arise when stabilization is implemented. Noting that the analytical separation is only applicable to integer nodes, we are thus unable to generate Benders cuts efficiently when stabilizing the solution, and an alternative method must be used instead. This means that an efficient cut loop is not ready to been obtained, and its effectiveness on the algorithm also needs more careful investigations. Secondly, the [SRM] formulation leveraged a closest assignment constraint (CAC) to transform the bilevel model into a MILP. The whole Benders approach was built upon this formulation. It is then interesting to explore other CACs since their strengths and sizes are different (Espejo et al., 2012) . It is possible to have a powerful Benders algorithm that is developed based on other CACs. Finally, since the operator needs to forecast customer preferences to anticipate the choices, there exists the possibility that the forecasted preference is subject to estimation errors. In this case, we should build a stochastic model or a robust optimization model, and we would like to leave the model formulation and the algorithm development for future research. This appendix presents a single-level reformulation model for [PUP] , which relies on the primaldual optimality conditions of the lower-level problem. Similar approaches can be found in Casas-Ramírez and Camacho-Vallejo (2017), Casas-Ramírez et al. (2018) . Here, we briefly describe the reformulation procedure. Given the operator's location decision x, the lower-level problem is min y∈R + i∈I j∈J π ij y ij (A.1a) st. j∈J y ij = 1 ∀i ∈ I (α i ) (A.1b) y ij ≤ x j ∀i ∈ I, j ∈ J (β ij ) (A.1c) The α i and β ij in the parentheses are the dual variables associated with the constraints. We can then rewrite the lower-level problem with its KKT conditions, i.e., j∈J y ij = 1 ∀i ∈ I (A.2a) y ij ≤ x j ∀i ∈ I, j ∈ J (A.2b) y ij ≥ 0 ∀i ∈ I, j ∈ J (A.2c) α i + β ij ≤ g ij ∀i ∈ I, j ∈ J (A.2d) where (A.2e) and (A.2f) are bilinear functions. Noting both y ij and x j will take 0/1 in the optimal solution and the maximum value of the matrix π is 1 by definition, these bilinear functions can exactly linearized by which is referred to as the primal-dual-reformulation model. In our preliminary computational experiment, we observe that setting y ij ∈ {0, 1} (which is feasible and equivalent) makes Gurobi run faster and enhances the numerical stability as well. This appendix presents the computational results of the RND instances under the 2-hour time limit in Table B .5. The RND dataset is generated with fixed random seed in Python files and can be provided upon request. Optimal location with equitable loads A lagrangean relaxation algorithm for the simple plant location problem with preferences A two-stage location problem with order solved using a lagrangian algorithm and stochastic programming for a potential use in covid-19 vaccination based on sensor-related data A matheuristic for solving the bilevel approach of the facility location problem with cardinality constraints and preferences The p-median bilevel problem under preferences of the customers. Recent Advances in Theory Solving the bilevel facility location problem under preferences by a stackelberg-evolutionary algorithm A strengthened formulation for the simple plant location problem with order Capacitated plant selection in a decentralized manufacturing environment: a bilevel optimization approach Solving the p-median bilevel problem with order through a hybrid heuristic Approximating solutions to a bilevel capacitated facility location problem with customer's patronization toward a list of preferences Benders decomposition for very large scale partial set covering and maximal covering location problems Network and discrete location: models, algorithms, and applications Grasp and hybrid grasp-tabu heuristics to solve a maximal covering location problem with customer preference ordering Closest assignment constraints in discrete location problems Benders decomposition without separability: A computational study for capacitated facility location problems Redesigning benders decomposition for large-scale facility location A facility location problem with clients' preference orderings Lower bounds for the uncapacitated facility location problem with user preferences. Groupe d'études et de recherche en analyse des décisions Facility location and scale decision problem with customer preference Branch-and-cut approach based on generalized benders decomposition for facility location with limited choice rule Exact approaches to the single-source network loading problem Robust bi-level programming for renewable energy location Metaheuristic methods for solving the bilevel uncapacitated facility location problem with clients' preferences Facility location and supply chain management-a review A variable neighborhood search for the budget-constrained maximal covering location problem with customer preference ordering Benders decomposition for the profit maximizing capacitated hub location problem with multiple demand classes New lower bounds for the facility location problem with clients' preferences The branch and cut method for the facility location problem with client's preferences This paper studied the exact solution approach for the P-median facility location problem with user preferences (PUP). By exploring the problem structure, we proved that in a CAC-based MILP model (Casas-Ramírez and Camacho-Vallejo, 2017), the high-dimensional binary variables used to model customer preferences and facility choices can be relaxed to continuous variables. Based on this, we proposed a branch-and-cut algorithm where Benders separation (i.e., the procedure of generating Benders cuts) at integer nodes of the searching tree was performed leveraging external solvers. However, such a standard Benders approach was not efficient enough. Therefore, we further proposed an acceleration technique to enhance the algorithm performance. Using a broad testbed, our computational experiments demonstrated that the proposed algorithm outperformed several benchmark approaches by a large margin and was able to handle large-scale instances satisfactorily. We also conducted sensitivity analysis and observed that when user preferences indeed exist, the operator must consider customer preferences and correctly anticipate the choices to avoid an unnecessarily high service provision cost. Furthermore, ignoring the preferences may lead to an ironic situation where opening more facilities could result in additional service costs.