An Approximation Algorithm for a Competitive Facility
Location Problem with Network Effects
Ling-Chieh Kung* and Wei-Hung Liao
Department of Information Management, National Taiwan University.
November 23, 2017
Abstract
When facilities are built to serve end consumers directly, it is natural that consumer demands are affected by the number of open facilities. Moreover, sometimes a facility becomes more attractive if other facilities around it are built. To capture these factors, in this study we construct a discrete location model for profit maximization with endogenous consumer demands and network effects. The effective demand is then a concave function of the sum of benefits of open facilities due to the diminishing marginal benefit effect. When the function is linear, we design a polynomial-time algorithm to find an optimal solution. When it is nonlinear, we show that the problem is NP-hard and develop an approximation algorithm based on demand function approximation, linear relaxation, decomposition, and sorting. It is demonstrated that the proposed algorithm has worst-case performance guarantees for some special cases of our problem. Numerical studies are conducted to demonstrate the average performance and general applicability of our algorithms.
Keywords: location, competitive facility location problem, network effect, approximation algorithm.
1
Introduction
The number and locations of facilities often affect consumers’ willingness of buying a product or using a service. In fact, in many cases consumer demands may be driven up only if
suf-*Corresponding author: [email protected] [email protected]
ficient facilities are built. Vehicle sharing systems provide a good illustration. As one of the most famous United States car sharing system nowadays, Zipcar is well known due to its new technologies and flexible renting plans. Nevertheless, one would sign up as a member of Zipcar only if she may easily find available cars nearby her location when she needs one. This can be achieved only if Zipcar sets up enough parking spaces in a city. A similar situation happens to public bicycle sharing system around the world. Almost all these systems can become a success only after enough rental sites have been built and enough bicycles have been supplied. As a non-transportation example, when a consumer wants to buy a personal computer or laptop, she would also evaluate how easy it is to find a warranty station when the product is broken. The distribution of warranty stations thus affects demands. Similar stories also exist for convenience stores, grocery stores, and charging or battery swapping stations for electronic vehicles.
Interestingly, opening a facility not only affects consumers directly but also changes how other facilities affect consumer demands. The first factor to consider is the interdependence among facilities. Consider the vehicle sharing business as an example. Many users of public bikes travel from subway stations to their offices, schools, and home. One would not find the service attractive even if there is one rental site just at the front door of her home. Instead, she will find it is much more attractive if there is another one at a subway station in the neighborhood. In general, while building a facility may attract consumers by itself, it may also create additional attractiveness together with other existing facilities. This is the upside of building a new facility. Nevertheless, there is also a downside: the marginal benefit of building facilities is diminishing. When one’s home/office has been surrounded by many public bike stations, building one more is not so attractive. All the aforementioned effects should be considered altogether when one makes facility construction decisions.
In this study, we investigate a profit-maximizing service provider’s facility location problem. Given a subset of locations, the service provider chooses a subset of locations to build facilities by considering two major types of effects: (1) the stand-alone benefit of a single facility and (2) the network benefit between a pair of facilities. To capture the diminishing marginal benefit property, we model a consumer’s willingness to use the service as a nondecreasing concave function of the sum of all benefits. Therefore, the sum of all benefits is converted into the consumer demand (and thus the service revenue) by a nondecreasing concave effective demand function. By considering the total service revenue and total cost of building facilities, the service provider decides where to build facilities to maximize her profit. The problem is thus
formulated as a nonseparable nonlinear integer program, whose objective function is to maximize a nonmonotone submodular function.
We start our analysis by considering the solvability of this problem and find that the shape of the effective demand function plays an important role. When the function is linear, the problem (which is still a nonlinear integer program with the presence of products of decision variables) can be solved in polynomial time. We solve this problem by reducing it to the maximum flow problem. When the function is nonlinear, however, the problem is at least weakly NP-hard even if there is no network benefit. This is shown by a reduction from the partition problem.
As one of the most common procedures to approach NP-hard problems is to develop heuris-tic algorithms, we propose one based on approximating the effective demand function, relaxing integer constraints, decomposing the problem into subproblems, and sorting facilities to make the construction decisions. The algorithm is thus named the approximation-relaxation-sorting-aggregation algorithm (ARSA). We prove that ARSA exhibits worst-case performance guaran-tees in a few special cases of our problem. This makes ARSA an approximation algorithm for the special cases (Williamson and Shmoys, 2011). To investigate its average-case performance, we compare ARSA with an exponential-time exact algorithm and a genetic algorithm in various scenarios through numerical experiments. We demonstrate that ARSA’s average-case perfor-mance is much better than its perforperfor-mance guarantees and the genetic algorithm. In many cases, it generates a solution that is close to an optimal solution.
Facility location problems with diminishing marginal benefit of building facilities are called competitive facility location problems (Karakitsiou, 2015). In these problems, facilities “com-pete” with each other in winning consumers. The competitive facility location problem has been studied from many perspectives. However, to the best of our knowledge, we are the first to develop an approximation algorithm while taking network effects into consideration. This is the major contribution of our study. Our study certainly has its limitations. In particular, we do not explicitly investigate how to quantify the form of the diminishing marginal benefit and estimate the stand-alone and network benefits. While this issue is ignored by almost all works in the literature of competitive facility location problems, we note that there is another stream of literature of developing data analytics methods to estimate consumer demands (some recent works include Hsieh et al. (2015), Tiwari and Kaushik (2014), among others). Interested readers may want to study these works. In Section 6, we briefly suggest how one may estimate parameters and apply our model and algorithm in practice. We hope one day we may connect
the two streams to make further contributions.
The remainder of this study is organized as follows. In Section 2, we discuss related works. In Section 3, we give a model formulation for our problem and conduct solvability analysis. The description of the ARSA algorithm and proofs of the worst-case performance guarantees are in Section 4. In Section 5, a numerical study is conducted to demonstrate the average-case performance and general applicability of our proposed algorithm. Concluding remarks and future works are provided in Section 6.
2
Literature review
Facility location problems have been widely studied in the past decades (Owen and Daskin, 1998; Daskin, 2013). Most traditional works are conducted by assuming that consumer demands are exogenous, i.e., not affected by the facility location decision. As the interrelationship between location decisions and demand endogeneity cannot be ignored in some situations, researchers are motivated to study competitive facility location problems (Karakitsiou, 2015).
In a static competitive facility location problem, a firm builds facilities to attract consumers (Aboolian et al., 2007a,b; Berman and Krass, 1998, 2002; Wu and Lin, 2003). Under the assumption that a consumer will only be served by one facility, facilities compete with each other. In these studies, the total market demand, as a function of the number of facilities, is generally assumed to be increasing for the market expansion effect and concave for the market cannibalization effect. In this study, we follow this stream and model the consumer demand by a concave function. Nevertheless, our model include the network effect between a pair of facilities, which is missing in most of the past works. In another stream of literature, researchers study sequential competitive facility location problems with a leader and a follower take turns to build facilities (Beresnev, 2014; K¨u¸c¨ukayadın et al., 2011, 2012; Mel’nikov, 2014). Each firm is assumed to be a profit maximizer, and building facilities is the key to win consumers from the competitor. In this study, we investigate a static facility location problem.
Most competitive facility location problems require one to maximize a submodular objec-tive function. Beside competiobjec-tive facility location problems, submodular function maximization appears in many kinds of problems, e.g., maximum coverage, expected utility maximization with discrete choices of a risk-averse decision maker, and combinatorial auctions with submod-ular utilities (Ahmed and Atamt¨urk, 2011). As submodular function maximization is NP-hard,
many studies appear to propose polynomial-time approximation algorithms, which are heuristic algorithms that possess worst-case performance guarantees (Williamson and Shmoys, 2011). For a general submodular function maximization problem, Nemhauser and Wolsey (1978) com-bine exhaustion and greedy search while Nemhauser et al. (1978) adopt greedy, local search, and linear programming relaxation to develop approximation algorithms. Unfortunately, their results only apply to nondecreasing submodular functions. As our problem has a nonmonotone objective function, it cannot be solved by their algorithms with a performance guarantee.
For nonmonotone submodular function maximization, a special case is the uncapacitated facility location problem. Several researchers have developed approximation algorithms for this problem, and the best performance guarantee known so far is 0.828 (Ageev and Sviridenko, 1999; Cornuejols et al., 1977a,b). For general nonmonotone submodular function maximization, Feige et al. (2011) develop several deterministic and randomized approximation algorithms. Established approximation factors include 25 for non-symmetric cases and 12 for symmetric cases. For our facility location problem, we show that our algorithm achieves an above 12 approximation factor under some explicit conditions.
To the best of our knowledge, the problem and method studied by Aboolian et al. (2007b) is the closest to our work. In their problem, a decision maker chooses several locations to build facilities to attract consumers at several markets. Each open facility gives each market a utility, which is negatively affected by the distance. The amount of realized demand at a market is a negative exponential function of the total utilities that the market obtains from all open facilities. The proportion of demand that goes to a facility is then the proportion of the total utility that comes from the facility. The decision maker’s problem is to maximize the total revenue subject to a fixed budget of opening facility. To solve the problem, they first show that the objective function is smoothly concave. They then approximate the smooth concave objective function by a piecewise linear function, which is larger than the original one by at most α%.1 This allows them to transform the original problem to an integer linear program, whose exact solution deviates from an optimal solution to the original problem by at most α%. Though this procedure can achieve any performance guarantee, it does not run in polynomial time. Therefore, they also developed a greedy algorithm having a performance guarantee 1 −1e by applying a result in Nemhauser et al. (1978). Because the objective function in their study is nondecreasing, both these two algorithms do not apply to our problem (whose objective function
1
is nonmonotone). Nevertheless, their strategy of approximating a given objective function by a piecewise linear one turns out to be a critical step in our proposed algorithm.
3
Model and solvability
3.1 The service provider’s problem
Suppose that there is a set of locations I = {1, 2, ..., n}, |I| = n, to build facilities. Let E be the set of all undirected edges [i, j] where i ∈ I, j ∈ I, and i < j. As |I| = n, we have |E| = n(n−1)2 . For facility location i, the decision variable xi = 1 means there is a facility in location i and
0 otherwise. Once a facility is built at location i, it will not only increase the demands but also affect the impact of other facilities. Let si ≥ 0 be a coefficient for the facility in location i
and tij ≥ 0 be a coefficient for the facilities between location i and j, we assume the effective
demand is g X i∈I sixi+ X [i,j]∈E tijxixj ! ,
where g(·) is a nondecreasing and concave function satisfying g(0) = 0. Note that the sum of benefits collected from open facilities is input into a concave function to generate the demand volume. This setting is also adopted by Aboolian et al. (2007b) (they adopt the term “utility” rather than “benefit”). To facilitate discussion, we refer to the function g(·) as the effective demand function, si as the stand-alone benefit, tij as the network benefit, and at some times use
“facility i” as an abbreviation of “the facility at location i.”
The stand-alone benefit si measures how a facility at location i may affect the demand by
itself. When si = 0, building a facility at location i does not drive up the sales itself; when
si > 0, however, the demand volume will be affected even if there exists no other facility. For
the case of building warranty stations for PCs and laptops, we should have si > 0 for all i.
On the contrary, for a public bike sharing system, most si would be close to 0.2 The network
benefit tij measures how building facilities i and j together may further bring up the demand.
For a bike sharing system, if locations i and j are not too far away, we should have tij > 0. In
this case, the distance between i and j can be one factor (but not the only factor) determining tij. For warranty stations, however, it should be reasonable to set tij = 0. Note that tij is
measuring how attractive it is to have both facilities i and j open; it has nothing to do with
2
the competition among facilities and splitting consumer demand to multiple facilities. The diminishing marginal benefit of building facilities is modeled by the effective demand function g(·), not tij. This is why in our problem tij ≥ 0 and should not be negative.
The service provider’s problem is to make the facility location decision so that the total profit is maximized. Her complete problem is
max xi∈{0,1} γg X i∈I sixi+ X [i,j]∈E tijxixj ! −X i∈I hixi
where γ > 0 is the price of the service and hi > 0 is the cost for building a facility at location
i. The objective function consists of two parts, the sales revenue and the construction cost. In this market, the effective demand is g(P
i∈Isixi+
P
[i,j]∈Etijxixj), the demand brought by
the stand-alone and network benefits of those open facilities. Multiplying the effective demand by the unit price γ results in the sales revenue. The service provider’s profit is then the sales revenue minus the construction costs P
i∈Ihixi. The only constraint is the binary constraint
for building facilities. Without loss of generality, we normalize γ to 1 to obtain
max xi∈{0,1} g X i∈I sixi+ X [i,j]∈E tijxixj ! −X i∈I hixi (1)
as our facility location problem in this study. For ease of exposition, we use z(x) to denote the objective value associated with a solution x = (x1, x2, ..., xn). Moreover, for location i we define
si− hi as its net stand-alone benefit. It turns out to play a critical role in the analysis of our
algorithms.
Table 1 lists all the notations mentioned above.
3.2 Effective demand function and solvability
The solvability of our facility location problem in (1) critically depends on the shape of the effective demand function g(·). For a special case that g(·) is negative exponential, Ahmed and Atamt¨urk (2011) prove in their Theorem 1 that the problem is (weakly) NP-hard even if tij = 0.
Below in Proposition 1 we extend their proof to any g(·) such that g(w) − kw has a unique maximizer in (0, ∞) for some k > 0, which is true for most concave g(·) that is nonlinear. Note that if g(·) is linear, g(w) − kw will not have a unique maximizer over the open interval (0, ∞). Therefore, whether our problem is NP-hard with a linear g(·) is not answered by Proposition 1.
Parameters I the set of candidate facility locations
E the set of all undirected edges connecting a pair of locations n The number of candidate facility locations
si the stand-alone benefit of facility i
tij the network benefit between facilities i and j
hi construction cost of facility i
g(B) the effective demand given the sum of all benefits B Decision variables
xi 1 if a facility is built at location i or 0 otherwise
Table 1: Notations
Proposition 1. Consider the problem defined in (1). Suppose that the function g satisfies the following condition: There exists some k > 0 such that the optimization problem
max
w≥0 g(w) − kw
has a unique positive optimal solution. In this case, the problem is NP-hard.
Proof. Let the unique optimal solution be w∗. Given a partition problem of integers a1, a2, ...,
and an that looks for a set S such thatPi∈Sai =Pi /∈Sai, we first normalize the integers so
thatP
i∈Sai = w∗. Then we construct an instance of our problem as
z∗= max xi∈{0,1} g X i∈I aixi ! −X i∈I kaixi,
i.e., we let si = ai, hi = kai, and tij = 0, where k > 0 is the value that makes g(w) − kw
having a unique maximizer in (0, ∞). We now claim that such a partition exists if and only if z∗ = g(w∗) − kw∗. If such a partition S exists, we may set xi = 1 for i ∈ S and 0 otherwise.
This results in the objective value g(P
i∈Sai) −Pi∈Skai = g(w∗) − kw∗. If our facility location
problem is solved to achieve z∗, as the maximizer w∗ is unique, it must be the case that we have selected a set S of facilities such thatP
i∈Sai = w∗. Such a set S then gives us a partition.
The above result states that our facility location problem is NP-hard as long as the problem exhibit diminishing marginal benefit. However, it leaves the case of constant marginal benefit
unanswered. Below we assume that g(w) = aw for some constant a > 0 and show a way to solve our problem in polynomial time.3
Consider the problem of maximizing
a X i∈I s0ixi+ X [i,j]∈E t0ijxixj ! −X i∈I h0ixi,
where a > 0 is a given constant. For this problem, we may without loss of generality normalize a to 1 by setting si = as0i, tij = at0ij, and hi = h0i. We then face an equivalent problem of
maximizing X i∈I sixi+ X [i,j]∈E tijxixj− X i∈I hixi= X i∈I (si− hi)xi+ X [i,j]∈E tijxixj
Suppose that for facility i1 we have si1− hi1 > 0, obviously xi1 = 1 in any optimal solution (as
tij ≥ 0). Therefore, we may without loss of generality assume that si− hi< 0 for all i ∈ I.
For our facility location problem in (1) where g(w) = w and si − hi < 0 for all i, we now
demonstrate how to reduce our problem to a maximum flow problem. To do so, we define a directed graph G = (V, A) as follows. The set of nodes V = {s} ∪ V1∪ V2∪ {t}, where s is the
source node, V1= E contains the set of location pairs, V2 = I contains the set of locations, and t
is the terminal node. The set of directed arcs A = A1∪A2∪A3, where A1 = {(s, [i, j]) | [i, j] ∈ V1}
is the set of links from s to V1, A2= {([i, j], i) | [i, j] ∈ V1, i ∈ V2}∪{([i, j], j) | [i, j] ∈ V1, j ∈ V2}
contains links from V1 to V2, and A3 = {(i, t) | i ∈ V2 contains links from V2 to t. Note that
each node [i, j] ∈ V1 results in exactly two arcs in A2 directing to the two locations i and j in
V2. Finally, the capacity is tij for arc (s, [i, j]) ∈ A1, hi− si for arc (i, t) in A3, and infinity
for each arc in A2. Note that all arc capacities are nonnegative. We then simply apply any
maximum flow algorithm to obtain a minimum s-t cut (S, T ), where S and T are the source and destination sets. The nodes in V2∩ S are the optimal set of facilities to build.
Proposition 2. Suppose that g(w) = w and si− hi< 0 for all i ∈ I. Let (S, T ) be a minimum
s-t cut of the graph G generated as above. An optimal solution x∗ for the problem in (1) can be obtained by setting x∗i = 1 if i ∈ V2∩ S and 0 otherwise.
Proof. For (S, T ) to be a minimum cut, it cannot contain any link in A2. Therefore, S cannot
3In practice, g(·) may be approximately linear when one is at the initial stage of building facilities when the
diminishing marginal benefit property is still weak. This is why in this study most of the efforts are devoted to the problem of a nonlinear g(·). The analysis for the linear case is included for completeness.
contain [i, j] ∈ V1 while either i ∈ V2 or j ∈ V2 is in T . Given this, the capacity of the cut is X [i,j]∈V1∩T tij+ X i∈V2∩S (hi− si) = X [i,j]∈V1 tij − " X [i,j]∈V1∩S tij+ X i∈V2∩S (si− hi) # = X [i,j]∈E tij − " X [i,j]∈E tijx∗ix ∗ j + X i∈I (si− hi)x∗i # .
Therefore, a minimum cut gives an optimal solution to our facility location problem.
In summary, our facility location problem is polynomially solvable if the marginal benefit is constant but NP-hard if it is diminishing. Below we will concentrate on the NP-hard case. In the next section, we design a polynomial-time heuristic algorithm for our problem and develop its worst-case performance guarantees in some special cases.
4
Algorithms and worst-case performance analysis
In this section, we will develop an algorithm for obtaining a feasible solution to our facility location problem defined in (1). We will introduce our proposed algorithm in Section 4.1, provide an illustrative example for it in Section 4.2, and prove its worst-case performance guarantees in some special cases in Section 4.3.
4.1 Approximation-relaxation-sorting-aggregation algorithm (ARSA)
Our algorithm has four major steps. First, regardless the functional form of the g(·) function, we approximate it by a kink function from above. The resulting problem can be transformed to an equivalent linear integer program. Second, we decompose the linear integer program into n subproblems (n is the number of candidate locations), relax the integer constraint, and solve the linear relaxation of each subproblems. Third, we construct a feasible solution to a subproblem through a simple sorting procedure. We finally aggregate all solutions from all subproblems to generate our solution to the original problem. Therefore, we call this algorithm the approximation-relaxation-sorting-aggregation algorithm (ARSA). Below we explain each step in details.
original g(·) function, ARSA first uniquely defines a kink function gK(w) = aw if w ≤ Ba B otherwise , (2)
where a = g0(0) is the slope of g(·) at x = 0 and B = g(P
i∈Isi+
P
[i,j]∈Etij) is the maximum
possible value of g(·) obtained by building all facilities. Note that an equivalent expression of gK(w) is gK(w) = min{aw, B}. It is clear that g(x) ≤ gK(x) for all x ∈ [0, B], i.e., gK(·) approximates g(·) from the above.4 We then replace g(·) by gK(·) to convert (Pori) to another nonlinear program.
Figure 1 provides a simple example of the approximation. Suppose that g(w) = 4(1−e−w/4) (plotted in Figure 1 as the solid curve). Suppose also that the maximum benefit that one may collect is 13, then the maximum possible revenue is B = g(13) ≈ 3.845. We then have a = g0(0) = 1 and gK(w) = min{w, g(13)}. In Figure 1, gK(w) is plotted as the dashed curve. We can see that it is indeed an approximation of g(w) from the above.
Figure 1: Approximating g(w) = 4(1 − e−w/4) by gK(w) = min{w, g(13)}
As gK(x) = min{ax, B}, the nonlinear objective function can now be linearized by intro-ducing a new variable p and two new constraints. We also replace xixj by a new variable yij
to eliminate products of decision variables. For any slope a, we can get an equivalent problem by setting the slope a to 1 and replace si and tij by asi and atij, respectively. Therefore, we
4
Aboolian et al. (2007b) adopt the same from-the-above approximation to establish their α-approximate algorithm. Unfortunately, while their approach works for their monotone objective function (which includes revenues only), it cannot be applied to our nonmonotone one (which includes revenues and costs).
normalize a to 1 without loss of generality. The equivalent linear integer program is max p −X i∈I hixi s.t. p ≤ B, p ≤X i∈I sixi+ X [i,j]∈E tijyij yij ≤ xi, yij ≤ xj ∀[i, j] ∈ E xi, yij ∈ {0, 1} ∀i ∈ I, [i, j] ∈ E.
We call this linear integer program (P ).
Relaxation. For the integer program (P ), ARSA then creates n subproblems by adding a constraint
X
i∈I
xi = k
to the kth subproblem, k = 1, 2, ..., n. We call the kth subproblem (Pk). The optimal solution
of (Pk) builds exact k facilities. Obviously, the best among the optimal solutions to (Pk),
k = 1, ..., n, is an optimal solution to (P ). Our strategy is to obtain a near-optimal solution for each (Pk) and show that the best among our near-optimal solutions for (Pk)s is indeed be
near-optimal to (P ). To generate an integer solution for (Pk), ARSA first relaxes the integer
constraints to obtain its linear relaxation, whose optimal solution can be found in polynomial time. Let xki be the value of xi in an optimal solution to the linear relaxation of (Pk).
Sorting. For each subproblem (Pk), ARSA now sorts candidate locations by xki in the
descending order. Let xki
1 ≥ x
k
i2 ≥ · · · ≥ x
k
in, where ties are broken arbitrarily, ARSA then
compares two solutions ¯x and ˆx, constructed in the following way. The first solution ¯x is constructed by simply choosing facilities i1, i2, ..., and ik. More precisely, we have ¯xi1 = · · · =
¯
xik = 1 and 0 otherwise. The second solution ˆx is constructed by choosing facilities i1, i2,
..., ik−1, and ik+1. By doing so, we have ˆxi1 = · · · = ˆxik−1 = ˆxik+1 = 1 and 0 otherwise.
ARSA compares the objective values achieved by the two solutions (with respect to the original function g(·), not gK(·)) and reports the better one to be a feasible integer solution xARSA-k to (Pk). By repeating this for all subproblems, we obtain n candidate solutions that are feasible
to (P ) as well as (Pori).
Aggregation. The last step is simple: Choose the best solution among these n candidates. Figure 2 sketches the flow of ARSA (when n = 4). Algorithm 1 is the pseudocode describing ARSA.
Figure 2: Basic flow of ARSA
4.2 An illustrative example
To better explain ARSA, in this section we provide a simple example for the illustration purpose. Suppose that a decision maker is given n = 3 locations to build facilities, as shown in Figure 3. The three locations are labeled as 1, 2, and 3. For the three locations, the stand-alone benefits are s1 = s2 = s3 = 3, construction costs are h1 = h3 = 1 and h2 = 2, and network benefits are
t12= t23= 2 and t13= 5. We also have g(w) = 4(1 − e−w/4), the one depicted in Figure 1, to
convert benefit w into revenue g(w). The optimal solution is to build facilities 1 and 3.
Figure 3: An illustrative example To solve this instance by ARSA, we do the four steps as follows.
Step 1: Approximation. As we have done in Figure 1, the kink function that we should adopt is gK(w) = min{w, 3.845} (cf. Figure 1). By replacing g(·) by gK(·), the original
problem (Pori) max
xi∈{0,1}
4 − e−(3x1+3x2+3x3+2x1x2+2x2x3+5x1x3)/4− (x
Algorithm 1 approximation-relaxation-sorting algorithm (ARSA)
1: Find the kink function gK(·) that approximate g(·). Replace g(·) by gK(·) to obtain the problem (P ).
2: Split the problem to n subproblems, one with an additional constraint P
i∈Ixi = k, k =
1, ..., n. Let the subproblems be (P1), ..., and (Pn).
3: for k from 1 to n do
4: Relax the integer constraints in (Pk).
5: Solve the relaxation of (Pk). Let xki be the value of xi in the optimal solution.
6: Sort locations so that xki
1 ≥ · · · ≥ x
k
in, where ties are broken arbitrarily.
7: Construct a solution ¯x such that ¯xi1 = · · · = ¯xik = 1 and 0 otherwise.
8: Construct a solution ˆx such that ˆxi1 = · · · = ˆxik−1 = ˆxik+1 = 1 and 0 otherwise.
9: if z(¯x) > z(ˆx) then
10: Report ¯x as the proposed solution xARSA-k for (Pk).
11: else
12: Report ˆx as the proposed solution xARSA-k for (Pk).
13: end if
14: end for
15: Report the best solution among xARSA-1, xARSA-2, ..., and xARSA-n.
can then be transformed to (P ) max p − (x1+ 2x2+ x3)
s.t. p ≤ 3.845, p ≤ 3x1+ 3x2+ 3x3+ 2y12+ 2y23+ 5y13
y12≤ x1, y12≤ x2, y23≤ x2, y23≤ x3, y13≤ x1, y13≤ x3
x1∈ {0, 1}, x2 ∈ {0, 1}, x3 ∈ {0, 1}, y12∈ {0, 1}, y23∈ {0, 1}, y13∈ {0, 1}.
Note that (P ) becomes a linear integer program.
Step 2: Relaxation. The second step is the relaxation step. We create three subproblems (P1), (P2), and (P3), where (Pk) is created by adding the constraint x1+ x2+ x3= k. We
More precisely, the relaxation of (Pk) is max p − (x1+ 2x2+ x3) s.t. x1+ x2+ x3 = k p ≤ 3.845, p ≤ 3x1+ 3x2+ 3x3+ 2y12+ 2y23+ 5y13 y12≤ x1, y12≤ x2, y23≤ x2, y23≤ x3, y13≤ x1, y13≤ x3 x1 ∈ [0, 1], x2∈ [0, 1], x3 ∈ [0, 1], y12∈ [0, 1], y23∈ [0, 1], y13∈ [0, 1].
Note that each of (Pk) is now just a linear program.
Step 3: Sorting. The third step is the sorting step. For each subproblem (Pk), we solve
its linear relaxation by any linear programming algorithm. Then we sort locations based on the values of xi in the obtained optimal solution and select either the first k locations
or the first k − 1 locations and the (k + 1)th location, depending on which solution gives us a higher profit (based on g(·), not gK(·)).
– For (P1), an optimal solution to its relaxation is x1 = (0.1207, 0, 0.8793), and the sorting results in the order 3, 1, and 2. Two candidate solutions to (P1) are then
¯
x = (0, 0, 1) and ˆx = (1, 0, 0). As z(¯x) = z(ˆx) = 1.1105, we break tie arbitrarily and report xARSA-1 = (0, 0, 1) as the solution to (P1).
– For (P2), an optimal solution to its relaxation is x2 = (1, 0, 1), and the sorting results in the order 1, 3, and 2 (ties are broken arbitrarily). Two candidate solutions to (P2) are then ¯x = (1, 0, 1) and ˆx = (1, 1, 0). As z(¯x) = 1.7443 and z(ˆx) = 0.4587, we report xARSA-2= (1, 0, 1) as the solution to (P2).
– For (P3), an optimal solution to its relaxation is obviously x3 = (1, 1, 1), the only feasible solution. We therefore report xARSA-3 = (1, 1, 1) as the solution to (P3), even though its associated profit z(xARSA-3) = −0.0444 is negative.
Step 4: Aggregation. We now compare the three candidate solutions (0, 0, 1), (1, 0, 1), and (1, 1, 1). As (1, 0, 1) results in the highest profit (1.7443 > 1.1105 > −0.0444), ARSA reports (1, 0, 1) as the proposed solution to the original problem (Pori). For this instance, ARSA finds an optimal solution.
4.3 Worst-case performance analysis of ARSA
Obviously, ARSA cannot always find an optimal solution. Therefore, it is worthwhile to in-vestigate how badly ARSA may perform. While we cannot analytically prove a worst-case performance guarantee for ARSA in general, in this section we will present two performance guarantees when ARSA is applied to two special cases of our facility location problem. The average performance of ARSA in the general case will be examined in Section 5.
To describe the first special case, we state three assumptions below.
Assumption 1. g(·) is the kink function gK(·) in (2) for a = 1 and some B > 0. Moreover, si< B for all i ∈ I.
This assumption restricts our attention to effective demand functions that are kink func-tions. In this case, the approximation step of ARSA can be omitted, and there is no precision loss due to the approximation. Moreover, as we mentioned above, a can be normalized to 1 without loss of generality. Note that under this assumption, it is without loss of generality to assume that there is no location whose si > B. If such a location exists, an optimal solution
will either contain only this location or does not contain it. We may thus safely remove this location for a while, solve the remaining problem, and at the end check whether selecting only this location is actually the best option.
Assumption 2. The net stand-alone benefits si− hi, i ∈ I, satisfy
max
i∈I {si− hi} > 0.
This assumption states that there is at least one facility which can bring profit itself. When this is true, a profitable option is to build only one facility. Under Assumption 2, we define a ratio that directly affects our worst-case performance guarantee.
Definition 1. For a given instance of our facility location problem in (1), we define its critical ratio as
r = −mini∈I{si− hi} maxi∈I{si− hi}
.
Note that r is readily available when an instance of our facility location problem is given. Also note that according to our definition and Assumption 2, r is negative if si− hi > 0 for all
i ∈ I or positive if there is at least one location whose net stand-alone benefit si− hi < 0. It is
Finally, we need all network benefits to be equal. One example in which the assumption may be valid is to build car/bike sharing stations where all drives/rides between any two locations are equally possible. The next assumption formalizes this requirement.
Assumption 3. There is a constant t ≥ 0 such that the network benefits tij, [i, j] ∈ E, satisfy
tij = t ∀[i, j] ∈ E.
Before proving that ARSA has a worst-case performance guarantee under our three as-sumptions, we first follow the idea of Caprara et al. (2000) to prove Lemma 1. It states that, under Assumption 3, an optimal solution to the relaxation of subproblem (Pk) has either zero
or two xk
is that are fractional.
Lemma 1. Suppose that Assumption 3 is satisfied. An optimal solution xLP −k to the relaxation of subproblem (Pk) has either zero or two fractional components. If there are two fractional
components, we have X i∈I sixLP −ki + k(k − 1) 2 t = B.
Proof. First, note that Assumption 3 and the additional constraint P
i∈Ixi = k together make
the subproblem (Pk) become
max p −X i∈I hixi s.t. p ≤ B, p ≤X i∈I sixi+ k(k − 1) 2 t X i∈I xi = k xi ∈ {0, 1} ∀i ∈ I,
where the total network benefit is exactly k(k−1)2 t regardless of which k locations are selected; the variable yij is therefore not needed anymore. For the linear relaxation of (Pk), the n + 1
variables require that at least n+1 constraints are binding at an optimal extreme point solution. If at least three xis are fractional, we will have at most n − 3 constraints binding in the set
of constraints xi ∈ [0, 1]. The maximum number of binding constraints is thus n, which is not
enough. If exactly two xis are fractional, the first three constraints must all be binding, and
the right-hand-side values of the first two constraints must be identical. BecauseP
i∈Ixi must
be an integer, there is no solution with only one fractional xi. Finally, no fractional variable is
Note that the relaxation of (Pk) is obtained by first adding the constraint Pni=1xi = k
and then relaxing the integer constraints, not the opposite. If we reverse the order, we will be unable to obtain the result in Lemma 1.
With Lemma 1, we now prove that ARSA is a 2+r1 -approximation algorithm for our facility location problem when our three assumptions hold.
Proposition 3. Suppose that Assumptions 1, 2, and 3 hold. For the problem defined in (1), let z∗ and z0 be the objective values of an optimal solution and the solution reported by ARSA, respectively. We then have zz∗0 ≥ 2+r1 , where r is defined in Definition 1.
Proof. Let zk∗ denotes the objective value of an optimal solution to subproblem (Pk) (with the
constraint P
i∈Ixi = k) and zkLP denote that to its relaxation. Note that
max k=1,...,n{z LP k } ≥ max k=1,...,n{z ∗ k} = z∗. (3)
Let xLP −k and xARSA−k be an optimal solution to the relaxation of (Pk) and the solution
reported by ARSA, our plan is to show that z0 = max
k=1,...,n{z(x
ARSA−k)} ≥ maxnz(xARSA−k), z(xARSA−1)o
≥ 1 2 + rz(x LP −k) = 1 2 + rz LP k (4)
for each k = 1, 2, ..., n, which implies that z0 ≥ 1
2+rmaxk=1,...,n{z LP
k }. Combining this and (3)
then completes the proof.
We now prove the second inequality in (4) for a given k. From Lemma 1, we know xLP −k has either two or no fractional value. If it has no fractional value, ARSA obviously selects an optimal solution to make z(xARSA−k) = z(xLP −k). Now suppose that there are two fractional variables. Without loss of generality, let the two variables be xLP −k1 and xLP −k2 such that s1 ≥ s2 and xLP −k1 = c = 1 − xLP −k2 for some c ∈ (0, 1). Let L0 be the set of k − 1 locations
with xLP −ki = 1. For the linear relaxation, from Lemma 1 we know zkLP = z(xLP −k) =X i∈I sixLP −ki + k(k − 1) 2 t − X i∈I hixLP −ki = X i∈L0 (si− hi) + c(s1− h1) + (1 − c)(s2− h2) + Tk, (5)
where the constant Tk = k(k−1)2 t is the total network benefit. ARSA will select the k − 1
former solution and L2 = L0∪ {2} be the latter. If location 2 is selected, because s1 ≥ s2, we have X i∈L2 si+ Tk < X i∈I sixLP −ki + Tk= B,
i.e., the total benefit obtained in the solution L2 does not exceed B. Therefore, we have
z(xARSA−k) ≥ X
i∈L0
(si− hi) + (s2− h2) + Tk, (6)
where the right-hand-side value is the objective value resulted from L2, which is just one of the
two candidates of xARSA−k. Combining (5) and (6), we then have the difference zkLP − z(xARSA−k) ≤ ch(s
1− h1) − (s2− h2)
i
< (s1− h1) − (s2− h2).
Obviously, s1 − h1 ≤ maxi∈I{si− hi}. Moreover, we have −(s2 − h2) ≤ r maxi∈I{si − hi}
according to Assumption 2. Therefore, we have z(xARSA−k) + (1 + r) maxi∈I{si− hi} ≥ zkLP.
As maxi∈I{si− hi} is the solution reported by ARSA for (P1), we have
z(xARSA−k) + (1 + r)z(xARSA−1) ≥ zkLP, which implies that max{z(xARSA−k), z(xARSA−1)} ≥ 2+r1 z(xLPk ).
Note that if there is at least one location such that si − hi < 0, we have r > 0 and thus
the performance guarantee 2+r1 < 12. On the contrary, if si− hi > 0 for all i ∈ I, we have
r < 0, which then implies that the performance guarantee is above one half. Roughly speaking, ARSA performs better (in the worst case) when the variability of net stand-alone benefit si− hi
decreases. Finally, note that if all si − hi are identical, one may easily evaluate the profit of
opening 1, 2, ..., and n facilities and find an optimal solution in polynomial time. When si− hi
are all the same, we have r = −1, which implies that the performance guarantee is 2+r1 = 1. In other words, ARSA does not fail to find an optimal solution for such a trivial case.
We next prove that ARSA is a max{12,2+r1 }-approximation algorithm for our facility loca-tion problem in another special case. We first state another assumploca-tion, which is more restricted than Assumption 3.
Assumption 4. The network benefits tij = 0 for all [i, j] ∈ E.
We now prove that ARSA has the performance guarantee max{12,2+r1 } under Assumptions 1 and 4. In other words, the performance guarantee is at least 12 but may be better if r < 0 (i.e., mini∈I{si− hi} > 0). Note that now Assumption 2 is not needed.
Proposition 4. Suppose that Assumptions 1 and 4 hold. For the problem defined in (1), let z∗ and z0 be the objective values of an optimal solution and the solution reported by ARSA, respectively. We then have zz∗0 ≥ max{12,2+r1 }.
Proof. Due to the similarity between this proof and the previous one, we only describe the different part. Because t = 0, the optimal solution to the linear relaxation does not contain any facility satisfying si− hi < 0. This immediately implies
z(xLPk ) − z(xARSA−k) ≤ ch(s1− h1) − (s2− h2) i < (s1− h1) + max{0, r} max i∈I {si− hi} ≤ (1 + max{0, r}) max i∈I {si− hi},
where the second inequality comes from the fact that s2− h2 ≥ 0. It then follows that
z(xARSA−k) + (1 + max{0, r})z(xARSA−1) ≥ zkLP and thus max{z(xARSA−k), z(xARSA−1)} ≥ max{1
2, 1
2+r}zLPk .
It is indeed unfortunate that the performance guarantees do not apply to our problem in general. Nevertheless, even if we only consider our algorithm as a heuristic algorithm for the general problem, we are still able to demonstrate its good average-case performance in Section 5 through numerical experiments. This shows the general applicability of our algorithm. Our analytical results in Section 4.3 and numerical results in Section 5 thus complement each other.
5
Numerical study
To understand how ARSA performs, we compare it with two other algorithms, an exact algo-rithm (either branch and bound or complete enumeration) and the genetic algoalgo-rithm. To test ARSA’s performance under different circumstances, the numerical study is done by incorporat-ing four factors. Below we will describe the experiment settincorporat-ing in Section 5.1, explain the exact and genetic algorithms in Section 5.2, and demonstrate the results of the numerical experiments in Section 5.3. A discussion about computation time is in Section 5.4.
5.1 Experiment setting
The first factor is the instance scale, measured by the number of candidate locations n, with three levels 20, 50, and 100. These three levels are labeled as small, medium, and large.
The second factor is the net stand-alone benefit si− hi. We first let hi ∼ U(1, 40), i.e.,
hi follows the uniform distribution between 1 and 40. We then consider two relationships
between si and hi: either si =
√
hi or si ∼ U(1, 40). The deterministic relationship si =
√ hi
is to generate cases where all locations have negative net stand-alone benefit.5 The alternative setting si ∼ U(1, 40) is to generate cases that these two parameters are unrelated. The net
stand-alone benefits may thus be of either sign. These two levels are labeled as negative and random.
The third factor is the distribution of network benefits. We consider three settings: tij = 0,
tij = t, where t ∼ U(0, 24), and tij ∼ U(0, 24). These three levels, which are labeled as zero,
identical, and random, are adopted to test the situations when Assumption 4 holds, Assumption 3 holds, or neither holds, respectively. The upper bound of tij when tij is not restricted to 0 is
chosen to be 24 intentionally. To see the reason, note that the expected value of si is
Z 40 1 √ x 1 39 dx ≈ 4.307 if si= √ hi or Z 40 1 x 1 39 dx = 20.5
if si ∼ U(1, 40). To generate both situations where the stand-alone benefits dominate network
benefits and the opposite, we should set the expected value of the network benefits to be around the average of 4.307 and 20.5. Therefore, we set the upper bound of tij to be 24.
Finally, the last factor is the shape of the g(·) function. The first shape we consider is a kink function, i.e., g(w) = min{w, K}. The second shape we consider is a smooth function, i.e., g(w) = K(1 − e−w/K). These two levels are labeled as kink and smooth. In either case, the parameter K is set so that it is optimal to build around 60% of facilities. Let µ be the expected value of si (which is 20.5 if si ∼ U(1, 40) or 4.307 if si=
√ hi), we set K = 0.6nµ if tij = 0 0.6nµ +(0.6n)(0.6n−1)2 t if tij = t 0.6nµ + 12 (0.6n)(0.6n−1) 2 if tij ∼ U(0, 24) .
The four factors together generate 36 scenarios. However, for the scenarios with si =
√ hi
and zero or identical tij, we will have hi > hj if and only if si− hi < sj− hj. In this case, it
5
When si =
√
hi, to make si− hi =
√
hi− hi < 0 for all i, we need hi > 1. This is why we set the lower
is always optimal to select facilities in the ascending order of hi until building one more facility
decreases the objective value. Therefore, we remove the corresponding twelve scenarios and discuss the remaining 24 ones below.
For each scenario, we generate 100 instances to test ARSA’s performance. Both the average performance across the 100 instances and the worst performance out of the 100 instances are recorded.
5.2 An exact algorithm and the genetic algorithm
To test the performance of ARSA, the first target for comparison is an exact integer solution (or its upper bound). Consider the situation when n = 20 first. For instances with a kink g(·), we linearize the piecewise-linear objective function and then apply CPLEX to find an optimal solution for the linear integer program. For those with a negative exponential g(·), we apply complete enumeration. Unfortunately, CPLEX and complete enumeration both fail to find an optimal solution when n = 50 and 100. In this case, we relax the integer constraints to generate a linear program if g(·) is kink or a convex program if g(·) is negative exponential. We then use CPLEX or MINOS to solve the linear or convex program respectively to obtain a deterministic upper bound of an optimal integer solution’s objective value.
We also compare ARSA with the genetic algorithm (GA), which is a meta-heuristic algo-rithm that is widely applied for solving NP-hard combinatorial problems (Hillier and Lieberman, 2014). We encode a solution to our facility location problem as a chromosome with n genes, one for each location. Gene i is 1 if a facility is built at location i and 0 otherwise. We adopt the objective value of a solution as the fitness value of a chromosome.
To start GA, we first generate 500 chromosomes to form the initial population. To randomly generate a chromosome, we first randomly select a probability p within 0.1, 0.2, ..., and 0.9. Then each gene is set to 1 with probability p. If a generated chromosome’s fitness value is negative, we redo the generation process until we obtain one with a positive fitness value. By doing so, all chromosomes in the population are guaranteed to have positive fitness values throughout the solution process.
With the initial population, we start to iterate. In each iteration, we select the two chromo-somes with the highest fitness values as parents. Then a number k ∈ {1, 2, ..., n − 1} is randomly generated to perform crossover. We generate two children, where child 1’s first k genes are from
parent 1, child 1’s last n − k genes are from parent 2, child 2’s first k genes are from parent 2, and child 2’s last n − k genes are from parent 1. Then each gene of each child is mutated to its opposite value with a mutation probability 0.02. Among the two parents and two children (after mutation), we select the best two to replace the two parents in the population. To conduct a fair comparison between ARSA and GA, we equalize the running times of GA and ARSA by having GA run 10000, 40000, and 100000 iterations for n = 20, 50, and 100.
The experiments are conducted on a personal computer with Windows 7, 12G RAM, and Intel i5-4570 3.2 GHz CPU. We use the Java programming language to implement complete enumeration, ARSA, and GA. The same language is used to invoke CPLEX and MINOS.
5.3 ARSA’s performance
Table 2 lists the complete result of our numerical experiments. The first four columns label the levels of the four factors. For instance scale, S, M, and L stand for small (n = 20), medium (n = 50), and large (n = 100); for stand-alone benefit, R and N stand for random (si∼ U(1, 40))
and negative (si =
√
hi); for network benefit, Z, I, and R stand for zero (tij = 0), identical
(tij = t, t ∼ U(0, 20)), and random (tij ∼ U(0, 20)); finally, for the effective demand function
g(·), K and S stand for kink (g(w) = min{w, K}) and smooth (g(w) = K(1 − e−w/K)).
For each of the 24 scenarios, we run ARSA and GA on the 100 cases, and for each algorithm we record the performance ratios zARSA/z∗ and zGA/z∗, where zARSA is the objective values of the solutions obtained by ARSA and zGA is that obtained by GA. Regarding z∗, it is the objective value of an optimal solution when n = 20 or an upper bound of that obtained by relaxing the integer constraints when n = 50 or 100. Each of the fifth and sixth columns then records the average of the 100 ratios for ARSA and GA, and each of the seventh and eighth columns records the minimum of the ratios. From Table 2, as the numbers in the fifth column are mostly close to 1, we conclude that ARSA typically obtains a near-optimal solution. Moreover, by comparing the fifth and sixth columns or the seventh and eighth columns, it is suggested that ARSA in general performs better than GA, both in the average and worst cases. To understand how ARSA performs under different scenarios, we conduct further summaries with respect to the four factors. The summaries are provided in Tables 3 to 6. Because the main reason of conducting the numerical study is to investigate the average-case performance of ARSA, below we will focus on the average of performance ratios. The worst-case performances are provided for completeness.
Instance Net stand-alone Network Effective Average Minimum scale benefit benefit demand zARSA/z∗ zGA/z∗ zARSA/z∗ zGA/z∗
S R Z K 0.983 0.758 0.891 0.620 S R Z S 0.971 0.817 0.807 0.664 S R I K 0.986 0.920 0.926 0.636 S R I S 0.980 0.925 0.877 0.618 S R R K 0.986 0.940 0.938 0.893 S R R S 0.922 0.970 0.832 0.936 S N R K 0.981 0.922 0.921 0.868 S N R S 0.857 0.947 0.714 0.868 M R Z K 0.991 0.504 0.952 0.373 M R Z S 0.623 0.333 0.533 0.224 M R I K 0.989 0.912 0.939 0.491 M R I S 0.836 0.811 0.452 0.264 M R R K 0.994 0.959 0.982 0.948 M R R S 0.897 0.895 0.876 0.877 M N R K 0.990 0.952 0.967 0.933 M N R S 0.889 0.889 0.865 0.867 L R Z K 0.993 0.349 0.961 0.258 L R Z S 0.621 0.220 0.556 0.149 L R I K 0.996 0.944 0.953 0.416 L R I S 0.877 0.846 0.462 0.206 L R R K 0.998 0.977 0.992 0.973 L R R S 0.926 0.915 0.923 0.906 L N R K 0.997 0.975 0.989 0.969 L N R S 0.930 0.918 0.924 0.907
First, we examine the impact of instance scale. By looking at Table 3, we realize that the good performance of ARSA is not prone to the instance scale. Note that while the performance ratios are calculated using an optimal solution when n = 20, they are calculated using an upper bound of an optimal solution when n = 50 or 100. This explains the drop in performance ratios from the small scale to the medium scale.
Instance Scale Average Minimum zARSA/z∗ zGA/z∗ zARSA/z∗ zGA/z∗ Small 0.958 0.900 0.714 0.618 Medium 0.901 0.782 0.452 0.224 Large 0.917 0.768 0.462 0.149
Table 3: Numerical results of problem scale
We next consider the impact of the relationship between the stand-alone benefit and con-struction cost. Table 4 shows that ARSA performs better when the relationship is deterministic rather than random. While this seems to be intuitive, we would like to emphasize that the deterministic setting si =
√
hi we set makes si− hi < 0 for all locations. As the net stand-alone
benefits are all negative, a consideration on network benefit is critical when searching for an op-timal solution. It is important to see that ARSA, which does take network benefit into account, performs well in this case. When si and hi becomes unrelated, naturally ARSA’s performance
becomes worse. Nevertheless, while the performance of GA drops significantly, that of ARSA only drops slightly. We conclude that ARSA performs much better than GA under the random relationship.
Net stand-alone benefit Average Minimum zARSA/z∗ zGA/z∗ zARSA/z∗ zGA/z∗
Random 0.921 0.778 0.452 0.149 Negative 0.941 0.934 0.714 0.867
Table 4: Impact of net stand-alone benefit
The impact of network benefit is illustrated in Table 5. Again, we see that ARSA outper-forms GA for all the three environments. It is somewhat surprising that ARSA peroutper-forms the worst when there is no network benefit. However, note that in average an ARSA-solution is above 85% as good as an optimal solution if the network benefits are identical and zero. This
is much better than the worst-case performance guarantees 2+r1 and max{12,2+r1 } obtained in Proposition 3 and 4. This allows decision makers to be more self-confident in applying ARSA to problems in practice.
Network benefit distribution
Average Minimum zARSA/z∗ zGA/z∗ zARSA/z∗ zGA/z∗ Zero 0.864 0.497 0.533 0.149 Identical 0.944 0.893 0.452 0.206 Random 0.947 0.938 0.714 0.867
Table 5: Impact of network benefit
Finally, we examine the impact of the effective demand function with Table 6. As the first step of ARSA is to use a kink function to approximate the true one, it is intuitive that ARSA’s performance is better when the function itself is a kink one. For a smooth g(·) function, we design ARSA by simply choosing a kink function to approximate the given function from above (cf. the definition of gK(·) in (2)). While this is naive and simple to implement, ARSA can still achieve above 85% of the solution quality of an optimal solution in average. Whether a better way of approximation exists deserves future investigation.
Function form
Average Minimum zARSA/z∗ zGA/z∗ zARSA/z∗ zGA/z∗ Kink 0.99 0.843 0.891 0.258 Smooth 0.861 0.791 0.452 0.149
Table 6: Impact of functional form
5.4 Computation time
To analyze the time complexity of ARSA, note that it is mainly governed by solving the linear relaxation of the subproblems. To solve a general linear program with m functional constraints by the simplex method, the time needed is roughly proportional to m3 (Hillier and Lieberman, 2014). For an n-location instance of our facility location, ARSA need to solve n LP problems, each has n + 3 functional constraints. Therefore, the time complexity of ARSA with the simplex method embedded is O(n4).
To verify the above time complexity analysis and provide a comparison between ARSA and finding an exact solution, we conduct another numerical study. For each level of instance scale n = 10, 20, ..., and 100, we consider the two levels of net stand-alone benefits and three levels of network benefits to form six scenarios. For each scenario, we again generate 100 instances. While we set up each instance by assuming a kink g(·) function, we exclude smooth g(·) functions from this numerical study to allow CPLEX to solve for an optimal solution. The average computation times per instance by each method for each instance scale are then recorded in Table 7 and plotted in Figure 4. The computation times of GA are intentionally set to be roughly the same as those of ARSA and are thus omitted.
n ARSA CPLEX 10 123.92 23.82 20 268.5 77.97 30 514.6 112.6 40 991.57 459.47 50 1936.12 4450.53 60 3585.87 49754.48 70 6574.62 – 80 11662.45 – 90 20846.22 – 100 30944.2 –
Table 7: Computation time comparison (in milliseconds)
From Table 7 and Figure 4, it can indeed be verified that the computation time of ARSA is proportional to n4 while that of CPLEX is proportional to 2n. In fact, if we fit these data into regression models, each with n4 and 2n as the independent variable and the computation time as the dependent one, the R2 values are both above 0.99. Though the O(n4) complexity of ARSA is not truly satisfactory, an instance with 100 candidate locations can be solved in around 30 seconds, and one with 300 locations is expected to be solved in around 40 minutes. As facility location problems are typically at the strategic level and do not need to be solved frequently, such a computation time should be acceptable in practice.
Figure 4: Average computing time of ARSA
6
Concluding remarks
In several service facility location problems such as building vehicle sharing sites, endogenous consumer demand and network effect play critical roles in finding an optimal way of building facilities. In this study, we capture these important effects into our facility location model. The solvability of our model depends on the shape of the concave effective demand function. When it is linear, we design a polynomial-time algorithm to find an optimal solution. When it is nonlinear, we show that it is NP-hard and design a polynomial-time heuristic algorithm to find a near-optimal solution. For two classes of special cases, our algorithm is proved to possess worst-case performance guarantees. Through a numerical study, our algorithm’s average-case performance is shown to be much better than the guarantees in general. It is also demon-strated that our algorithm can indeed find near-optimal solutions and outperforms the genetic algorithms in most cases.
To apply the model and algorithm in practice, practitioners need to put efforts in parameter estimation for si and tij (hj is the construction cost of building facility j and can be directly
traveling between locations i and j through mass transportation and private car. If the number is high, tij can be set high accordingly. Similarly, a questionnaire-based survey can help reveal
people’s willingness to use the two facilities at a pair of locations. In the questionnaire, one may provide an interviewee all the potential locations and ask her/him to select, say, five pairs of locations that she/he travels in between the most. By combining the answers from a sufficiently large number of interviewees, we will discover the potential needs of each pair of facilities and thus suggest the values for tij. If there is an old system existing (e.g., when one copies a
public bike system from one city to another city), the transaction data from the old system may certainly be used in estimating the attractiveness of the new system. The above approaches may all be applied to the estimation of si with some modifications. In fact, because one needs to
have both facilities built to travel in between, for transportation systems typically tij dominates
si. If a researcher has only limited budget/time to do parameter estimation for a transportation
system, we suggest to set si= 0 and focus on the estimation of tij only.6
There are several ways to extend this research. First, it is unclear whether there can still be a worst-case performance guarantee if the assumptions made in this study are relaxed. In particular, the case with different network benefits among locations is definitely worth of investigation. This study can serve as a good starting point. It is also a promising direction to find a way of approximating a smooth effective demand function by a kink function without losing the performance guarantee. Second, as our problem is only proved to be weakly NP-hard, a pseudopolynomial algorithm which uses dynamic programming may be designed to find an exact solution, at least for some special cases. Such a pseudopolynomial exact algorithm may complement our polynomial approximation algorithm. Finally, most (if not all) competitive facility location studies ignore the issue of parameter estimation, and most (if not all) demand estimation studies ignore the competitive nature among facilities. We hope to help connect these two streams in the future.
6Some recent works adopt various data analytics techniques to do the estimation; see, e.g., Hsieh et al. (2015),
Tiwari and Kaushik (2014), and the references therein. Interested readers may want to study these works. Nevertheless, most of these works only estimate the consumer demand at a location before a facility is built there. The endogeneity that new demand will be created after building a facility is widely considered difficult to estimate.
Acknowledgment
We thank the editor-in-chief Jos´e Oliveira and four anonymous reviewers for their detailed comments and many valuable suggestions that significantly enhanced the quality of this work. We also thank Zuo-Jun (Max) Shen for his helpful advising at the beginning stage of this work and Bertrand M.T. Lin and Kwei-Long Huang for their constructive comments at the final stage. Finally, we thank Po-Hsuan Chiang, Chien-Lin Chang, and Tzu-Hsiang Chien for helping us refining the numerical study. All remaining errors are our own.
References
Aboolian, R., O. Berman, D. Krass. 2007a. Competitive facility location and design problem. European Journal of Operational Research 182(1) 40–62.
Aboolian, R., O. Berman, D. Krass. 2007b. Competitive facility location model with concave demand. European Journal of Operational Research 181(2) 598–619.
Ageev, A., M. Sviridenko. 1999. An 0.828 approximation algorithm for uncapacitated facility location problem. Discrete Applied Mathematics 93(2-3) 149–156.
Ahmed, S., A. Atamt¨urk. 2011. Maximizing a class of submodular utility functions. Mathemat-ical Programming 128(1) 149–169.
Beresnev, V.L. 2014. On the competitive facility location problem with free choice of supplier. Automation and Remote Control 75(4) 668–676.
Berman, O., D. Krass. 1998. Flow intercepting spatial interaction model: a new approach to optimal location of competitive facilities. Location Science 6(1–4) 41–65.
Berman, O., D. Krass. 2002. Locating multiple competitive facilities: spatial interaction models with variable expenditures. Annals of Operations Research 111(1–4) 197–225.
Caprara, A., H. Kellerer, U. Pferschy, D. Pisinger. 2000. Approximation algorithms for knapsack problems with cardinality constraints. European Journal of Operational Research 123(2) 333– 345.
an analytic study of exact and approximation algorithms. Management Science 23(8) 789– 810.
Cornuejols, G., M. Fischer, G. Nemhauser. 1977b. On the uncapacitated location problem. Annals of Discrete Mathematics 1 163–178.
Daskin, M. S. 2013. Network and Discrete Location: Models, Algorithms, and Applications. Wiley, USA.
Feige, U., V.S. Mirrokni, J. Vondrk. 2011. Maximizing non-monotone submodular functions. SIAM Journal on Computing 40(4) 1133–1153.
Hillier, F., G. Lieberman. 2014. Introduction to Operations Research. 10th ed. McGraw Hill, USA.
Hsieh, H.-P., C.-T. Li, S.-D. Lin. 2015. Estimating potential customers anywhere and any-time on location-based social networks. The Proceedings of the 30th European Conference on Machine Learning and Principles and Practice of Knowledge Discovery in Databases (ECML/PKDD’15). Porto, Portugal, 576–592.
Karakitsiou, A. 2015. Modeling Discrete Competitive Facility Location. Springer, Germany. K¨u¸c¨ukayadın, H., N. Aras, ˙I.K. Altinel. 2011. Competitive facility location problem with
attrac-tiveness adjustment of the follower: a bilevel programming model and its solution. European Journal of Operational Research 208(3) 206–220.
K¨u¸c¨ukayadın, H., N. Aras, ˙I.K. Altinel. 2012. A leader-follower game in competitve facility location. Computers and Operations Research 39(2) 437–448.
Li, H.L., C.T. Chang, J.F. Tsai. 2002. Approximately global optimization for assortment prob-lems using piecewise linearization techniques. European Journal of Operational Research 140(3) 584–589.
Mel’nikov, A.A. 2014. Randomized local search for the discrete competitive facility location problem. Automation and Remote Control 75(4) 700–714.
Nemhauser, G.L., L.A. Wolsey. 1978. Best algorithms for approximating the maximum of a submodular set function. Mathematics of Operations Research 3(3) 177–188.
Nemhauser, G.L., L.A. Wolsey, M.L. Fisher. 1978. An analysis of approximations for maximizing submodular set functionsi. Mathematical Programming 14(1) 265–294.
Owen, S.H., M.S. Daskin. 1998. Strategic facility location: A review. European Journal of Operational Research 111(3) 423–447.
Tiwari, S., S. Kaushik. 2014. User category based estimation of location popularity using the road GPS trajectory databases. Geoinformatica: An International Journal 4(2) 20–31. Williamson, D.P., D.B. Shmoys. 2011. The Design of Approximation Algorithms. Cambridge
University Press, London, UK.
Wu, T.-H., J.-N. Lin. 2003. Solving the competitive discretionary service facility location prob-lem. European Journal of Operational Research 144(2) 366–378.
