The facility location problem is one of the most fundamental problems investigated in operations research. In this problem, a decision maker needs to decide locations of a limited number of facilities (factories, retail centers, power plants, service centers, etc.), and determine coverage of demand from different sites by the located facilities. The objective is to minimize the facility setup cost and the cost of production and delivery. This problem provides a basic framework to formulate related problems in resource allocation, supply chain management and logistics, etc. The facility location models in the stochastic and robust (distributionally-robust) optimization framework have also received sufficient investigation. In the stochastic optimization framework, the customer demand can be random parameters following a known or partially known probability distributions, while in the robust (distributionally-robust) optimization framework customer demand (the probability distribution of customer demand) can have a certain level of uncertainty (ambiguity). Incorporating randomness or uncertainty in customer demand estimation is a sensible and realistic model improvement, as the decision maker has no perfect prediction on the demand in practice. Recently, there is a trend of research on distributionally-robust optimization problems with decision-dependent ambiguity, in which the ambiguity set is specified by decision dependent parameters. Introducing decision dependency has certain merit in situations when model parameters are naturally functions of decision variables. For example, in a basic pricing and revenue management problem the demand of a product can be a decreasing function of the price to be determined. The notion of modeling decision-dependency in the robust (distributionally-robust) optimization framework can be applied to investigate a variety of specific problems.
In this paper, we investigate a distributionally-robust service center location problem with decision-dependent customer demand. In the model, customers need to physically access a service center in order to receive service. This problem is motivated from clinic or medical test center (such as for COVID-19 screen test) allocation. The objective is to maximize the total utility gain of all customers who have decided to receive service from the opened service centers. The decision dependence is decoupled from the ambiguity set in our model, which is different from [1] . Specifically, we impose the decision dependency on the demand itself but not on the uncertainty of demand, and the decision-dependent demand is induced from a maximum attraction principle based on a ranking of opened service centers in the neighborhood. Our investigation shows that the decoupling approach leads to high computational efficiency of solving the problem and it can be more data-driven in practice. Furthermore, the novel approach of modeling decision dependency introduce additional combinatorial properties to the model which is of independent interest.
The contribution of this paper is summarized as follows:
1. A novel approach of modeling decision dependency of the customer demand has been established for the service center location problem, with a possible extension on the modeling of decision dependency. Tractable reformulations and solution strategies of this problem in the capacitated or incapacitated cases are investigated. 2. A learning model has been established to estimate utilities and demand based on data from survey. 3. Submodularity of the objective function and polyhedral properties of the polytope induced by the demand and capacity constraints have been well studied for independent interest. 4. The numerical experience with this model shows that mid-size instances can be solved very efficiently due to decoupling of decision dependency from the ambiguity set.
Facility location problems with non-deterministic demand have received plentiful investigation in the framework of stochastic optimization and robust (distributionally-robust) optimization. In the stochastic optimization framework, customer demand are independent random parameters with known probability distributions and the problem can be formulated as a two-stage stochastic program in which the location vector is in the first-stage decision that needs to be made before realization of customer demand [2, 3] . In a special case when the customer demand rates and service rates of each facility are assumed to follow exponential probability distributions, the problem of minimizing long term average cost can be reformulated as a deterministic optimization problem with corresponding Poisson rates to characterize the demand and service. This setting has been applied to an allocation problem of ATMs [4] . The probability distribution of demand can also be used to define chance constraints to ensure a required service level under possible stockout and supply disruption [5] [6] [7] [8] .
In the robust (distributionally-robust) optimization framework for facility location problems, the information of demand is partially known to the decision maker, and the goal is to find a robust optimal location vector that optimizes the objective after a worst-case realization of demand information [9] . Following this direction of research, Baron and Milner [10] investigated a robust multi-period facility location problem with box and ellipsoid sets of uncertainty. Gourtani et al. [11] investigated a distributionally-robust two-stage facility location problem with an ambiguity set defined corresponding to the mean and covariance matrix of a random parameter for expressing the demand. Facility location problems with demand uncertainty have been studied with a variety of novel application background, which includes but not limited to medication coverage and delivery under a large-scale bio-terror attack [5] , medical equipment (defibrillators) location problem to reduce cardiopulmonary resuscitation (CPR) risk [12] , humanitarian relief logistics [13] , and hazardous waste transportation [14] , etc.
In a lot of real world problems, the uncertainty of parameters can interplay with the decision to be made. This behavior is well observed especially in a sequential (multi-stage) decision-making process, in which information about system parameters are gradually revealed and the decisions made up until the current stage can reshape the uncertainty in future [15] . It motivates the research on multi-stage stochastic optimization with decision-dependent uncertainty. Solution strategies based on Lagrangian duality and novel branchand-bound methods are developed to solve this family of problems [16] [17] [18] , and it has a rich application in oil-chemical industrial [19] [20] [21] [22] . An approximation scheme [23] is proposed to tackle the high complexity in solving the multi-stage problems with decision-dependent information discovery. Decision-dependent uncertainty has also been considered in stochastic optimization problems such as resource management [24] , stochastic traffic assignment modeling [25] , and robust network design [26, 27] .
Imposing decision-dependent uncertainty for robust (distributionally-robust) optimization has received great attention in recently years [1, [28] [29] [30] [31] . Nohadani and Sharma [28] studied robust linear programs with decision dependent budget-type uncertainty and its generalization with a polyhedral uncertainty set. This concept is demonstrated in a robust shortest-path problem, where the uncertainty is resolved progressively when approaching the destination. Noyan et al. [30] investigated a family of distributionally-robust optimization problems with an ambiguity set defined using earth mover's distances (including total variation distance and the Wasserstein metric) with decisiondependent parameters such as the nominal probability distribution and the radius, and focused on understanding which settings can lead to tractable formulation. Royset and Wets [32] provided a variational principle analysis for optimization under stochastic ambiguity, which gives fruitful tools for analyzing solution quality and price of robustness. Some novel applications are studied in [33] for robust software partition, and in [34] for radiation therapy design.
The work of Basciftci et al. [1] and our previous work [35] are mostly related to the work of this paper. Basciftci et al. [1] investigated a distributionallyrobust facility location problem with decision-dependent ambiguity, where the ambiguity set is defined using disjointed lower and upper bounds on the mean and variance of the candidate probability distributions for customer demand, and these bounds are defined as linear functions of the location vector to admit a tractable MILP reformulation. The model considered in [35] assumed deterministic customer demand and imposed a linear decision dependency on the moment bounds for candidate probability distributions of utility, which leads to a MISOCP reformulation.
We consider a problem of allocating a set of service centers (facilities) in a region to meet the customer demand from the region. Similarly to the traditional facility location problem, the region consists of customer sites and candidate locations for the service centers. Each customer site has certain demand that needs to be fulfilled by a service center from its neighborhood. Especially, we consider the case that the demand of each customer site is a random parameter which depends on the locations of service centers. In most literature on decision-dependent robust (distributionally-robust) optimization, the decision-dependency is formulated as a linear function. In this problem, we consider a special form of decision dependency that is based on an assumption of maximum attraction principle, which will be clear in Section 2.1.
We first formulate a deterministic model of this problem, in which we treat the demands as deterministic parameters that depend on the locations of ser-vice centers and introduce the notion of maximum attraction principle. The notations used in this model is given by the following list:
S the set of customer sites; F the set of candidate locations for service centers; b j the cost of opening a service center at location j ∈ F ; B the budget of allocating service centers; C j the service capacity of a service center located at j ∈ F ; u ij the utility gain obtained by a customer from site i ∈ S who gets service from the service center at j ∈ F ; y j the binary decision variable of opening a service center at location j ∈ F ; D i (y) the demand from customer site i that depends on the decision vector y; x ij the demand flow from customer site i ∈ S to the service center at j ∈ F . The deterministic model is formulated as follows:
The first three constraints represent the limit of budget, capacity and the amount of demand, respectively. Note that in general the demand from site i ∈ S can depend on the pattern of service center locations which is represented by the decision vector y. In this paper, we investigate a specific type of location dependent demand which is simple but sensible in depicting customers' behavior in practice. This type of location dependency of demand is referred as the maximum attraction principle in this paper. We now give a detailed description of this principle. First, each customer site i is associated with a preferable subset F i of candidate facility locations. When there is only one service center opened at j ∈ F i , it will attract D ij amount of demand from i. If there are multiple service centers opened at locations in F i , the maximum demand that can be attracted from i is equal to the maximum D ij for j ∈ F i that has a facility. This principle is formally described by the following equation:
where e j is the |F i |-dimensional vector with the j-th entry being 1 and other entries being 0, and the set F i (y) is defined as F i (y) := {j ∈ F i | y j = 1}. The maximum attraction principle matches with our intuition from practice. This principle first assumes that locations of service centers that are not within the preferable location set F i are not attractive to customers from site i at all. The most natural way of establishing F i is based on the distance from site i to the candidate locations. This is a reasonable assumption because customers usually will not consider visiting a service center that is beyond a certain distance from their living place. The maximum attraction principle further assumes that within the preferable opened service centers F i (y), there exists one service center j * (or multiple centers) that is (are) most attractive to the customers from site i in the sense that j * ∈ argmax j∈Fi(y) D ij , and the presence of multiple service centers in F i (y) including j * attracts the same amount of demand from site i as the presence of just a single service center at j * . This assumption is also consistent with our intuition from practice. If every service center is identical in the sense of scale and service quality, the most attractive one to site i is likely the one that is most close to i. Furthermore, customers who are willing to visit further service centers are also willing to visit the one that is most close to their living place. Note that the maximum attraction principle can be extended to incorporate the impact of the number of opened service centers in F i on the demand. This extension is discussed in Section 2.2. We also define the subset S j as S j := {i ∈ S | j ∈ F i } for all j ∈ F. Based on the maximum attraction principle, the service center location problem (1) is written as the follows:
To go one step further, we can linearize the term max j∈Fi(y) D ij by introducing some continuous auxiliary variables q ij to form it as a convex combination of D ij for j ∈ F i . After this transformation, we obtain the following equivalent formulation of (3):
Notice that in the above reformulation, the auxiliary variables {q ij | j ∈ F i } are used select the most attractive candidate location driven by the sense of maximizing the objective. The constraint q ij ≤ y j ensures that only opened service centers in F i are involved in the maximum attraction principle. Since j∈Fi q ij ≤ 1, the model will set q ij * = 1 and q ij = 0 for all j ∈ F i \ {j * } to relax the constraint j∈Fi x ij ≤ j∈Fi D ij q ij as possible.
In the model setting, we implicitly assume that the customers are willing to corporate with the decision maker to maximize the total utility gain. Although this assumption is highly impractical, the rationality of this model depends on what metric we use to measure the system performance. If the goal is to estimate what is expectation of total utility in practice, then the following questions should be addressed and the corresponding aspects should be properly modeled: What is the service policy used by each service center (FIFO or some other policies)? How to characterize customers' behavior and the mechanism of competing for the limited service capacity? If their most favored service center does not have any capacity, are they willing to accept the service from the less preferred locations with less utility gain? Incorporating all these factors into the model can easily make it very complicated, and meanwhile a large amount of customers' information are needed to drive this approach of modeling. On the other side, if the goal is to access what is the maximum potential utility that can be achieved by the system in the most ideal situation, then the (DDSL) model can be used to give an estimation. A possible modification of modeling the decision-dependent demand is to add a perturbation term to the demand based on the number of opened service centers. Specifically, we can modify the constraint i∈Fi x ij ≤ max j∈Fi D ij to be i∈Fi x ij ≤ max j∈Fi D ij + a i ( j∈Fi y j ), where a i is a parameter that measures the influence of opening one more service center on the demand. In this way of modeling, it is assumed that the marginal increase of demand may depend on the number of opened service centers. In reality, when potential customers see more chain stores are opened in the neighborhood, they may have higher intention to try one of them. In this case, the parameter a i can be a random parameter with certain level of ambiguity in the distributionallyrobust extension of the model.
The deterministic service center location model with decision dependent demand can be further extended to a distribtuionally-robust two-stage stochastic program after imposing an ambiguity set on the pairwise demand D ij . This extension is motivated by the fact that estimation of the demand parameters D ij could be inaccurate. Therefore, we can treat D ij as random parameters with an unknown joint probability distribution, and apply the distributionallyrobust framework on this service center location problem with uncertainty. We assume a finite support of the joint probability distribution of the pairwise random demand vector D := {D ij | i ∈ S, j ∈ F i }. Let the finite support be based on |Ω| samples written as D ω = D ω ij i ∈ S, j ∈ F i for all ω ∈ Ω. In this case, any probability distribution of D can be represented as a |Ω|-dimensional vector. We define a nominal probability distribution µ 0 of D as
In the vector representation, we write it as µ 0 = [µ ω 0 : ω ∈ Ω]. The ambiguity set of candidate joint probability distribution of D is defined based on the total variation distance between two probability distributions. Specifically we consider an ambiguity set of the following form:
Based on the above definition of ambiguity set, the distributionally-robust two-stage stochastic extension of the service center location model can be formulated as follows: max
where the recourse function Q(y, D ω ) for scenario ω is given by
Using a standard technique for two-stage stochastic programming, we can decompose (DRO-FL) into a master problem and scenario sub-problems. The original problem (DRO-FL) can be solved iteratively. In each iteration, we solve the current master problem in the space of y and pass the current master solution to each scenario sub-problem. After solving each scenario sub-problem for the fixed first-stage solution, we can generate a valid inequality for each scenario using the optimal dual values associated with constraints of the scenario sub-problem. Then we aggregate the valid inequalities from all scenario sub-problems using the worst-case measure of scenarios to get a single cut for the master problem which is an optimality cut. The optimality cut is added to the master problem in the next iteration. Specifically, the master problem at iteration n can be represented as max η
where µ (k) is the iteration based worst-case probability measure on scenarios at iteration k ∈ [n − 1]. The way of determining this worst-case probability measure is given in (10) . At iteration n we solve the master problem (Master) and obtain the current optimal first-stage solution y (n) . This y (n) is input into every second-stage scenario sub-problem. The second-stage linear program is solved and let α
ω,ij ≥ 0 for all i ∈ S, j ∈ F i are optimal dual values corresponding to the constraints (7b), (7c), (7d) and (7e), respectively. Using the standard technique of Bender's decomposition [cite ref], we can obtain the following valid inequality on the value function Q(y, D ω ):
The strong duality implies that when evaluating Q(y, D ω ) at y (n) , we get
The worst-case probability measure µ (n) is an optimal solution of the following linear program: min
where Q(y (n) , D ω ) is the value function evaluated at the first stage solution y (n) and scenario ω. Once the current worst-case probability measure is obtained, the following inequality will be added to the first-stage master problem:
The algorithm and convergence property for solving (DRO-FL) are given in Appendix A. It can be shown that for a given first-stage solution y, the scenario sub-problem can be solved using a greedy algorithm.
We consider a special case of (DRO-FL) in which each candidate service center has sufficient capacity for service. In this case, the capacity constraints (7b) can be removed from the second-stage scenario problems. We will show that this leads to a simplified scenario problem that admits a closed form optimal solution and optimal objective for the second-stage problem with a mild regularity condition on the parameters. Then after dualizing the inner minimization problem over the probability measure on scenarios, (DRO-FL) can be reformulated as a mixed 0-1 linear program. The single-stage reformulation result is given by the following theorem.
The utility gain and demand are consistent if for every ω ∈ Ω, i ∈ S and F ⊆ F i there exists a j ∈ F such that u ij = max j ∈F u ij and
The consistency condition for utility gain and demand says that the location that attracts the most demand over other locations should also correspond to the highest utility.
Theorem 2.1 Suppose the utility gain and demand are consistent. In the uncapacitated case, the distributionally-robust service center location problem (DRO-FL) with decision dependent demand based on the maximum attraction principle and the ambiguity set (6) can be reformulated as the following mixed 0-1 linear program:
Proof Without capacity constraints, it is easy to see that the optimal objective value of the scenario problem is given by
where we use the assumption that the utility gain and demand are consistent. The problem (DRO-FL) becomes the following:
The terms involved in the inner minimization problem of (14) can be written as min
Taking the dual of the above linear program with respect to the probability measure µ, we obtain the following inner problem:
To linearize the term max j∈Fi u ij D ω ij y j , we can introduce binary indicator variables s ij , and reformulate the first constraint as follows i∈S j∈Fi
where we implicitly use the consistency condition of utility and demand which implies argmax j∈Fi u ij D ω ij y j is scenario independent and hence the variable s ij . Incorporate with the outer maximization yields the reformulation (12).
We establish a regression model based on the maximum attraction principle for estimating utility value u ij and demand D ij that are input parameters to the service center location model (DRO-FL). As a by-product, this regression model can also be used to estimate the size parameter d of the ambiguity set (6) . The regression model requires samples of survey among potential customers from a site i ∈ S on their ratings of multiple candidate service center locations. Before formulating the regression model, we first introduce how the samples of survey are collected for fitting the regression model. Suppose N residents have been randomly selected from the site i, and they are viewed as potential customers of service centers under planning. Each of them is asked to give a score in the range {0, 1, . . . , q} to each candidate service center location in F i . The score that a customer is assigned to a location j ∈ F i is taken as the potential utility gained by the customer if going to the service center at j. Score value 0 means the customer is unwilling to get service from the corresponding location. Since the regression model structure is identical for each i ∈ S, we omit the customer site index i and re-write F i as G in the following modeling and analysis. Suppose the indices of location in G are labeled as G := {1, 2, . . . , g}, the N residents are labeled as {1, 2, . . . , N }, and the score assigned to a candidate location j ∈ G by the resident k is denoted as a kj .
The maximum attraction principle described in (2) may not be satisfied in reality. But it is possible to establish a utility-demand estimation model that approximately meets the maximum attraction principle. For example, we can verify whether the samples from survey satisfy the maximum attraction principle by grouping the customers who have taken the survey as follows:
where [N ] := {1, . . . , N }. In words, V j is the set of residents who are willing to go to a service center located at j. According to (2) , the samples satisfy the maximum attraction principle exactly if there exists a permutation σ on the indices in G such that the following inclusive condition holds:
A logic behind the maximum attraction principle is that if two candidate locations have similar features, then if a customer is attracted by one location, the customer should also be attracted by the other location with high chance. In this case, the difference in the attractability of the two candidate locations viewed by the customer is reflected in the score assigned to the two locations by the customer. The scores in this case are both non-zero indicating that the customer is willing to visit any of them. Only in the case that two candidate locations have some substantial differences (i.e., one is too far away from the customer site or one is located at a bad community), a customer will be willing to visit one location (assigning a non-zero score to it) while unwilling to visit the other one (assigning a zero score to it).
We establish an inclusive chain model to represent the set-level realization of the maximum attraction principle. First we build the subsets (16) using collected samples, and the we sort the indices in G such that
where σ is a permutation on G that makes this condition hold. For clarity we assume that |V 1 | ≤ |V 2 | ≤ . . . ≤ |V g | without loss of generality. The inclusive chain model assumes that the score assigned to the candidate locations in G by a customer from a fixed location should match with one of the following patterns:
where in the i-th pattern, first i − 1 scores are all zero and the remaining g − i + 1 scores are all non-zero. By convention, the (g + 1)-th pattern is just [0, . . . , 0]. In a probabilistic flavor, the inclusive chain model (ICM) for the score vector can be formally established as a statistical model presented as
(ICM) where ξ is the random score vector, B is an indicator random variable that selects the pattern matching with the score vector, p j is the probability of score vector matching with the j-th pattern, and π (i) jr r = 1, . . . , q gives the probability distribution of the random score Q Proof We prove it contradiction. Suppose there exist subsets V i and V j (with i < j) such that V i \ V j = ∅ with some positive probability. Suppose k is the customer who is willing to visit V i but not V j , and let ξ k be the score vector of this customer. The definition of V i implies that ξ i ≥ 1 but ξ k / ∈ V j implies ξ k j = 0. On the other side, for every pattern vector ζ we should have 1{ζ i > 0} ≤ 1{ζ j > 0} almost surely, which implies that 1{ξ k i > 0} ≤ 1{ξ k j > 0} almost surely. But this contradicts to ξ k i ≥ 1 and ξ k j = 0.
We consider the problem of fitting the inclusive chain model with the collected samples (score vectors). The first step is to determine the rank (level of attracability) of candidate locations G involved in the inclusive chain. An empirical method by simply sorting the cardinality of V i can be used to achieve this. The probability guarantee of this method which will be discussed in a moment. For now, assume that the order has been identified, and suppose |V 1 | ≤ |V 2 | ≤ . . . |V g | without loss of generality. So empirically, this implies that the rank is 1, 2, . . . , g sorted by level of attractability from low to high. Based on this information, we can establish an inclusive chain model as (ICM), but this model is not capable to handle score vectors in which there is at least one zero-value entry between two non-zero value entries, i.e., [0, . . . , 0, 5, 7, 0, 0, 7, 6, 9] . We call a score vector defective if it can not match with any pattern vector in (18) valid for the inclusive chain model. Note that to convert the defective score vector [0, . . . , 0, 5, 7, 0, 0, 7, 6, 9] into a qualified pattern one can either add two 1 score entries to the middle ([0, . . . , 0, 5, 7, 0, 0, 7, 6, 9] → [0, . . . , 0, 5, 7, 1, 1, 7, 6, 9]) or remove the first two nonzero entries ([0, . . . , 0, 5, 7, 0, 0, 7, 6, 9]→[0, . . . , 0, 7, 6, 9] ). The first approach has added a total 2 units of score while the second one has removed a total 12 units of score, and it can be verified that the first approach is the one that makes the minimum total score change among all possible changes that can convert this score vector into a qualified pattern. Based on this rule, we can define the defective score units for an arbitrary score vector as follows.
Definition 3.1 For a given arbitrary score vector v, the defective score quantity is given as
w is a qualified pattern of the inclusive chain model.
Observation 1 The defective score quantity has an upper bound g − 1.
Let H s be the event that the defective score quantity of score vector is s, and m s is the probability of this event. After introducing these additional parameters, the inclusive chain model can then be adjusted to incorporate all defective score vectors as follows (AICM) The above model is referred as the adjusted inclusive chain model (AICM).
As mentioned at the beginning of Section 3.2, a simple sorting method can be used to identify the rank of the candidate locations in G for the AICM. We sort the cardinality of subsets V j (j ∈ G) defined as
where a k is the score vector of the k-th customer participated in the survey, and a kj is its j-th element. The following proposition provides a probability guarantee on identifying the order of inclusive chain model given that samples obey the AICM. 
Proof We first show that with the given assumptions, the rank of the inclusive chain model can be identified almost surely as the total number of samples N goes to infinity. In the analysis, we let σ be the permutation that gives the rank of locations, i.e., σ(1) < . . . < σ(g) is the rank of locations from low to high. Let N j be the number of samples that match with the j-th pattern
Based on the probability setting of AICM, and the law of large numbers we should have the following asymptotic relation hold:
This implies that N j > 0 almost surely as N → ∞. By the definition of V j , we should have the following bounds for every
For
> 0 by assumption. Therefore the rank of the inclusive chain can be identified almost surely as N goes to infinity. Now we prove the non-asymptotic result. Note the probability P (E) can be lower bounded as follows
Hoeffding's inequality gives the following bound
Then it follows that
To avoid redundant notations in the following analysis, we assume that candidate locations are ranked as 1 < . . . < g based on attractability. To fit the AICM with samples, a likelihood function and some regularity terms need to be established. Note that without consideration of any regularity, all model parameters can be estimated naturally using a frequency count approach based on the given samples. For example, an empirical estimation of p j can bê
jr is the number of samples that have Q (i) j = r. However, this estimation can lead to overfitting of new samples. We now focus on establishing a likelihood function and introduce some regularity terms to prevent overfitting. First, the parameter m s can be estimated as m s = #{k∈[N ]: ξ k ∈Hs} N for all s ∈ {1, . . . , g − 1}. All defective samples can be further used to estimate other parameters with more information. If a sample ξ k is defective, it can be converted into a qualified score vector ξ k by imposing a defective score quantity (Definition 3.1) to ξ k , and hence the original samples {ξ k } N k=1 can then be converted into a set of qualified score vectors {ξ k } N k=1 . We can evaluate the probability of generating ξ k by the model using the model parameters, and denote this probability as P (ξ k |p, π). Then the (logarithmic) likelihood function of generating {ξ k } N k=1 can then be formulated as
Note that maximizing the above likelihood function with only probabilitynormalization constraints leads to a frequency estimation of all parameters. Some regularity terms can be added to balance this trend. First, for score vectors in the same pattern, parameters that lead to smaller variance of score on a given candidate location are preferred, and hence the following regularity terms can be added (with some regularity coefficients):
∀i ∈ {1, . . . , g}, ∀j ∈ {i, . . . , g}.
Second, consider any two consecutive locations j and j + 1 in the inclusive chain. The two locations are both included in patterns H 1 , . . . , H j . In the pattern H i (1 ≤ i ≤ j), the mean score difference between the two locations is given by
For every pattern that contains locations j and j + 1, we expect the above mean score difference to be similar. To enforce this similarity, we can add the following regularity terms (with some regularity coefficients):
After incorporating the above regularity terms, the parameters can be estimated by solving the following constrained minimization problem of the loss function.
min L(p, π) := 2λ 1 g(g + 1)
s.t. 
where λ 1 and λ 2 are regularity coefficients that can be tuned using a crossvalidation approach. Note that we add a negative sign to the logarithmic likelihood function when it is incorporated into the loss function L(p, π). Notice that (24) is a nonconvex optimization problem with continuous parameters and simple constraints, and it can be solved to local optimality very efficiently using numerical optimization solvers.
Once the regression problem (24) is solved to give an estimation of the AICM parameters, the AICM can be used to estimate the (mean) utility value, guide the generation of demand samples. First, the utility value u ·,j (the customer site index is ignored since the AICM is established for any given customer site) can be estimated as follows:
wherep andπ are the estimation of AICM parameters. Generating a number of joint demand samples and use them to construct the ambiguity set can be a challenging problem in practice. There are two major obstacles in achieving this. (a) Scaling factor: since only a limited amount of survey samples are collected, how to scale the demand estimation based on the samples to roughly estimate the total demand from a customer site? (b) Demand catenation: how to catenate the estimated demand samples for each customer site to get joint demand vectors (for all customer sites) as samples that can be directly input to the ambiguity set? We try to give a few guidelines for dealing with these problems in this paper. For (a), the total potential demand from a customer site could be estimated via a different channel which is not the focus of this paper. Some investigation on the demand estimation and especially the healthcare service demand estimation is conducted in [36] [37] [38] [39] . The survey of scoring on candidate locations are only used to quantify the attractability of each candidate location that can lead to different estimated demand for different candidate locations when the total potential demand is given. In the next paragraph for addressing (b), we provide a estimation of demand associated with each candidate location based on assuming the total potential demand is N (the number of collected samples). In the case that we have a separate estimation A (can be a random variable) of total potential demand the scaling factor can be defined as A N , and it can be used to scale up estimated demand associated with different candidate locations. In the following discussion for (b), we assume that the demand samples generated (using Algorithm 1-2) will be multiplied by the scaling factor (could be customer site dependent), and omit this factor in these algorithms.
For (b), one approach is to assume the distributions of demand D i := {D ij | j ∈ F i } for different customer sites i are independent. In this case we can generate demand samples for each i and combine them (in an outer-product manner) to construct samples for the random vector D := {D ij | i ∈ S, j ∈ F i }. For example, if we have constructed |Ω i | samples for D i which are denoted as {D ω i | ω ∈ Ω i }, then we can get samples for D as {[D ωi i : i ∈ S] | ω i ∈ Ω i ∀i ∈ S}. In this manner, the number of samples for D can be as large as i∈S |Ω i |. To reduce the sample size for D, one can select a desired number N of combinations of samples for each D i at random. Another way is to assume the demand from different costumer sites are strongly correlated. As a consequence, we can divide the samples for each D i into several bands according to the magnitude of the demand. Then we select samples from the same band for different D i 's to construct a sample for D. [40] . Now let us focus on how to generate samples for each D i . The pseudo code for generating a sample for D i is given in Algorithm 1. Note that by repeatedly call this algorithm, one can generate any given number of samples for D i . A pseudo-code for generating a sample for D is given as Algorithm 2.
Algorithm 1 An algorithm for generating a mixed-defective sample for D i (i ∈ S) following AICM.
1: Input data: score vector samples Y i = { ξ k } N k=1 that have been collected from the survey at the customer site i. 2: Determine the rank of ICM based on the score vectors in Y i . Sort the indices in F i (from low to high) based on this rank. 3: Based on this rank, partition Y i into qualified subset Y q i of samples and defective subset (24) to get an optimal estimation of probability parameters denoted as p * , π * and substitute p * , π * into (ICM). 5: Generate N i.i.d. score vectors denoted as {ξ k } N k=1 such that each score vector is generated using the way given as follows: with probability ρ sample ξ k following the fitted (ICM), and with probability 1 − ρ sample ξ k uniformly from Y d i .
Algorithm 2 An algorithm for generating a mixed-defective sample for D following AICM. 
An important property of a classic capacitated facility location model is that the objective as a function of location vector is submodular (Section 4 of [41] ). In Section 4.1, we will investigate if the submodularity still holds for models presented in this work. We have shown that the objective function in the uncapacitated case is submodular if the demand-utility consistency condition (Definition 2.1) holds, while for other uncapacitated cases or the capacitated cases the objective function is not submodular in general which is demonstrated using examples. The submodularity results are summarized in Table 1 . In section 4.2, we investigate a new family of valid and facet-defining inequalities of the intersection polytope induced from the demand and capacity constraints for the deterministic problem. incon. cap.
con./incon.
Consider the baseline deterministic model of service center location problem with decision dependent demand (DDSL). This problem can written as
where f (y) is the optimal value of (DDSL) for the given location vector y.
We will show that f (y) is a submodular function for the uncapacitated case which satisfies the condition of demand-utlity consistency. This indicates that (DDSL) in this case is equivalent to a submodular function maximization problem with a knapsack constraint. An immediate consequence is that there exists a (1 − e −1 ) approximation algorithm for solving this problem in polynomial running time [41, 42] .
Proposition 4.1 For the uncapacitated case of (DDSL) that satisfies the condition of demand-utlity consistency, the objective function f : {0, 1} |F | → R + is submodular.
Proof For a given y, the value f (y) is given by the optimal objective of the following linear program:
With the demand-utility consistency condition, the optimal objective has an analytical form
Since summation preserves submodularity, it suffices to show each h i (y) is submodular. Consider the function h(y) := max j∈[n] v j y j of y ∈ {0, 1} n . Let I(y) = {j ∈ [n] | y j = 1}. Let y 1 and y 2 be any two binary vectors and suppose j * ∈ argmax j∈[n] v j (y 1 ∨ y 2 ) j . If j * ∈ I(y 1 ) ∩ I(y 2 ), then f (y 1 ) + f (y 2 ) = 2v j * y j * = f (y 1 ∨ y 2 ) + f (y 1 ∧ y 2 ). If j * ∈ I(y 1 ) \ I(y 2 ), then f (y 1 ) + f (y 2 ) = v j * y j * + f (y 2 ) ≥ f (y 1 ∨ y 2 ) + f (y 1 ∧ y 2 ). A similar result holds for the case of j * ∈ I(y 2 ) \ I(y 1 ). This shows that h(y) is submodular, and hence f (y).
The following example indicates that the objective function is not submodular in general for uncapacitated case without consistency condition nor the capacitated case. We can get a sense from Example 4.1(b) on why the submodularity is violated for the capacitated problem. Consider the marginal increment of objective brought by F 3 (the third candidate location). Notice that both F 1 and F 3 cannot attract much demand, however F 3 has the largest capacity and the second largest utility value. On the other side, the candidate location F 2 can attract the largest amount of demand (and also with highest utility value) but it has small capacity. So when F 2 is opened, the extra demand attracted by it can be assigned to F 3 . This implies that the marginal increment in objective induced by F 3 is larger in the case that F 2 is opened than F 2 is not opened. A way of restoring submodularity for our model is to formulate a further specification of customers' profile of preference. More specifically, for a customer site i we construct a collection of distinct subsets G = {G k } m k=1 where G k ⊆ F i . This induces a partition of customers into m groups, customers in the k-th group are only willing to go to the locations at in G k . But in this case, the problem is equivalent to a classic capacitated facility location problem and the feature of decision-dependent demand is lost.
Understanding the polyhedral structure of the convex hull induced by its feasible set can be of an independent interest in the theory of mixed-integer programming. Notice that the valid and facet defining inequalities corresponding to the knapsack constraint (DDSL.2) and the capacity constraint (DDSL. 3) have been well studied [43] [44] [45] . However, the demand attraction constraints (DDSL.4) and (DDSL.5) introduced in this research is a signature set of novel constraints that brings new features to the polyhedral structure when competing with the capacity constraints. We focus on investigating the polyhedral structure induced by the demand attraction constraints and the capacity constraints. Specifically, the polytope of interest is defined as follows:
For simplicity in the polyhedral structure analysis, we assume F i = F for all i ∈ S. To motivate insightful valid inequalities, we consider the quantity i∈I j∈J x ij for any subsets I ⊆ S and J ⊆ F. The capacity and demand attraction constraints imply that i∈I j∈J x ij ≤ j∈J C j y j and i∈I j∈J x ij ≤ i∈I max j∈F D j y j . The upper bound of the quantity i∈I j∈J x ij depends on which constraint is more restrictive. Figure 1 illustrates a critical situation that if any one of the service centers in J is not opened the total attracted demand is greater than the total capacity, while if all service centers in J are opened, the total attracted demand is less than the total capacity. Based on this critical situation, we can derive a facet-defining valid inequality of Z(I, J), where Z(I, J) is the convex hull of [y, x, q] with indices restricted on the subsets I and J (see (34) for the formal definition). This inequality can also be facet-defining for Z under some mild conditions. The theoretical results on the polyhedral properties of Z are organized as follows: Theorem 4.1 provides some mild conditions such that a facet-defining inequality of Z(I, J) is also a facet-defining inequality of Z (trivial lifting of facets from a lower dimensional space to a higher one). The proof of Theo- 
We provide a sequence of technical lemmas that will be used to prove the facet-extension results given in Theorem 4.1. Fig. 1 : An illustration of the competition between the capacity constraints and the demand attraction constraints. In all subfigures, three customer sites are represented as circles which are denoted as S1, S2, S3 (from top to bottom) and two candidate service center locations are represented as squares which are denoted as F1, F2 (from top to bottom). The number in a customer site (circle) indicates the attracted demand, and the number in a service center (square) indicates the capacity. In case (a) only F1 is opened which attracts 2 units of demand from S1 and S2 respectively and 1 unit from S3. In this case the total attracted demand is greater than the total capacity. In case (b) only F2 is opened and the situation is similar to (a). In case (c) two service centers are opened, which attracts 2 units of demand from all customer sites. The total attracted demand is less than the total capacity.
Proof Suppose S = {v k } n k=0 are n+1 affinely independent vectors. If span(S) ≥ n, then there must exist n linearly independent vectors in S. We now suppose that span(S) ≤ n − 1. Specifically, suppose span(S) = m, with m ≤ n − 1. Then there exist m linearly independent vectors in S. We assume that {v k } n k=n−m+1 are linearly independent without loss of generality. It follows that v i ∈ span{v k } n k=n−m+1 , ∀i ∈ {0, 1, . . . , n − m}. Therefore, we have span{v k − v 0 } n k=1 ⊆ span{v k } n k=n−m+1 , and hence dim{v k − v 0 } n k=1 ≤ m, which contradicts to that {v k − v 0 } n k=1 are linearly independent.
are affinely independent vectors in R n .
Proof Consider any coefficients {a i } k i=2 and {b i } m i=1 that satisfy the following equation
By linearly independence, we have j∈J b j u j = 0. Since
, we must have j∈J b j u j = 0, which leads to a contradiction. Therefore, we have proved that a i = 0 ∀i ∈ {2, . . . , k} and b j = 0 ∀j ∈ [m], which shows that Suppose {u k + v k } n k=1 are n affinely independent vectors. We now show that span{u k + v k : ∀k ∈ [n]} ∩ span{v k : ∀k ∈ [l]} = {0}. Then by Lemma 4.2, we can deduce that {u k + v k } n k=1 and {v k } l k=1 are n + l affinely independent vectors. Assume that span{u k + v k : ∀k ∈ [n]} ∩ span{v k : ∀k ∈ [l]} = {0}. Then there exist coefficients {α k } n k=1 and {β k } l k=1 which are not all zero, such that
holds by the given hypothesis. It follows that n k=1 α k v k = 0 and hence α k = 0 ∀k ∈ [n] and β k = 0 ∀k ∈ [l], which leads to a contradiction. Proof Part (a) can be easily verified. We now focus on proving Part (b). Without loss of generality, we prove for the case that n l ≥ 4 for all l ∈ [k]. It is easy to verify that all vectors in ∪ k l=1 U l are linearly independent, and by Part (a), the vectors in V are affinely independent. We claim that span(V ) ∩ span{∪ k l=1 U l } = {0 N }. If the claim holds then by Lemma 4.2, the N vectors in ∪ k l=1 U l ∪ W are affinely independent. To prove this claim, we label the vectors in U l as u li for l ∈ [k] i ∈ [n l − 1], and label the vectors in V as v l for l ∈ [k]. We can verify that {v l } k−1 l=1 are linearly independent and v k = − 
Consider the e ns i entry of w for s ∈ [k] and i ∈ {2, . . . , n l − 1}. Evaluating this coefficient based on l∈[k−1] β l v l implies that the coefficient should be 0. On the other side, evaluating this coefficient based on l∈[k] i∈[n l −1] α li u li implies that it should be equal to α s,i − α s,i−1 . Therefore, we must have α s,i = α s,i−1 for s ∈ [k] and i ∈ {2, . . . , n l − 1}. For l ∈ [k], let a l = α l,i for i ∈ {2, . . . , n l − 1}. By matching the coefficient corresponding to the e n1 1 entry of w, we get a 1 = β 1 . By matching the coefficient corresponding to the e n l n l entry of w, we get a l = 0 for l ∈ [k]. This shows that w = 0, which leads to a contradiction.
Lemma 4.5 Let l, m, n be three positive integers. Let u k ∈ R l ⊗ R n and v k ∈ R m ⊗ R n for all k ∈ [n] be arbitrary vectors with dimension l × n and m × n, respectively. Let u ∈ R l ⊗ R n , and v ∈ R m ⊗ R n be two arbitrary vectors. Let {w i } n i=1 be n linearly independent vectors in R n with w i j ∈ {0, 1} ∀i, j ∈ [n]. Let J i = {j ∈ [n] : w i j = 1}. Consider the following sets of vectors in the R l×n ⊕ R m×n ⊕ R n space:
with the convention n + 1 → 1, 
are n linearly independent binary vectors in R n , for each j ∈ [n], there exists an i ∈ [n] (depending on j) such that j ∈ J i . Therefore, there exists a mapping π : [n] → [n] such that π is surjective and j ∈ J π(j) . We define the following sets of vectors:
Clearly, the sets 
Since we have U * ∪V * ∪Y * ∪W ⊆ S, it implies that S should contain (l+m+1)n affinely independent vectors. 
with coefficients α j , β ij and γ is a facet defining valid inequality of Z(I, J) and it is valid for Z. If there exists a set V of (2|I| + 1)|J| affinely independent points in Z(I, J) that are on the plane defined by this inequality and V satisfies the following conditions:
(1) For any point [y, x, q] ∈ V , it holds that y j ∈ {0, 1} ∀j ∈ J;
(2) ∀j ∈ J there exists an point [y (j) , x (j) , q (j) ] ∈ V such that y (j) j = 1 and i∈I x (j) ij < C j (these points don't have to be different with respect j); (3) For each i ∈ I, there exists a point [ỹ (i) ,x (i) ,q (i) ] ∈ V such that j∈J q (i) ij < 1 (these points don't have to be different with respect i).
Then the inequality (34) defines a facet of Z.
Proof We first show that Z is of full dimension. Since Z is defined by (2|S| + 1)|F| variables, we will prove dim Z = (2|S| + 1)|F| by constructing (2|S| + 1)|F| + 1 affinely independent points in Z as follows:
[0, 0, 0], [e j , 0, 0] ∀j ∈ F, [e j , 0, e ij ] ∀i ∈ S ∀j ∈ F, [e j , e ij , e ij ] ∀i ∈ S ∀j ∈ F, where is a sufficient small positive constant. In particular, < min{D ij , C j }. It is easy to verify the above (2|S| + 1)|F | + 1 points are affinely independent points in Z, and hence dim Z = (2|S|+1)|F |. Using a similar argument, we can show that dim Z(I, J) = (2|I| + 1)|J|. We now focus on proving the theorem.
Part (a) We refer the inequality j∈J α j y j + i∈I j∈J (x ij + β ij q ij ) ≤ γ as ineq in the proof. Since ineq defines a facet of Z(I, J), there exist (2|I|+1)|J| affinely independent points in Z(I, J) that are on the plane defined by ineq. We denote this set of the (2|I|+1)|J| affinely independent points as V . We will construct additional (2|S| + 1)|F| − (2|I| + 1)|J| linearly independent points of Z that are on the plane defined by ineq. Let [y, x, q] be an arbitrary point in V , we first construct the following set of points:
with > 0 sufficient small. We can verify that all points in V 1 are linearly independent points of Z that are on the plane defined by ineq, and |V 1 | = (2|S|−2|I|+1)(|F|−|J|). Let [y (j) , x (j) , q (j) ] ∈ V for every j ∈ J be the vector satisfying the condition (2) of (a). We construct the following set of points:
with sufficient small > 0. One can verify that all points in V 2 are linearly independent points of Z that are on the plane defined by ineq, and V 2 contains 2(|S| − |I|)|J| linearly independent points. Finally, we construct the following set of points: (29) . Suppose F i = F for all i ∈ S and D ij > 0 for all i ∈ S, j ∈ F. The following hold:
(a) For any subsets of indices I ⊆ S, J ⊂ F and index j ∈ F \ J the following inequality is valid for Z i∈I j∈J∪{j }
(b) Let I ⊆ S and J ⊂ F be two subset of indices, if j∈J C j ≤ i∈I max j∈J D ij , then the following inequality is valid for Z i∈I j∈F
Proof Part (a): It suffices to show the validness for the two cases: y j = 0 and y j = 1. If y j = 0, then we have
where we use y j = 0, y j ≤ 1 and the constraint (DDSL.3) in deriving (38) .
For the case of y j = 1, we have
where we use y j = 1, and the constraint (DDSL.4) in deriving (39) . Part (b): For a feasible point [y, x, q], let J = {j ∈ F \ J | y j = 1}. For any i ∈ I, let π(i) = argmax j∈J∪J D ij where ties are broken arbitrarily. Let I j = {i ∈ I | π(i) = j} for every j ∈ J ∪ J , and by definition the subsets {I j } j∈J∪J form a partition of I. If |J | ≤ 1, then (37) reduces to (36) . Therefore without loss of generality, we assume |J | ≥ 2. The expression on the left can be relaxed as i∈I j∈F
Notice that we have the following inequalities hold i∈I j∈J∪J
Furthermore, the second term in (37) is evaluated as
(42) The first term in (42) can be decomposed as
where T 1 and T 2 represent the first and second term on the right side of (41), respectively. The term T 1 is evaluated as
where we use the fact that for any j * ∈ J, i ∈ I j * and j ∈ F \ J, we have D ij * = max j∈J∪{j } D ij = max j∈J D ij . The term T 2 can be lower bounded as:
Substituting (44) and (45) into (43) yields the following inequalities:
where we use the fact that {I j * | j * ∈ J ∪ J } form a partition of I. Substituting (46) into (42) yields
where we use the assumption i∈I max j∈J D ij − j∈J C j ≥ 0. Substituting (41) and (47) 
i∈I j∈F
This shows that (37) is valid for Z. 
Then the following inequality defines a facet of Z i∈I j∈J
Proof We will first provide an intuition of how the coefficients of the inequality (50) are determined given the condition in (49) . Since we restrict the decision variables induced by the subsets I and J, a facet-defining valid inequality is likely in the form of i∈I j∈J
We need to determine the coefficients a j and b. Notice that there are |J| + 1 coefficients to be determined, and we will set up |J| + 1 linear equations corresponding to the following |J| + 1 situations: (1) opening service centers at all locations in the subset J, (2) opening service centers at all locations in J \ {j} for all j ∈ J. Notice that (2) contains |J| different situations referred as situation (2.j). In the first situation, all locations in J have service centers. So the maximum value of i∈I j∈J x ij is i∈I max j∈J D ij , i.e., the aggregated demand constraint is binding. Similarly, in the situation (2.j 0 ), locations in J \ {j 0 } have service centers, and the maximum value of i∈I j∈J x ij is j∈J\{j0} C j , i.e., the aggregated capacity constraint is binding. These observations help us set up the following linear equations to determine a j and b:
By solving the above linear equation system, we can determine the values of all coefficients:
which leads to the inequality (50). Now we give a formal proof of that (50) is a facet-defining valid inequality of Z. To see the validness, we divide the analysis into four cases: (a) y j = 1 for all j ∈ J; (b) y j = 1 and y j = 0 for all j ∈ J \ {j}; (c) there exists a J 0 ⊂ J satisfying 2 ≤ |J 0 | ≤ |J| − 1 and y j = 1 for all j ∈ J \ J 0 and y j = 0 for all j ∈ J 0 ; (d) y j = 0 for all j ∈ J. The validness of (50) clearly holds in the cases (a) and (b) by the way of determining the coefficients. For the case (d), (50) becomes i∈I j∈J
It is valid because in this case we have i∈I j∈J x ij = 0 and C j − i∈I max j ∈J D ij > 0 by assumption. It remains to verify the validness in the case (c). Substitute y j = 1 for j ∈ J \ J 0 into (50), and notice that we have
where the assumed condition is used to get the first inequality above. This proves the validness in this case.
To show that (50) defines a facet of Z, we will first show that it defines a facet of Z(I, J) (see (34) ), and apply Theorem 4.1 to conclude that it also defines a facet of Z via verifying the three conditions in 
It suffices to show that the above coefficients are different from the coefficients in (50) by a positive constant. Note that all points in V should satisfy (53) as an equality. We are going to propose some specific points from V and plug them into (53) to determine the coefficients. First, we notice that there must exist a vector [ j∈J e j , x 0 , q 0 ] ∈ V such that x 0 ij > 0, 0 < q 0 ij < 1 for all (i, j) ∈ I × J, and i∈I x 0 ij < C j for all j ∈ J. It can be deduced that the points [ j∈J e j , x 0 + e ij1 − e ij2 , q 0 ] for all i ∈ I, j 1 , j 2 ∈ J, j 1 = j 2 , and [ j∈J e j , x 0 + e i1j − e i2j , q 0 ] for all i 1 , i 2 ∈ I, j ∈ J are also in V for some sufficient small . Substituting these points into (53) as an equality, we conclude that α ij1 = α ij2 and α i1j = α i2j , which implies that α ij = α ∀(i, j) ∈ I × J for some α > 0. Without loss of generality, we can set α = 1 since the inequality is invariant under multiplication of a positive constant. Similarly, the exists a point [ j∈J\{j0} e j , x 1 , q 1 ] ∈ V such that 0 < q 1 ij < 1 for all (i, j) ∈ I × (J \ {j 0 }). It can be deduced that the points [ j∈J\{j0} e j , x 1 , q 1 + e i1j1 ] ∀(i 1 , j 1 ) ∈ I × (J \ {j 0 }) are all in V . Substituting this point into (53) as an equality, we conclude that γ ij = 0 for all (i, j) ∈ I × (J \ {j 0 }). Since j 0 is arbitrary, we conclude that γ ij = 0 for all (i, j) ∈ I × J. Therefore, (53) becomes i∈I j∈J
Clearly, since the previously considered point [ j∈J e j , x 0 , q 0 ] ∈ V should satisfy i∈I j∈J x 0 j = i∈I max j∈J D ij (this can be obtained by substituting the point into (50) as an equality). Then substituting this point into (54) as an equality gives an equation in the same form of (52a). Similarly, substituting the point [ j∈J\{j0} e j , x 1 , q 1 ] ∈ V into (54) as an equality gives an equation in the same form of (52b). Therefore the to be determined coefficients [β, γ] are the solution of (52), which shows that the facet defining valid inequality (54) of Z(I, J) is exactly (50). Now we show that (50) is also a facet defining inequality of Z using Theorem 4.1. To achieve this, we first need to show that (50) is valid for Z. The validness of (50) for Z is based on the following two observations: if a subset J ⊆ F \ J of service centers are opened, the maximum attracted demand from a customer site i is at most max j∈J∪J D ij ≤ max j∈F D ij . But max j∈F D ij = max j∈J D ij by assumption, which implies that once all places in J have service centers, the aggregated demand constraint i∈I j∈J x ij ≤ i∈I max j∈J D ij is valid and it is independent of whether or not other places having service centers. Similarly, in the case that places at J \{j 0 } have service centers, the aggregated capacity constraint dominates the aggregated demand constraint which is also independent of service center locations at other places due to the fact max j∈J\{j0} D ij ≤ max j∈J∪J \{j0} D ij for any J ⊆ F \ J. It remains to show that there exist (2|I| + 1)|J| affinely independent vectors in Z(I, J) that satisfy the three conditions in Theorem 4.1. Since we have already proved that (50) is a facet defining inequality of Z(I, J), there must exist a set V (2|I| + 1)|J| affinely independent vectors in Z(I, J) that satisfy (50) as an equality. Furthermore, we can make sure that V includes the following two vectors: [ j∈J e j , x 0 , q 0 ] and [ j∈J\{j0} e j , x 1 , q 1 ] where the two vectors are introduced in the previous analysis of this proof. Notice that [ j∈J e j , x 0 , q 0 ] satisfies the condition (2) of Theorem 4.1 for all j ∈ J, and [ j∈J\{j0} e j , x 1 , q 1 ] satisfies the condition (3) of Theorem 4.1 for each i ∈ I. This concludes the proof.
The numerical investigation consists of two major parts: the numerical analysis of the AICM and the computational performance of the (DRO-FL) model. In Section 5.1, we provide numerical analysis of the AICM, where we simulate certain amount of score vector samples and use them to fit the AICM. In Section 5.2, we show the computational performance of solving generated instances of the (DRO-FL) model, and provide sensitivity analysis of the robust optimal locations with respect to the sample size and the choice of radius in the ambiguity set.
For the numerical analysis of AICM, we focus on a specific customer site with 4 candidate service center locations labeled as L 1 ∼ L 4 (the rank of attractability is unknown). We consider 6 possible scores {0, 1, . . . , 5}. We simulate N i.i.d. score vectors with some defected ones and use them to fit the AICM. The generation of each score vector is based on the following procedures. First, we simulate a 4-dimensional real vector v following the Gaussian distribution N (w, Σ) with w = [1, 3, 4.5, 2] corresponding to the mean score of L 1 ∼ L 4 respectively, and Σ = diag[0.8 2 , 1.5 2 , 1, 1.5 2 ] as the covariance matrix. Next, each entry of v is rounded to the closest integer from {0, 1, . . . , 5}. Finally, with probability 0.8 we accept this score vector, while with probability 0.2 we randomly select an entry of v, set it to be zero and accept the resulting score vector. We follow the steps instructed in Section 3.2 to fit the AICM. An essential step is to solve the nonlinear optimization problem (24) . This problem is implemented in Julia 1.5.2 [46] with optimization package JuMP dev [47] as the modeling interface and Ipopt 3.13.2 [48] as the solver. For simplicity, we let λ 1 = λ 2 = λ in (24) . We use a 4-fold cross validation approach to select an appropriate regularization parameter λ. The parameter selection consists of two rounds. The first round is a rough selection, where the candidate λ (1) is selected from {0, 0.1, 0.2,. . . , 5} based on the out-of-sample performance. Suppose λ (1) * is the selected value after the first round. In the second round, we select an optimal λ * from the set λ (1) * ± 0.01k k = 0, 1, . . . , 10 based on the out-of-sample performance. The out-of-sample performance is measured by the averaged logarithm of likelihood. The first-round selection results are given in Figure 2 for three different sample sizes N = 50, 100, 200. The second-round selection results are given in Table 3 , from which we see that the optimal λ is 2.98, 2.98 and 0.37 for N = 50, 100, 200, respectively. The curves in Figure 2 indicate that for a small sample size (N = 50) a larger regularization parameter lead to a better out-of-sample performance, while for a larger sample size (N = 200) the best out-of-sample performance is achieved at a smaller regularization parameter. This behavior is intuitive since training the model with a small sample size can lead to over fitting if the model parameters are not sufficiently regularized. When more samples are available to train the model, the regularization can be less restrictive.
Finally, we fit the AICM using N samples and the optimal λ parameter. The probability parameters p and π from AICM can be used to estimate the probability distribution on scores for each candidate location. The results are given in Figure 3 for the 4 candidate locations. The probability distribution on score values reveal that the rank of attractability from high to low is: Location 3 > Location 2 > Location 4 > Location 1, which matches with the underline stochastic model that generates the scores. The mean scores corresponding to the four locations are given in Table 2 . (a) Probability distribution on scores of location 1.
(b) Probability distribution on scores of location 2.
(c) Probability distribution on scores of location 3.
(d) Probability distribution on scores of location 4. 
We investigate the computational performance of solving numerical instances of the (DRO-FL) model. To generate numerical instances, we create a 100×100 two-dimensional square as the map. The customer sties are points inside the square, and the coordinates of each customer site are uniformly drawn from the range [0, 100] 2 . We let F = S, which means each customer site can be a candidate location of a service center. We sort the set of distances W = {d(i, j) | ∀i ∈ S, ∀j ∈ S, i < j} from small to large, and let d c be the 5% quantile of the sorted list. For every customer site i, we select all customer sites (including i) that are within d c distance from i to form the set F i . We let the utility completely depend on the distance. Specifically, we set u ii = 5 for all i ∈ S. We also let the first element in W correspond to utility 5, the last element in W correspond to utility 0.5, and any other elements correspond to utility values obtained using a linear interpolation. For a sample ω, the demand D ω ij is drawn from the normal distribution N (120u ij , (12u ij ) 2 ). The capacity of each candidate service center is set to be 1000. We assume the cost of opening each service center is identical (an unit cost), and the budget B in this case is equivalent to the maximum number of service centers that can be opened. We have generated 12 instances of (DRO-FL), and each instance is characterized by the number |S| of customer sites (also candidate locations) and the budget B (number of maximum service centers that can be allocated). For each instance, we have tried three different sample-size options (N = 100, 200, 500), and the total variation distance is set to be 0.2 in these calculation. The (DRO-FL) model is implemented and solved using a Gurobi Python interface (Python 3.7 and Gurobi 9.0). For every instance, the master problem and all scenario sub-problems are solved using a single 2.50GHz CPU. The computational results for the 12 instances are shown in Table 4 . It is shown that all instances have been solved to optimality within 120 seconds, and within 12 master iterations (see Algorithm 3). The range of the robust optimal objective is less than 0.7% across three different sample-size options for each instance, and we have also observed that the robust optimal solution remains the same for different sample-size options.
We further investigate the impact of the total variation distance d on the robust optimal decision and the objective value. In this study, we focus on the two instances (|S|, B) = (60, 20) and (|S|, B) = (60, 10). For a given sample size N (with N = 100, 200, 500), we tune the total variation distance d from 0 to 0.4 with an increment 0.05 and investigate the change of optimal objective of (DRO-FL). The results are shown in Figure 4 . Obviously, the objective should decrease with d as a larger d value can allow the worse scenarios get more weight and hence makes the objective decrease. With d increasing from 0 to 0.4, the decrement of the objective is 1.5% (1000 units) for the case B = 20, and 1.9% (700 units) for the case B = 10. The relative decrement is about one magnitude smaller compared to the 10% relative deviation in generating samples for the demand. It is also observed that the optimal decision of (DRO-FL) does not change with the sample size or the total variation distance. Table 4 : Computational performance of solving 12 numerical instances of (DRO-FL). In all instances, the total variation distance for the ambiguity set is set to be 0.2. All instances are solved to optimality. The impact of total variation distance on the robust optimal objective of (DRO-FL).
The distributionally-robust service center location problem investigated in this research possesses decision-dependent demand induced naturally from a maximum attraction principle and the number of opened service centers in a neighborhood. The ambiguity set considered in this work is a decision-independent one defined using the total-variational distance while many other metics such as the Wasserstein distance, Phi-divergences, etc. can also be applied to define the ambiguity set with minor modification in the formulation. The modeling approach in this work decouples the endogenous impact (decision dependency) and distributional ambiguity. As a result, the model is highly computationally efficient for mid-and large scale instances compared to models with an endogenous uncertainty (ambiguity) set. In practice, such decoupling is likely to be more data driven, as there is a higher chance to identify data-based evidence that supports decision-dependent demand than decision-dependent uncertainty of demand.