Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
108,733 characters · 21 sections · 44 citation commands
Competitive Facility Location with Market Expansion and Customer-centric Objective
\newcolumntype{C}{>{\arraybackslash}p{4em}}
\affil[1]{\it School of Computing and Information Systems, Singapore Management University} \affil[2]{\it Interdisciplinary Centre for Security, Reliability and Trust, University of Luxembourg} \affil[3]{\it SLSCM and CADA, Faculty of Data Science and Artificial Intelligence, College of Technology, National Economics University, Hanoi, Vietnam } \affil[*]{\it Corresponding author, [email removed]}
{\bf Keywords:} Competitive facility location, market expansion, customer satisfaction, inner-approximation.
The problem of facility location has been a key area of focus in decision-making for modern transportation and logistics systems for a long time. Typically, it involves selecting a subset of potential sites from a pool of candidates and determining the financial investment required to establish facilities at these chosen locations. The goal is usually to either maximize profits (such as projected customer demand or revenue) or minimize costs (such as operational or transportation expenses). A critical factor in these decisions is customer demand, which significantly influences facility location strategies. In this study, we focus on a specific class of competitive facility location problems, where customer demand is predicted using a random utility maximization (RUM) discrete choice model McFaTrai00 and the aim is to locate new facilities in a market already occupied by a competitor train2009discrete,BenaHans02,mai2020multicut. Here, it is assumed that customers choose between available facilities by maximizing the utility they derive from each option. These utilities are typically based on attributes of the facilities, such as service quality, infrastructure, or transportation costs, as well as customer characteristics like age, income, and gender. The use of the RUM framework in this context is well supported, given its widespread success in modeling and predicting human choice behavior in transportation-related applications Mcfadden2001economicNobel,BenABier99a.
In the context of competitive facility location under the RUM framework, to the best of our knowledge, existing studies generally assume that the maximum customer demand that can be captured by each facility is fixed and independent of the availability of new facilities entering the market. However, this assumption is limited in many practical scenarios. Intuitively, the total market demand is likely to expand when more facilities are built. Moreover, most of the existing studies focus solely on maximizing the total expected captured demand, ignoring factors that account for customer satisfaction.
An example that highlights the importance of considering such factors is when a new electrical vehicles (EV) company plans to build electric vehicle charging stations to compete with other competitors (such as gas stations or public transport). A critical consideration for the firm is that adding more EV stations in the market could likely expand the EV market, attracting more customers from competitors Sierzchula2014,Li2017. Additionally, in certain cases, building more EV stations in urban areas might help generate more profit by attracting more customers, but this may not be the best long-term strategy. Customers from non-urban areas would have less access to these facilities and may lose interest in adopting EVs, which may hinder broader EV adoption Gnann2018,Bonges2016. Thus, for a long-term, sustainable development strategy, the company would need to balance overall profit with customer satisfaction.
Motivated by this observation, in this paper, we explore two new considerations that better capture realistic customer demand and balance both the company’s profit and customer satisfaction. Specifically, we assume that the maximum customer demand (i.e., the total number of customers that existing and new facilities can attract) is no longer fixed but modeled as an increasing market expansion function of customer utility. We also introduce a term representing total customer utility to account for customer satisfaction in the main objective function. The resulting optimization problem is highly non-convex, and to the best of our knowledge, no existing algorithm can solve it to optimality, or guarantee near-optimal solutions. To address these challenges, we have developed innovative solution algorithms with theoretical support that can guarantee near-optimality under both concave and non-concave market expansion functions. Our key contributions are detailed as follows:
\paragraph{Paper Outline:} The paper is structured as follows: Section (ref) introduces the problem formulation. Section (ref) discusses the submodularity of the objective function in the context of concave market expansion functions. In Section (ref), we present our inner-approximation solution method. Section (ref) addresses our approaches for the facility location problem with a general non-concave market expansion function. Section (ref) presents the numerical results, while Section (ref) concludes the paper. Additional proofs and further details not covered in the main body are provided in the appendix.
Notation: Boldface characters represent matrices (or vectors), and $a_i$ denotes the $i$-th element of vector $\textbf{a}$. We use $[m]$, for any $m\in \mathbb{N}$, to denote the set $\{1,\ldots,m\}$.
Competitive facility location under random utility maximization (RUM) models has been a topic of interest in Operations Research and Operations Management for several decades. This area of research differentiates itself from other facility location problems through the use of discrete choice models to predict customer demand, drawing from a well-established body of work on discrete choice modeling train2009discrete. In the context of competitive facility location (CFL) under RUM models, most studies adopt the Multinomial Logit (MNL) model to represent customer demand. Notably, BenaHans02 were among the first to introduce the CFL problem under the MNL model, utilizing a Mixed-Integer Linear Programming (MILP) approach that combines a branch-and-bound procedure for small instances with a simple variable neighborhood search for larger instances.
Subsequent contributions include alternative MILP models proposed by Zhang2012 and Haase2009. haase2014comparison conducted a benchmarking study of these MILP models, concluding that Haase2009's formulation exhibited the best performance. Freire2015 enhanced Haase2009's MILP model by incorporating tighter inequalities into a branch-and-bound algorithm. Additionally, Ljubic2018outer developed a Branch-and-Cut method that combines outer-approximation and submodular cuts, while mai2020multicut introduced a multicut outer-approximation algorithm designed for efficiently solving large instances. This method generates outer-approximation cuts for groups of demand points rather than for individual points.
A few studies have also explored CFL using more general choice models, such as the Mixed Multinomial Logit (MMNL) model Haase2009,haase2014comparison. However, applying the MMNL model typically requires large sample sizes to approximate the objective function, leading to complex problem instances. dam2022submodularity,dam2023robust incorporated the Generalized Extreme Value (GEV) family into CFL and proposed a heuristic method that outperforms existing exact methods. mendez2023follower investigated CFL under the Nested Logit (NL) model, proposing exact methods based on outer-approximation and submodular cuts within a Branch-and-Cut procedure. Recently, le2024competitive explored CFL under the Cross-Nested Logit model, considered one of the most flexible discrete choice models in the literature. In their work, the authors demonstrated that, although the objective function is not concave, it can be reformulated as a mixed-integer concave program, allowing the use of standard exact methods like outer-approximation.
In all the aforementioned studies, the market size is assumed to be fixed and independent of the customer’s total utility. Additionally, these works focus solely on maximizing expected captured demand, neglecting factors related to customer satisfaction. On the other hand, because the objective function in most cases can either be shown to be concave or reformulated as a concave program, outer-approximation methods mai2020multicut,duran1986outer have remained the state-of-the-art approaches. Our work, therefore, makes a significant advancement in this literature by introducing a novel problem formulation that accounts for both market dynamics and customer satisfaction. Furthermore, we propose a new near-exact approach based on inner-approximation, which guarantees smaller approximation errors compared to traditional outer-approximation methods.
Our work and the general context of choice-based competitive facility location are related to a body of research on competitive facility location where customer behavior is modeled using gravity models DREZNER2002, ABOOLIAN2007a, ABOOLIAN2007b, ABOOLIAN2021, lin2022locating. These models, in their classical form without market expansion and customer objective components, share a similar objective structure with the CFL problem under the MNL model. Market expansion perspectives have also been considered in this line of work ABOOLIAN2007a, ABOOLIAN2007b, lin2022locating. However, since these studies rely on different customer behavior assumptions, the form of the customer's total utility significantly differs from the total utility function under the discrete choice models considered in our work. Moreover, while these works are restricted to concave market expansion functions, our work considers both concave and non-concave functions, allowing broader applications. In terms of methodological developments, while prior work employs outer-approximation approaches to handle the nonlinear concave demand function, we explore a new type of approximation based on “inner-approximation”. This approach not only offers smaller approximation gaps but also allows efficient solving of problems with general non-concave market expansion functions.
In the classic facility location, decision-makers aim to establish new facilities in a manner that optimizes the demand fulfilled from customers. However, accurately assessing customer demand in real-world scenarios is challenging and inherently uncertain. In this study, we study a facility location problem where discrete choice models are used to estimate and predict customer demand. Among various approaches discussed in demand modeling literature, the Random Utility Maximization (RUM) framework train2009discrete stands out as the most prevalent method for modeling discrete choice behaviors. This method is grounded in the random utility theory, positing that a decision-maker's preference for an option is represented through a random utility. Consequently, the customer tends to opt for the alternative offering the highest utility. According to the RUM framework McFa78,FosgBier09, the likelihood of individual $n$ choosing option $i\in S$ is determined by $P(u_{ni}\geq u_{nj},;\forall j\in S)$, implying that the individual will select the option providing the highest utility. Here, the random utilities are typical defined as $u_{ni}=v_{ni} + \epsilon_{ni}$, where $v_{ni}$ represents the deterministic component, which can be calculated based on the characteristics of the alternative and/or the decision-maker and some parameters to be estimated, and $\epsilon_{ni}$ represent random components that are unknown to the analyst. Under the popular Multinomial Logit (MNL) model, the probability that a facility located at position $i$ is chosen by an individual $n$ is computed as $P_n(i|S) = \frac{e^{v_{ni}}}{\sum_{i\in S}e^{v_{ni}}}$, where $S$ is the set of available facilities.
In this study, we consider a competitive facility location problem where a “newcomer” company plans to enter a market already captured by a competitor (e.g., an electrical vehicle (EV) company is aiming to break into the transportation market, which is currently dominated by companies offering gasoline-powered vehicles or other EV brands.). The main objective is to secure a portion of the market share by attracting customers to their newly opened facilities. To forecast the impact of these new facilities on customer demand, we employ the RUM framework, which assumes that each customer assigns a random utility to each facility (both the newcomer's and competitors') and makes decisions aimed at maximizing their personal utility. Consequently, the company's strategy revolves around selecting an optimal set of locations for its new facilities to maximize the anticipated customer footfall.
To describe the mathematical formulation of the problem, let $[m]$ be the set of available locations, $[N]$ be the set of customer types available in the market, whereas a customer's type can be defined by geographic locations. Moreover, let $v_{ni}$ be the utility of facility located at location $i\in [m]$ associated with customer type $n\in [N]$, and ${\mathcal{S}}^c$ be the set of competitor's facilities. We also denote $q_n$ be the {maximum customer expenditure} in zone $n\in [N]$. Given a location decision $S\subseteq [m]$, i.e., set of chosen locations and under the MNL choice model, the choice probability of a new facility $i\in [m]$ is given as: \[ P(i\Big|S \cup {\mathcal{S}}^c ) = \frac{e^{v_{ni}}}{\sum_{j\in S} e^{v_{nj}} +\sum_{j\in {\mathcal{S}}^c} e^{v_{nj}}}. \] The competitive facility location problem, in its classical form, can be formulated as:
The above formulation has been widely employed in the context of choice-based facility location BenaHans02,FLO_Hasse2009MIP,Ljubic2018outer,mai2020multicut. This formulation, however, presumes that the total demand for customer type $n$ (that is, $q_n$) remains constant, regardless of an increase in demand as more facilities become available in the market. Additionally, this formulation does not consider customer satisfaction, which is likely to enhance with the availability of more facilities in the market. To address these shortcomings, let us consider the following customer's expected utility function as a function of the chosen locations $S$, under the assumption that customers make choices according to the MNL model train2009discrete:
The function \(\phi_n(S)\) represents the expected utility experienced by customers of type \(n\) when the available facilities in the market are those in the set \(S \cup \mathcal{S}^c\). This function is commonly referred to as the expected consumer surplus associated with the choice set \(S \cup \mathcal{S}^c\). It captures the inclusive value of the choice set of available facilities, reflecting the combined attractiveness of all available alternatives within it train2009discrete,daly1978logsum.
It is to be expected that the total demand of customers would be a increasing function of $\phi_n(n)$, since an increase in customer utilities should be likely to attract more customers to the market. With this consideration, we introduce the following formulation that enable us to capture both market expansion and customer-centric values in the objective function.
where $g(t)$ is an increasing function that reflects the impact of customers' expected utilities on market expansion (namely, customers' expenditures), and $\alpha_n$ represent specific parameters. These parameters, $\alpha_n$, are scalar values that help quantify the balance between the firm's captured demand and the expected utility for customers. An increase in $\alpha_n$ would enhance customer satisfaction, but might negatively influence the firm's captured demand, and vice versa. Furthermore, a location solution $S$ that boosts the customer's expected utility $\phi_n(S)$ will also attract more customers, thereby expanding the overall market via the increasing function $g_n(\cdot)$. For notational simplicity, we include only a basic cardinality constraint on the number of open facilities, \( |S| \leq C \), while noting that our approach is general and capable of handling any linear constraints.
It is convenient to formulate (ref) as a binary program. To simplify notation, let's first denote $V_{ni} = e^{v_{ni}}$ and $U^c_n = \sum_{j\in {\mathcal{S}}^c} V_{nj}$. We then reformulate (ref) as the following nonlinear program:
We refer to the problem as the maximum capture problem with market expansion (ME-MCP). By further letting $z_n = U^c_n + \sum_{i\in [m]} x_iV_{ni}$ \footnote{Previous works typically assume that $U^c_n = 1$ for all $n \in [N]$ for ease of notation, without loss of generality mai2020multicut,dam2022submodularity. This is possible because we can divide both the numerator and denominator of each fraction in (ref) to normalize $U^c_n$ to one. However, this approach is not applicable in our context as it would affect the total expected utility $U^n_c + \sum_{i \in S} V_{ni}$.}. We now write rewrite (ref) in a more compact form as follows:
In the context of choice-based competitive facility location, without the market expansion term $g(\log(z_n))$ and the customer-centric term $\alpha_n \log(z_n)$, existing solutions typically rely on the objective function being concave in $\mathbf{x}$ and submodular, enabling exact solutions via outer-approximation methods, or rapid identification of good solutions with approximation guarantees through the use of greedy location search algorithms Ljubic2018outer,mai2020multicut,Dam2021submodularity. This approach prompts the question of whether such convexity and submodularity properties remain preserved in our new model with the market expansion and customer-centric terms. We will investigate this matter further in the next section.
To effectively and reasonably address market expansion, it is reasonable to assume that the market-expansion function $g(t)$ exhibits an increasing behavior in $t$, as an increase in customers' utilities typically fosters market growth. Additionally, it is essential that $\lim_{t\rightarrow \infty} g(t) = 1$, ensuring that total demand does not surpass the maximum customer expenditure, i.e. $q_n$. Commonly utilized function forms in the literature of market expansion include $g(t) = \frac{t}{t+\alpha}$ and $g(t) = 1 - \alpha e^{-\beta t}$ ABOOLIAN2007a,lin2022locating, both of which exhibit concavity in $t$. Thus, in the subsequent section, our primary focus will be on solving the facility location problem under concave market expansion functions $g(t)$, followed by an exploration for addressing the problem under more general, non-concave market expansion functions.
In this section, we focus on the setting that the market expansion function $g(t)$ is concave, delving into the question of under which conditions the overall objective function is concave and submodular, enabling the use of some efficient outer-approximation and local search algorithms. Specifically, we will first establish conditions for the market expansion function $g(\cdot)$ under which the objective function $F(\textbf{z})$ is concave in $\textbf{z}$. We will further show that under these conditions, the objective function $ {\mathcal{F}}(S)$ (the objective function defined in terms of a subset selection) is monotonically increasing and submodular. As a result, (ref) can be conveniently solved by outer-approximation or local search methods. We further leverage the fact that $F(z)$ is univariate to explore an inner-approximation mechanism which allows us to approximate (ref) by a MILP with arbitrary precision. We will theoretically prove that this inner-approximation approach always yields small approximation errors, as compared to an outer-approximation approach.
From the formulation in (ref), we first consider function $\Psi_z(z_n)$, for any $n\in [N]$, defined as follows: \[ \Psi_n(z_n) = q_n g(\log(z_n)) \left( \frac{ z_n - U^c_n}{z_n}\right) + \alpha_n \log(z_n). \] This is an univariate function of $z_n$, depending on the market expansion function $g(t)$. In the following theorem, we first state conditions under which $\Psi_n(z_n)$ the objective function ${\mathcal{F}}(\textbf{z})$ are concave in $\textbf{z}$.
Given two popular forms $g(t) = \frac{t}{t+\alpha}$ and $g(t) = 1 - \beta e^{-\alpha t}$, for $\alpha, \beta > 0$, Proposition (ref) establishes conditions for $\alpha$ and $\beta$ that ensure $\Psi_n(z_n)$ exhibits concavity with respect to $z_n$.
The proposition can be verified straightforwardly. The concavity of $\Psi_n(z_n)$ implies that the objective function in (ref) is also concave, enabling exact methods such as an outer-approximation algorithm duran1986outer,mai2020multicut to be applied. Second, leveraging the concavity, we can further demonstrate that the objective function of (ref), when defined as a subset function, is submodular. To prove this result, let us consider the objective function defined as a set function in (ref), which can be written as: \[ {\mathcal{F}}(S) = \sum_{n\in [N]} q_n g\left(\log\left(U^c_n+\sum_{i\in S} V_{ni}\right)\right) \left( \frac{\sum_{i\in S} V_{ni}}{U^c_n+ \sum_{i\in S} V_{ni}}\right) + \sum_{n}\alpha_n \log \left(U^c_n+\sum_{i\in S} V_{ni}\right). \] The following theorem demonstrates that the conditions used in Theorem (ref), which ensure that \({\mathcal{F}}(\textbf{x})\) is concave with respect to \(\textbf{x}\), are also sufficient to guarantee that \({\mathcal{F}}(S)\) is submodular.
The proof, which explicitly leverages the concavity of ${\mathcal{F}}(\textbf{x})$ to verify submodularity, is provided in the appendix. A direct consequence of the submodularity shown in Theorem (ref) is that a simple polynomial-time greedy algorithm can always return $(1-1/e) \approx 0.6321$ approximation solutions. Such a greedy algorithm can be executed by starting from the null set and adding locations one at a time, choosing at each step the location that increases ${\mathcal{F}}(S)$ the most. This phase finishes when we reach the maximum capacity, i.e., $|S| = C$. This greedy procedure can run in $(mC \tau)$ time, where $\tau$ is the computing time to evaluate ${\mathcal{F}}(S)$ for a given subset $S\subseteq [m]$. Due to the monotonicity and submodularity, if $\overline{S}$ is a solution returned by the above greedy procedure, then it is guaranteed that ${\mathcal{F}}(\overline{S}) \geq (1-1/e) \max_{S,~ |S|\leq C}{\mathcal{F}}(S)$ Nemhauser1978. We state this result in the following corollary.
In this section, we discuss two methods—exact or near-exact—for solving (ref), taking advantage of the concavity property outlined in Theorem (ref). Specifically, we will briefly introduce the outer-approximation method, widely recognized in the literature for addressing mixed-integer nonlinear programs with convex objectives and constraints duran1986outer,mai2020multicut. Additionally, we explore an approximation approach that allows solving (ref) to near-optimality (with arbitrary precision) by approximating it by a MILP.
The outer-approximation method duran1986outer,mai2020multicut,fletcher1994solving is a well-known approach for solving nonlinear mixed-integer linear programs with convex objective functions and convex constraints. A multi-cut outer-approximation algorithm can execute by building piece-wise linear functions that outer-approximate each nonlinear component of the objective functions (or constraints). In the context of the ME-MCP, this can be done by rewriting (ref) equivalently as:
Since $\Psi_n(z_n)$ is concave in $z_n$, it is well-known that, for any $\overline{z}_n\geq 0$, $\Phi(z_n) \leq \Psi_n(\overline{z}_n) + \Psi'_n(\overline{z}_n)(z_n - \overline{z}_n)$, for any $z_n >0$. This implies that, for any $\overline{z}_n >0$, the following inequality is valid for the ME-MCP: $\theta_n \leq \Psi_n(\overline{z}_n) + \Psi'_n(\overline{z}_n)(z_n - \overline{z}_n)$. Such valid inequalities are typically refereed to as outer-approximation cuts. It then follows that one can replace the nonlinear constraints (ref) by sub-gradient cuts $\theta_n \leq \Psi_n(\overline{z}_n) + \Psi'_n(\overline{z}_n)(z_n - \overline{z}_n)$ for all $\overline{z}_n$ in the feasible set. The multi-cut outer-approximation is an iterative cutting-plane procedure where, at each iteration, a master problem is solved, with (ref) being replaced by linear cuts. After each iteration, let $(\overline{\pmb{\theta}},\overline{\textbf{z}},\overline{\textbf{x}})$ be a solution candidate obtained from solving the master problem. The algorithm then checks if the nonlinear constraints (ref) are feasible within an acceptance threshold (denoted as $\epsilon$), i.e., if $\overline{\theta}_n \leq \Psi_n(\overline{z}_n) +\epsilon$, for a given threshold $\epsilon>0$. If this condition holds true, the algorithm terminates and returns $(\overline{\pmb{\theta}},\overline{\textbf{z}},\overline{\textbf{x}})$. Otherwise, outer-approximation cuts based on $(\overline{\pmb{\theta}},\overline{\textbf{z}},\overline{\textbf{x}})$ are generated and added to the master problem to proceed to the next iteration. It can be shown that the above procedure is guaranteed to terminate after a finite number of iterations and return an optimal solution to (ref) duran1986outer.
In the approach described above, outer-approximation is performed as an iterative cutting-plane process, where outer-approximation (or sub-gradient) cuts are iteratively added to a master problem, which takes the form of a MILP. In the context of the competitive facility location problem, an outer-approximation approach can also be implemented differently using a piecewise linear approximation method ABOOLIAN2007a, ABOOLIAN2007b. In this approach, the univariate and concave functions $\Phi_n(\overline{z}_n)$ are approximated by piecewise linear concave functions. Specifically, the approximation of $\Phi_n(\overline{z}_n)$ is achieved by constructing sub-gradient cuts based on a set of breakpoints within the range of $\overline{z}_n$. This method is also referred to as the Tangent Line Approximation (TLA). Compared to the cutting-plane method mentioned earlier, this approach offers several advantages. Notably, it allows the nonlinear program in (ref) to be reformulated as a single MILP (with arbitrary precision) without introducing additional binary variables. This MILP can then be solved in one step to obtain a near-optimal solution, whereas the cutting-plane method requires solving a sequence of MILPs.
While outer-approximation in the form of cutting-plane methods has demonstrated state-of-the-art results in the context of competitive facility location without market expansion mai2020multicut, Ljubic2018outer, it is no longer an exact method when market expansion considerations are introduced, particularly when the market expansion function is non-concave. Therefore, in the following, we investigate outer (and inner) approximation approaches in the form of piecewise linear approximations. As mentioned earlier, this approach offers the advantage of approximating (ref) as a single MILP. This not only simplifies the problem structure but also provides a practical and efficient way to handle non-concave market expansion functions.
In the aforementioned outer-approximation (OA) approach, the mixed-integer nonlinear problem is tackled by approximating each concave component $\Psi_n(z_n)$ with concave piece-wise linear functions in $z_n$, enabling the solution of the ME-MCP through a sequence of MILPs. While achieving state-of-the-art performance in the context of the MCP, this outer-approximation approach is incapable of handling non-concave objective functions, becoming heuristic when the objective function is no-longer concave. In this section, we explore an alternative approach, called piece-wise linear inner-approximation (PWIA), which facilitates solving the ME-MCP by constructing piece-wise linear functions that inner-approximate $\Psi_n(z_n)$ internally. Our PWIA approach offers two advantages. First, as demonstrated later, such an inner-approximation function always yields smaller approximation errors compared to its outer-approximation counterpart. Second, as elucidated in the following section, under a general non-concave market expansion function, PWIA allows us to approximate the ME-MCP (with arbitrary precision) via MILPs, rendering it convenient for near-optimal solutions.
To facilitate our later exposition, let us first introduce formal definitions of piece-wise linear inner and outer approximations as below:
Now, given a concave piece-wise linear approximation function $\Gamma(t) = \min_{k\in[K]}\{a_k t + b_k\}$, let $\{(t_1,\Gamma(t_1));\ldots;(t_H,\Gamma(t_H))\}$ be $H$ “breakpoints” of $\Gamma(t)$, i.e. points where the function transitions from one linear segment to another within its piece-wise structure, such that $L = t_1<t_2<\ldots<t_H=U$ . Such breakpoints can be founds by considering all the intersection points of all pairs of the linear functions $\{a_k t + b_k,~k\in [K]\}$ and select points $(t^*,\Gamma(t^*))$ such that $\Gamma(t^*) \leq a_kt^*+b_k$ for all $k\in [K]$. The piece-wise linear function $\Gamma(t)$ can be equivalently represented as: \[ \Gamma(t) = \min_{h\in [H-1]} \left\{\Gamma(t_h) + \frac{\Gamma(t_{h+1}) - \Gamma(t_{h})}{t_{h+1} - t_h} (t-t_h) \right\}. \] It can be seen that $H$ is the minimum linear segments necessary to represent $\Gamma(t)$ in $[L,U]$. We are now ready to state our result saying that, given any piece-wise linear function $\Gamma(t)$ that outer-approximate a concave function, there are always another piece-wise linear approximation function that inner-approximates that concave function with the same number of necessary line segments, but yields smaller approximation errors. approximation gaps. We state this result in the following theorem.
The proof can be found in the appendix, which highlights that the inequality in (ref) is active (i.e., equality holds) only when the concave function $\Phi(t)$ exhibits uniform curvature across the interval \([L, U]\). This condition occurs if $\Phi(t)$ is either a linear function or takes the shape of a circle. The theorem and its proof further imply that for any piecewise linear outer-approximation of $\Phi(t)$, it is always possible to construct breakpoints within \([L, U]\) that yield a piecewise linear inner-approximation with a smaller approximation gap and the same number of line segments.
Later, we will demonstrate that such a piecewise linear approximation enables reformulation of the original problem as a MILP, with its size generally proportional to the numbers of line segments. Thus, the use of an inner-approximation proves to be more advantageous compared to its outer-approximation counterpart, particularly in terms of computational efficiency.
We begin by presenting a MILP approximation of (ref), where the nonlinear components are approximated using inner-approximation techniques. Following this, we discuss an approach to optimally select the linear segments for the inner-approximation, aiming to minimize the size of the resulting MILP formulation while ensuring a certain level of approximability.
Now we show how to approximate (ref) as a MILP using inner-approximation functions. We first let $L_n$ and $U_n$ be an upper bound and lower bound of $z_n$ in its feasible set. Such bounds can be estimated quickly by sorting $V_{ni}$, $i\in [m]$, in ascending order and select $C$ first elements for the lower bound, and $C$ last elements for the upper bound. This is possible because, if $\sigma_1,\ldots,\sigma_m$ is a permutation of $(1,\ldots,m)$ such that $V_{n\sigma_1}\leq \ldots \leq V_{n\sigma_m}$, then the following always holds true: \[ 1+\sum_{i =1}^{C}( V_{n\sigma_i}) \leq 1+ \sum_{i\in [m]} x_iV_{ni} \leq 1 + \sum_{i =m-C+1}^{m}( V_{n\sigma_i}) \] for all $\textbf{x} \in\{0,1\}^m$ such that $\sum_{i}x_i = C$. We can then select the lower and upper bounds for $z_n$ as $L_n = 1+\sum_{i =1}^{C}( V_{n\sigma_i}) $ and $U_n = 1+ \sum_{i =m-C+1}^{m}( V_{n\sigma_i})$.
To construct piece-wise linear functions that inner-approximate each component $\Psi_n(z_n)$ of the objective function, we split $[L_n; U_n]$ into $K_n$ sub-intervals $[c^n_k;c^n_{k+1}]$ for $k \in [K_n]$, where $c^n_k,~ k\in [K_n+1]$, are breakpoints such that $L_n = c^n_1 < c^n_2<\ldots< c^n_{K_n+1} = U_n$. We define the following piece-wise concave linear function: \[ \Gamma_n(z) = \min_{k\in [K_n-1]} \left\{\Psi_n(c^n_k) + \frac{\Psi_n(c^n_{k+1}) - \Psi_n(c^n_{k})}{c^n_{k+1} - c^n_k} (z-c^n_k) \right\},~ \forall n\in[N]. \] We then approximate each concave function $\Psi_n(z_n)$ by $\Gamma_n(z_n) $, resulting in the following mixed-integer nonlinear problem:
We then can see that (ref) can be reformulated as a MILP with no additional binary variables (Proposition (ref) below).
The proposition is obviously verified. The next theorem provides a performance guarantee for a solution returned by (ref).
Theorem (ref) tells us that we can obtain an $(N\epsilon)$-approximation solution if we select piece-wise linear functions such that $\max_{n\in [N]}\max_{z \in [L_n,U_n] } \left|\Psi_n(z) -\Gamma_n(z) \right| \leq \epsilon$. It is clear that this can be always achievable for any $\epsilon>0$ by selecting sufficiently small intervals, because $$\lim_{\max_{n\in [N],k\in [K_n+1]} |c^n_{k+1}- c^n_{k}|\rightarrow 0}\max_{z \in [L_n;U_n] } |\Phi_n(z) -\Gamma_n(z)| = 0.$$ However, increasing the number of breakpoints also results in the growth of the size of the approximate MILP (ref). Since we aim to optimize the size of (ref), in the following, we demonstrate how to select the breakpoints in a manner that minimizes $K_n$ while ensuring an approximation guarantee.
In this section, we explore an approach for minimizing the number of breakpoints while ensuring that the piece-wise linear approximation functions $\Gamma_n(z_n)$ remain within an $\epsilon$-neighborhood of the true objective functions $\Psi_n(z_n)$. To minimize the number of breakpoints, we would need to expand the gap between any consecutive breakpoints as much as possible, while guaranteeing that the approximation errors do not exceed a given threshold. That is, from any breakpoint $a\in [L_n, U_{n}]$ and given $\epsilon>0$, we need to find a next breakpoint $b>a$ such that \[ \max_{z\in [a,b]}\left|\Psi_n(z) - \Gamma_n(z)\right| \leq \epsilon, \] recalling that \[ \Gamma_n(z) = \left(\Psi_n(a) +\frac{\Psi_n(b) -\Psi_n(a) }{b-a} (z - a)\right), ~\forall z\in [a,b]. \] Since we want to minimize the number of line segments, we will need to choose $b$ in such a way that the gap $|b-a|$ is maximized. We then introduce the following problem to this end:
Let us define, for ease of notation, denote:
For solving (ref), we first introduce the following lemma showing some important properties of the above functions:
We now discuss how to solve (ref) using the monotonicity and convexity of $\Theta_n(t)$ and $\Lambda_n(t|a)$. We first write this problem as: \[ \max \left\{t\in [a, U_n]\Bigg|~ \Lambda_n(t|a)\leq \epsilon\right\}. \] Since $\Lambda_n(t|a)$ is (strictly) increasing in $t$ and $\Lambda_n(a|a) = 0$, the above problem always yields a unique optimal solution that can be found by a binary search procedure. Briefly, such a binary search can start with the interval $[l,u]$ where $l=a$ and $u = U_n$. We then check if $\Lambda_n(u)\leq \epsilon$ then return $t^* = u$ as an optimal solution. Otherwise we take middle point $r = (u+l)/2$ and compute $\Lambda_n(r|a)$. If $\Lambda_n(r|a)<\epsilon$ we update the interval as $[r,u]$, otherwise we update the next interval as $[l,u]$. This process stops when $u-l \leq \delta$ for a given threshold $\delta$. It is known that this procedure will terminate after ${\mathcal{O}}(\log(1/\delta))$ iterations.
Now, having an efficient method to solve (ref), we describe below our method to (optimally) calculate breakpoints for the inner-approximation:
We characterize the properties of the breakpoints returned by the above procedure in Theorem (ref) below (the proof is given in appendix):
Theorem (ref) establishes that the proposed procedure generates an optimal number of breakpoints. Specifically, there exists no other piecewise linear inner-approximation function with fewer breakpoints that achieves the same or smaller approximation gap compared to the one generated by the procedure. This result is intuitive, as the procedure optimizes each new breakpoint at every step. Consequently, for any smaller set of breakpoints, there will always be at least one pair of consecutive points where the approximation gap exceeds $\epsilon$.
The second part of Theorem (ref) highlights two important (and non-trivial) aspects. First, the breakpoint-finding procedure always terminates after a finite number of steps. Second, the number of steps (or generated breakpoints) is in $\mathcal{O}(1/\sqrt{\epsilon})$ and is generally proportional to the marginal value of the second-order derivative of $\Psi(z_n)$. This implies that the number of breakpoints increases to infinity as $\epsilon$ approaches zero. Moreover, the number of breakpoints will be larger if the concave function $\Psi(z_n)$ has high curvature and smaller if $\Psi(z_n)$ has low curvature (i.e., closer to a linear function). In the special case where $\Psi(z_n)$ is linear, the upper and lower bounds satisfy $L^{\Psi}_n = U^{\Psi}_n = 0$, and only one breakpoint is needed ($K_n = 0$), which aligns with expectations.
Our analysis thus far heavily relies on the assumption of concavity for the market expansion function. While such an assumption has been widely utilized in the literature and enables us to derive neat results (such as the concavity and submodularity of the objective function), aiding in efficiently solving the nonlinear optimization problem, it also presents certain limitations that may inaccurately capture market growth dynamics.
Specifically, the concavity assumption implies that the total demand of customer type $n$, calculated as $q_n g(u)$ (where $u$ represents the total expected customer utility offered by available facilities), grows rapidly when $u$ is small and gradually converges to $q_n$ as $u$ approaches infinity. However, this behavior may not be realistic as the addition of a few new facilities to the market would not immediately impact market growth. Conversely, it would be more realistic to assume that total demand grows slowly when $u$ is small and accelerates when a significant number of additional facilities are introduced to the market (resulting in a notable increase in $u$). To further illustrate this remark, Figure (ref) below depicts the market growth behavior under two popular concave functions $g_1(t) = \frac{t}{t+\alpha}$ and $g_2(t) = 1-e^{-\alpha t}$ (as mentioned previously) and a non-concave function (i.e., sigmoidal function $g_3(t) =\frac{1}{1+e^{-\alpha t}}$). We can observe that both \( g_1(t) \) and \( g_2(t) \) grow rapidly as \( t \) increases from 0, slowing down only when \( t \) becomes sufficiently large. Mathematically, this is because both \( g_1(t) \) and \( g_2(t) \) are concave, resulting in decreasing gradients with respect to \( t \). In contrast, \( g_3(t) \) exhibits smaller growth rates when \( t \) is small and increases faster as \( t \) becomes larger. Consequently, \( g_3(t) \) would better reflect the influence of customer utility on the market size in practical scenarios.
To address the aforementioned limitation of the concavity assumption, in this section, we will consider the ME-MCP with a general non-concave market-expansion function. We will first present a general method to approximate the ME-MCP via a MILP with arbitrary precision. We then demonstrate that, by identifying intervals where the objective function is either concave or convex, we can utilize the methods outlined earlier to optimally compute the breakpoints, thereby reducing the size of the MILP approximation. Furthermore, we will show that certain binary variables can be relaxed, further enhancing the efficiency of the approximate MILP.
We now consider the case that the market-expansion function $g(t)$ is not concave. As a result, the objective function $\Phi_n(z_n)$ is no-longer concave in $z_n$. We propose to approximate the non-concave function $\Phi_n(z_n)$ by a piece-wise linear function and show that (ref) can be approximated by an MILP with an arbitrary precision.
First, let us assume that $g(t)$ is twice-differentiable in $z$, implying that $\Psi_n(z)$ is also twice-differentiable in $z$ for all $n\in [N]$. By taking the second derivative of this function and find solutions to $\Psi''_n(z) = 0$, one can identify intervals in which $\Psi_n(z)$ is either convex or concave. This allows us to well optimize the line segments and reduce the number of additional binary variables. That is, assume that we can split $[L_n;U_n]$ into some sub-intervals such that $\Psi_n(\cdot)$ is either concave or convex in each sub-interval. For each sub-interval, if $\Psi_n(z_n)$ is concave,we can use the method above to further split it into smaller interval $[c^n_k; c^n_{k+1}]$ in such a way that the gap between $\Psi_n(z)$ and the piece-wise linear function $\Gamma_n(z)$ is less than $\epsilon$ for any $z\in [c^n_k; c^n_{k+1}]$. On the other hand, if $\Psi_n(z_n)$ is concave, we show in Appendix (ref) that one can use methods similar to those described in Subsection (ref) to optimize the number of intervals $[c^n_k; c^n_{k+1}]$. We will describe this in detail later in the section. However, before that, let us show how to approximate the ME-MCP with a non-concave market-expansion function by using a MILP with such breakpoints.
Let us assume that after this procedure we also obtain a sequence of sub-intervals $\{c^n_1,\ldots,c^n_{K_n+1}\}$ such that within each interval $[c^n_k,c^n_{k+1}]$, $k\in [K_n]$, the gap between $\Psi_n(z_n)$ and the linear function $\Gamma_n(z)$, defined as: \[ \Gamma_n(z) = \Psi_n(c^n_k) + \frac{\Psi_n(c^n_{k+1}) - \Psi_n(c^n_{k})}{c^n_{k+1} - c^n_k} (z-c^n_k), \] is not larger than an $\epsilon >0$. We now can approximate $\Psi_n(z)$ via the following piece-wise linear function:
where \[ \gamma^n_k = \frac{\Psi_n(c^n_{k+1}) - \Psi_n(c^n_{k}) }{c^n_{k+1}-c^n_{k}},~\forall n\in [N],k\in [K_n]. \] We now represent the condition \( z \in [c^n_k, c^n_{k+1}] \) using a binary variable \( y_{nk} \) and a continuous variable \( r_{nk} \). The binary variable \( y_{nk} \) satisfies \( y_{nk} \geq y_{n,{k+1}} \) for all \( k \in [K_n-1] \), and the continuous variable \( r_{nk} \) lies in the interval \([0,1]\) for all \( n \in [N] \) and \( k \in [K_n] \). Additionally, we require \( r_{nk} \geq y_{nk} \) and \( r_{n,k+1} \leq y_{nk} \) for all \( n \in [N] \) and \( k \in [K_n-1] \). This setup is to ensure that if \( y_{nk}=1 \), then \( r_{nk}=1 \); otherwise, if \( y_{nk}=0 \), then \( r_{nk'} = 0 \) for \( k' = k+1, \ldots \). The binary variables \( y_{nk} \) indicate the interval \([c^n_k, c^n_{k+1}]\) where \( z_n \) belongs, and the continuous variable \( r_{nk} \) captures the part \( z_n - c^n_k \). Using these variables, any \( z \in [L_n, U_n] \) can be expressed as \( z_n = \sum_{k \in [K_n-1]} (c^n_{k+1} - c^n_k) r_{nk} \). Moreover, the approximate function \( \Gamma_n(z) \) can be written as: \[ \Gamma_n(z) = \Psi_n(L_n) + \sum_{k \in [K_n-1]} \delta^n_k (c^n_{k+1} - c^n_k) r_{nk} \] We then can approximate the ME-MCP by the following piece-wise linear problem:
This case differs from the concave market expansion scenario in that additional binary variables are required to construct the MILP approximation of the facility location problem. This raises concerns when a highly accurate approximation is needed, as the number of additional binary variables is proportional to the number of breakpoints used to form the piecewise linear function. In the subsequent section, we will demonstrate that some of these additional variables can be relaxed, resulting in a significantly simplified MILP approximation formulation.
Before discussing this relaxation, we state the following theorem showing that solving (ref) provides a solution with the same performance guarantees as solving (ref) in the concave market expansion case considered earlier.
As mentioned earlier, in the case of general market expansion functions, we can minimize the number of breakpoints by dividing the range $[L_n, U_n]$, for any $n \in [N]$, into sub-intervals where the objective function $\Psi_n(z_n)$ is either concave or convex in $z_n$. We can then apply the techniques described in Section (ref) (for concave intervals) and in Appendix (ref) (for convex ones) \footnote{This application is generally straightforward, as a convex function can be viewed as the inverse of a concave function.} The detailed steps are described as follows:
From Theorem (ref) and Appendix (ref), we can see that within any interval $[a, b]$ where $\Psi_n(z)$ is either concave or convex, the above procedure guarantees that the number of breakpoints is minimized. Moreover, Proposition (ref) below states that this procedure always terminates after a finite number of iterations and provides an upper bound on the number of breakpoints generated.
The proof can be found in the appendix where we leverage the second-order Taylor expansion to establish the bound. Similar to the case of concave market expansion, the number of breakpoints generated by the [Finding Breakpoints] procedure is always finite and bounded above by $\mathcal{O}(U^\Psi_n / \sqrt{\epsilon})$. Consequently, a higher number of breakpoints (and thus a larger MILP size) will be required if the desired accuracy $\epsilon$ is small or if the curvature of $\Psi_n(z)$ within \([L_n, U_n]\) is high. Conversely, fewer breakpoints will be needed if the curvature is low or the approximation accuracy requirement is less stringent.
As described in the previous section, the breakpoints are generated by dividing the interval $[L_n, U_n]$ (for any $n \in [N]$) into sub-intervals where $\Psi_n(z)$ is either concave or convex. The main problem (ref) can then be approximated by (ref), whose size is proportional to the number of breakpoints. Specifically, (ref) requires $\sum_{n} K_n$ additional binary variables. According to Proposition (ref), the number of additional binary variables is proportional to $\frac{1}{\sqrt{\epsilon}}$, which increases as $\epsilon$ approaches zero. In the following, we show that the number of breakpoints (and thus the number of additional binary variables) can be significantly reduced by relaxing part of the additional binary variables.
Let define ${\mathcal{K}}_n\subset [K_n]$ such that $\Psi_n(z)$ is concave in $[c^n_{k}; c^n_{k+1}]$ for all $k\in {\mathcal{K}}_n$. We have the following theorem stating that all the binary variables $y_{nk}$ for all $k\in {\mathcal{K}}_n$ can be safely relaxed.
The proof can be found in the appendix. Theorem (ref) indicates that some of the additional variables associated with regions where the function $\Psi_n(z_n)$ is concave can be relaxed. Specifically, if $\Psi_n(z_n)$ is concave across the entire interval \([L_n, U_n]\), all the additional binary variables can be relaxed, as in the case of the concave market expansion scenario discussed earlier.
Since it is expected that the market expansion function \(g(t)\) will always increase and converge to 1 as \(z\) approaches infinity, \(\Psi_n(z_n)\) remains concave when \(z_n\) is sufficiently large. This allows a significant portion of the additional variables to be safely relaxed, thereby improving the overall computational efficiency of solving the MILP approximation.
This section presents the experimental results to assess the performance of three solutions methods introduced in Section (ref). In particular, the first Subsection (ref) describes the benchmark datasets and experimental settings. The next Subsection (ref) present a sensitivity analysis for choosing the error threshold $\epsilon$ in the Piece-wise Inner-approximation method. Subsection (ref) provides the computational results under the concave market expansion setting. Finally, Subsection (ref) presents the results on the general non-concave market expansion functions.
We utilize three benchmark datasets in our experiments, all of which are widely used in prior work in the context of competitive facility location Ljubic2018outer, mai2020multicut.
For each test instance, the number of open locations $H$ is varied over $\{2, 3, \ldots, 10\}$. The utility associated with the customer zone $n$ and location $i$ is given by $v_{ni} = - \theta c_{ni}$ for $i\in[m]$ and $v^c_{ni} = - \gamma\theta c_{ni}$ for $i\in S^c$, where $S^c$ is randomly sampled from $[m]$ with $|S^c| = \ceil{m / 10}$, $\theta \in \{1, 5, 10\}$ and $\gamma \in \{0.01, 0.1, 1\}$ for the HM14 and OR Lib datasets, and $\theta \in \{0.5, 1, 2\}$ and $\gamma \in \{0.5, 1, 2\}$ for the NYC. Combining all the parameters results in total of 1215 test instances of the HM14 dataset, 891 instances of the \texttt{ORLib} dataset, and 81 instances of the \texttt{NYC} dataset. For the objective function, we set the trade-off parameter $\lambda$ to 1.
For comparison, since there are no direct solution methods capable of solving the problem under consideration, we adapt state-of-the-art methods developed in the existing literature. Specifically, we include the following three approaches for comparison:
The experiments are implemented by C++ and run on Intel(R) Xeon(R) CPU E5-2698 v3 @ 2.30GHz. All linear programs are carried out by IBM ILOG CPLEX 22.1, with the time limit for each linear program being set to 5 hours. Each method under consideration (PIA, OA and LS) is given a time budget of 1 hours.
We begin by conducting an experiment to analyze the practical impact of the parameter \(\epsilon\) on the performance of the PIA method. For this purpose, we select a concave market expansion function \(g(t) = 1 - e^{-t}\) and run the PIA method on three instances of the HM14 dataset, where the number of customers is fixed at \(N = 50\) and \(m \in \{25, 50, 100\}\). The value of \(\epsilon\) is varied from \(1\text{E-5}\) to \(1.0\). For each value of \(\epsilon\), we measure and report the runtime and the percentage error of the corresponding solution relative to the solution obtained with \(\epsilon = 1\text{E-5}\). Here, we assume that setting \(\epsilon\) to \(1\text{E-5}\) will generally yield optimal solutions. The percentage errors and runtimes are plotted in Figure (ref).
For smaller values of \(\epsilon\) (e.g., \(1\text{E-5}\), \(1\text{E-4}\), and \(1\text{E-3}\)), the Error (%) remains consistently zero for all problem sizes \(m\). This indicates that PIA achieves almost-optimal solutions when \(\epsilon\) is set to a sufficiently small value. However, this accuracy comes at the cost of increased computational time. For instance, when \(m = 50\), the runtime is 6.7 seconds for \(\epsilon = 1\text{E-5}\) and reduces to 2.4 seconds for \(\epsilon = 1\text{E-4}\), showing that even small increases in \(\epsilon\) can lead to significant improvements in efficiency.
As \(\epsilon\) increases to \(1\text{E-2}\), \(1\text{E-1}\), and \(1\text{E+0}\), a slight error begins to emerge, particularly for larger problem sizes. For example, at \(m = 50\), the error increases to 0.031% when \(\epsilon = 1\text{E-1}\). Similarly, at \(m = 100\), the error increases to 0.066% for both \(\epsilon = 1\text{E-1}\) and \(\epsilon = 1\text{E+0}\). Despite the presence of these small errors, the runtime decreases significantly. For \(m = 25\), the runtime reduces from 4.8 seconds (\(\epsilon = 1\text{E-5}\)) to just 0.6 seconds (\(\epsilon = 1\text{E+0}\)). This demonstrates that larger values of \(\epsilon\) lead to coarser approximations that accelerate computation but slightly compromise accuracy.
The results also reveal the scalability of the PIA method with respect to the problem size \(m\). As \(m\) increases from \(25\) to \(100\), the runtime increases, particularly for smaller values of \(\epsilon\). For instance, at \(\epsilon = 1\text{E-5}\), the runtime grows from 4.8 seconds for \(m = 25\) to 6.69 seconds for \(m = 50\). However, for larger values of \(\epsilon\), such as \(1\text{E-1}\) or \(1\text{E+0}\), the runtime remains relatively low even as \(m\) increases. This suggests that the computational burden of PIA can be effectively mitigated by selecting a larger \(\epsilon\) when slight errors are acceptable.
Overall, the results demonstrate that the choice of \(\epsilon\) is critical in balancing solution accuracy and computational efficiency. Smaller values of \(\epsilon\) are suitable for applications requiring high accuracy, as they ensure optimal solutions at the cost of longer runtimes. On the other hand, larger values of \(\epsilon\) significantly reduce runtime while maintaining near-optimal solutions, making them ideal for scenarios where computational speed is prioritized. Based on these analyses, we select \(\epsilon = 1\text{E-2}\) for our comparison results, as it appears to ensure an almost optimal solution for the PIA method while maintaining a reasonable size for the approximation problem.
In this section, we present the numerical results obtained by three solution methods for addressing the facility location problem with concave market expansion. The market expansion function is selected as \( g(t) = 1 - e^{-\alpha t} \) with \(\alpha = 1\), a popular choice in prior studies related to market expansion in competitive facility location problems ABOOLIAN2007a,lin2022locating. The results for three datasets are reported in Table (ref), where each row contains results for instances grouped by \((N, m)\).
Three evaluation criteria are considered: (1) the number of instances solved to optimality within the time budget, (2) the number of instances where the corresponding method achieves the best solution among the three methods, and (3) the average computing time in seconds required to confirm the optimality of the solution. In this setting, since the objective function is concave (Theorem (ref)), both PIA and OA serve as exact (or near-exact) methods, while LS remains heuristic. Therefore, the number of solved instances is only reported for PIA and OA.
The results generally show that PIA emerges as the most efficient method, consistently solving all instances across all datasets and configurations. This reliability is evident in both the HM14 and cap datasets, where PIA solves all 81 and 324 instances, respectively, achieving the best objective value in every case. Furthermore, its computational efficiency is particularly noteworthy, especially for larger datasets such as cap_abc, where PIA completes the task in 13.35 seconds. This combination of reliability, optimality, and efficiency positions PIA as the most favorable method for solving these optimization problems.
OA closely mirrors the performance of PIA in terms of solution quality and reliability. It also solves all instances across the datasets and achieves the best objective value in every case. However, OA tends to require slightly higher computational times, particularly for larger datasets. For example, in the NYC dataset, OA’s computing time (5547.00 seconds) is significantly higher than that of PIA (1433.80 seconds). Despite this, OA remains a viable choice for scenarios where computational cost is less of a concern, given its ability to consistently deliver high-quality solutions. This observation aligns with the fact that OA has been recognized as a state-of-the-art approach for competitive facility location problems under fixed market sizes mai2020multicut.
LS, on the other hand, provides a contrasting performance profile. While LS often requires less computational time compared to PIA and OA, as demonstrated in the NYC dataset (973.86 seconds), it does not guarantee optimal or near-optimal solutions. This limitation is reflected in the “-” entries under the “Number of solved instances” column for LS. These entries highlight that LS, being a heuristic method, prioritizes computational speed over solution quality. Although LS can occasionally match the best objective values achieved by PIA and OA, such occurrences are less consistent. As a result, LS is less suitable for applications where solution quality or optimality is critical.
The scalability of PIA and OA across increasing problem sizes further underscores their suitability for large-scale instances. As the values of \(N\) and \(m\) grow, both methods maintain their ability to solve all instances while achieving the best objective values. In contrast, LS, despite its computational efficiency, struggles to balance scalability and solution quality, particularly in larger datasets.
In summary, the results highlight that PIA stands out as the most reliable and efficient method, particularly for scenarios requiring optimal solutions. OA offers a strong alternative, especially for smaller datasets, though it may incur higher computational costs for larger problems. LS, with its emphasis on computational speed, is best suited for applications where solution quality is less critical, and computational resources are limited.
In this experiment, we evaluate the performance of PIA under a general non-concave market expansion function. The market expansion function is defined as \( g(t) = \frac{1}{1 + e^{-\alpha(t - \beta)}} \), where \(\alpha = 5\) and \(\beta = 4\). The results are presented in Table (ref), using the same format as in the previous experiment. Since both OA and LS are heuristic methods in this setting, we report the number of solved instances only for PIA.
Similar to the previous experiment, PIA emerges as the most efficient method. It consistently solves all instances across the different datasets and configurations, as reflected in the “No. of solved instances” column. Unlike OA and LS, PIA guarantees optimal or near-optimal solutions. This highlights its ability to handle the complexity of the solution space, particularly in cases where other methods fail. For example, in the HM14 and cap datasets, PIA solves all instances while achieving the best objective values. This makes PIA the preferred choice for problems requiring both practical accuracy and reliability.
The analysis of computing times provides further insights into the trade-offs between solution quality and efficiency. While ensuring optimal or near-optimal solutions, PIA maintains competitive computing times across all problem sizes. For instance, in the cap_abc dataset with \(N = 1000, m = 100\), PIA completes the task in 36.43 seconds, which is slower than OA (1.56 seconds) but significantly faster than LS (57.17 seconds). OA often demonstrates shorter computational times, particularly for smaller datasets, but this efficiency comes at the cost of reduced robustness. Notably, the unusually fast runtime of OA in the second dataset coincides with its poor solution quality, which can be attributed to invalid cutting planes introduced at the early stages of the algorithm. LS, on the other hand, achieves the fastest computing times in some cases, such as the NYC dataset, but its inability to guarantee solution quality undermines its overall performance.
Scalability is another critical factor. PIA demonstrates strong scalability as the problem size increases, maintaining its ability to solve all instances even for large datasets. For example, in the NYC dataset (\(N = 82341\)), PIA successfully solves all instances, achieving the best objective values while maintaining a reasonable computational cost. In contrast, OA and LS struggle to scale effectively, with performance deteriorating as the problem size increases. This issue is particularly pronounced in the larger datasets, such as cap_abc and NYC, where neither OA nor LS matches the robustness and reliability observed in PIA.
The table also highlights interesting results regarding solution quality. For the HM14 and NYC datasets, all methods achieve comparable solution quality; however, PIA stands out as the fastest method in these instances. In the second dataset, ORlib, PIA demonstrates superior performance by providing the best solutions for all 324 test instances in the cap_10 and cap_13 cases. In contrast, OA solves only 268 and 288 instances, while LS solves 308 and 276 instances, respectively. For the large cap_abc dataset within ORlib, PIA solves 222 out of 324 instances to optimality. Despite this limitation, the number of best solutions found by PIA remains comparable to those obtained by OA and LS, further reinforcing its overall reliability and efficiency.
In summary, the results clearly demonstrate that PIA is the most effective method for solving the facility location problem under general non-concave market expansion functions. PIA guarantees near-optimal solutions while maintaining competitive computational efficiency and strong scalability. While OA and LS offer faster runtimes in specific cases, their inability to consistently solve instances and ensure solution quality limits their applicability. For problems requiring reliability, accuracy, and scalability, PIA remains the method of choice.
The slope of the market expansion function reflects how the market grows as the total consumer surplus increases. In the context of concave market expansion with \( g(t) = 1 - e^{-\alpha t} \), this behavior is captured by the parameter \(\alpha\). Since \(g(t)\) is an increasing function of \(\alpha\), and the second-order derivative of \(g(t)\) with respect to \(t\) is given by \( g''(t) = \alpha^2 e^{-\alpha t} \), it decreases exponentially to zero as \(\alpha\) increases. Intuitively, when \(\alpha\) is large, the market expansion function increases more rapidly towards 1 and exhibits lower curvature. On the other hand, when \(\alpha\) is smaller, the market expansion function adheres to a higher curvature. Moreover, the bounds reported in Theorem (ref) indicate that functions with lower curvature require fewer breakpoints in the PIA, and vice versa. Consequently, a higher \(\alpha\) leads to a lower curvature of the objective function, resulting in a smaller approximation problem. As a result, the PIA method is expected to run faster when \(\alpha\) is larger. To experimentally illustrate this, we conduct a series of experiments with varying values of \(\alpha\).
To this end, we choose the concave market expansion function \( g(t) = 1 - e^{-\alpha t} \) and vary the parameter \(\alpha \in \{1, 5, 10, 20, 50, 100\}\). We select the ORlib dataset (excluding cap_abc due to its large size) for this experiment since it is most sensitive to the concavity of \(g(t)\). The results are plotted in Figure (ref). As expected, when \(\alpha\) is small, PIA requires more computational time, whereas it becomes faster as \(\alpha\) increases. This aligns well with the intuition discussed earlier: larger values of \(\alpha\) reduce the curvature of the objective function, thereby requiring fewer breakpoints for the PIA method. In contrast, the runtimes of OA and LS are not significantly affected by changes in \(\alpha\) -- their overall runtimes remain stable when \(\alpha\) increases.
In this paper, we studied a competitive facility location problem with market expansion and a customer-centric objective, aiming to capture the dynamics of the market while improving overall customer satisfaction. The novel problem formulation, to the best of our knowledge, cannot be directly solved to near-optimality by any existing approach, particularly under a general non-concave market expansion model.
To address these challenges, we first demonstrated that under concave market expansion, the objective function exhibits both concavity and submodularity. This allows the problem to be solved exactly using an outer-approximation approach. However, this property does not hold under a general non-concave market expansion function. To overcome this limitation, we proposed a new approach based on an inner-approximation method. We showed that our PIA approach consistently yields smaller approximation gaps compared to any outer-approximation counterpart. Furthermore, the inner-approximation program, in addition to being able to achieve arbitrarily precise solutions, can be formulated as a MILP.
We further strengthened the proposed approach by developing an optimal strategy for selecting breakpoints in the PIA, minimizing the size of the approximation problem for a given precision level. Additionally, we showed how to significantly reduce the number of binary variables in the case of non-concave market expansion by examining regions of the objective function where it behaves either as convex or concave.
Experiments conducted for both concave and non-concave market expansion settings demonstrate the efficiency of the proposed PIA approach in terms of solution quality, solution guarantees, and runtime performance. We also provided an analysis of the impact of the approximation accuracy threshold \(\epsilon\) and the slope parameter \(\alpha\) of the market expansion function on the performance of PIA.
Future research will focus on developing an advanced version of PIA that returns exact solutions or extending the proposed PIA approach to other variants of the competitive facility location problem. For instance, it could be applied to models involving more complex choice behaviors, such as the nested logit or multi-level nested logit models train2009discrete,mai2017dynamic.