European Journal of Operational Research 181 (2007) 598–619 www.elsevier.com/locate/ejor Discrete Optimization Competitive facility location model with concave demand Robert Aboolian a, Oded Berman b, Dmitry Krass b,* a b College of Business Administration, California State University San Marcos, San Marcos, CA 92096, USA Rotman School of Management, University of Toronto, 105 St. George Street, Toronto, Ont., Canada M5S 3E6 Received 1 March 2004; accepted 6 October 2005 Available online 18 September 2006 Abstract We consider a spatial interaction model for locating a set of new facilities that compete for customer demand with each other, as well as with some pre-existing facilities to capture the ‘‘market expansion’’ and the ‘‘market cannibalization’’ effects. Customer demand is assumed to be a concave non-decreasing function of the total utility derived by each customer from the service offered by the facilities. The problem is formulated as a non-linear Knapsack problem, for which we develop a novel solution approach based on constructing an efficient piecewise linear approximation scheme for the objective function. This allows us to develop exact and a-optimal solution approaches capable of dealing with relatively largescale instances of the model. We also develop a fast Heuristic Algorithm for which a tight worst-case error bound is established.  2006 Elsevier B.V. All rights reserved. Keywords: Location; Integer programming; Competitive facility location models; Non-linear Knapsack problem; Alpha-optimal solutions; Greedy heuristics; Worst-case bounds 1. Introduction The current paper has a dual objective: to introduce an interesting and important competitive location model, and to develop an efficient solution technique, which is applicable to a class of non-linear non-separable Knapsack models—which includes our competitive location model. Consider the problem of locating m new retail or service facilities controlled by the same entity (‘‘company’’) in an area where some other facilities providing comparable product or service may already exist. The new facilities, along with the pre-existing ones (some of which may be owned by the competitors, while others may belong to the company) will compete for customer demand with each other. We would like to find an optimal set of locations for the new facilities so that the total customers demand seen at all of the facilities controlled by the company is maximized. We assume that customer demand is concentrated at a set of n discrete points and there is a finite set of potential locations for the new facilities. * Corresponding author. Fax: +1 416 978 5433. E-mail address: krass@rotman.utoronto.ca (D. Krass). 0377-2217/$ - see front matter  2006 Elsevier B.V. All rights reserved. doi:10.1016/j.ejor.2005.10.075 R. Aboolian et al. / European Journal of Operational Research 181 (2007) 598–619 599 In order to be able to determine the total demand seen at a particular facility, it is first necessary to make certain modeling assumptions with respect to customer behavior, namely how do customers decide which of the available facilities to patronize and how much can they be expected to spend there. We adopt the general framework of ‘‘spatial interaction models’’, which originated with the work of Huff [18] (sometimes referred to as the ‘‘gravity’’ or ‘‘Huff-type’’ models in the location literature, and ‘‘brand choice’’ or ‘‘market share attraction’’ models in the marketing literature). Under this framework it is assumed that customer choice rules are probabilistic, i.e., that the customers split their demand between several (possibly all) of the available facilities, with the frequency (probability) of a visit to a particular facility decreasing with the distance from the facility and increasing with the ‘‘attractiveness’’ of the facility. We refer the reader to Berman and Krass [6,5], and Ghosh, McLafferty and Craig [15] for a detailed discussion of this class of competitive location models. The ‘‘standard’’ spatial interaction model assumes that customers’ expenditures are constant and known in advance—and are thus independent of the number and location of service facilities. This implies that in the presence of some existing facilities, the only effect of locating new facilities is the ‘‘cannibalization’’ of some of the market share by the new facilities from the existing ones. However, this assumption ignores the fact that in many locational decisions one of the key considerations is ‘‘market expansion’’—i.e., attracting new demand from customers that were ‘‘underserviced’’ by the existing facilities. This effect may cause the total customer demand in a given area to increase when new facilities are located there. In order to capture both the ‘‘cannibalization’’ and the ‘‘market expansion’’ effects, we follow the approach of Berman and Krass [5], who introduced the concept of ‘‘spatial interaction model with variable expenditures’’. Under this framework, the total expenditure (demand) of each customer is a function of the total utility derived by the customer from all of the available facilities (both new and pre-existing). While this approach does capture both the ‘‘cannibalization’’ and ‘‘market expansion’’ effects, it runs into severe technical difficulties, due to the complexity of the resulting optimization models when variable expenditures are introduced (in fact, even the constant expenditure case considered in [6] is computationally challenging). For the ‘‘general’’ form of the expenditure function considered in [5] (which assumed a ‘‘step function’’ form that can be used to approximate an arbitrary function), the only case that could be analyzed to any degree was that of the ‘‘monopolistic competition’’ (where the company owns all of the available facilities); even in that case the results were quite limited. In this paper we assume that the total demand functions (i.e., demand seen at each facility) are concave and non-decreasing functions of the total utility derived by each customer from the service offered by the facilities. We note that this appears to be a stronger assumption than simply assuming that the customer expenditure functions are concave and non-decreasing—since the demand observed at a particular facility also depends on the market share captured by this facility, which is itself a non-linear function of the utility derived by each customer. However, we show that it is, in fact, sufficient to assume that the customer expenditure functions are concave and non-decreasing functions of the utility derived by each customer. Concave non-decreasing demand assumption appears to be quite mild—it is a standard assumption in much of the marketing and economics literature—see, for example, Buzzell [7], Lodish, Montgomery and Webster [21], Little and Lodish [20]. However, we show that this assumption does come at a price—in some sense, it implies that the market expansion effect can never overtake the cannibalization effect. The concavity of the demand function allows us to formulate our model as a non-linear integer Knapsack problem with a particular form of the objective function. Non-linear Knapsack models are known to be difficult even with a separable objective function—see Hochbaum [16]; Bretthauer and Shetty [8] provide a good recent review. Unfortunately the Competitive Facility Location Model with concave demand falls within the sub-category of non-separable convex Knapsack problems, for which results are much more limited—see Pardalos et al. [27], Caprara et al. [9], and Dussault et al. [14]. An interesting related work was done by Klastorin [19] who considered a non-separable Knapsack model in the context of the media planning; in fact his problem is a special case of the class of problems for which we develop a general solution approach described below. The solution methodology developed in [19] is heavily dependent on the special functional form of his objective function and is not directly extendable to the significantly more general form considered in the current paper. 600 R. Aboolian et al. / European Journal of Operational Research 181 (2007) 598–619 Due to complexity of the objective function, we are not aware of any existing methods capable of solving moderate or large instances of our problem. This motivated us to develop a new approach that can be applicable to other problems as well. Our approach—which we call the Tangent Line Approximation (TLA)—is based on over-approximating the non-linear objective function with a piecewise linear function and is applicable to problems with any number of constraints (i.e., it contains Knapsack problems as a special case). While piecewise linear approximations are well-studied in the context of separable Knapsack problems, we could find no evidence of previous application of this approach for non-separable Knapsack problems. There are important differences between the TLA procedure and the standard linearization techniques in the optimization and numerical analysis literature. The standard approach, outlined in Bazaraa et al. [3] is to use a coarse grid initially and then to use a finer grid around the best solution obtained with the coarse grid—this procedure continues until the desired error has been approached (a variation of this approach is used in Meyer [22,23] where new grid points are generated only when necessary). This approach can lead to a large number of grid points—which increases the complexity of the approximating problem; moreover several approximation models may be required. In contrast, the TLA procedure employs an ‘‘adaptive grid’’ which is optimized to have the minimum number of grid points required to achieve a certain pre-specified level of the relative error; only one approximating model has to be solved. The use of the relative error is another distinguishing characteristic of the TLA—the standard numerical analysis approaches (e.g., see Cheney [10] or Nurenberger [26]) focus on defining an approximation scheme that achieves a certain level of the absolute error—which is of limited usefulness in our context (and, in fact, in most other optimization contexts, as absolute tolerance levels are hard to set a priori). The TLA procedure allows us to provide exact and a-optimal solutions to realistic-size instances of our competitive location model (by ‘‘a-optimal’’ we mean that the maximum allowable deviation from optimality a may be specified in advance by the user). We believe that the TLA procedure may have applications far beyond the specific location problem considered in the current paper. The plan for the remainder of the paper is as follows: we start by formally introducing the spatial interaction models with concave demand function in Section 2. Some structural results relating to model behavior, and an illustrative example of the model are presented there. The model is formulated as a non-linear Integer Programming model with concave objective function. In Section 3 we outline an efficient algorithmic approach (the ‘‘Tangent Line Approximation’’) for a general class of concave integer programming problems, which includes our problem as a special case (many technical details of the algorithm are relegated to the Appendix). In Section 4, we apply the TLA procedure to develop efficient exact and a-optimal solution approaches for our competitive location model. An alternative heuristic procedure is suggested in Section 5. While the maximum error for the heuristic cannot be chosen in advance, we do derive a tight worst case error bound. The results of a series of computational experiments aimed at testing the performance of the exact, a-optimal and the heuristics procedures are presented in Section 6. Some concluding remarks are presented in Section 7. To close the current section, we mention briefly some of the related work in competitive facility location. The problem of finding an optimal location for a single new facility in continuous space in the framework of (constant-expenditure) spatial interaction models is discussed in Drezner [13] and Plastria [28]. A corresponding multi-facility location model is investigated in continuous space using spatial interaction models in Drezner [12]. Following up on the initial work by Huff [18], discrete multi-facility location models of this kind have been at the heart of research in the retail sector (see e.g., [2,15]). In these papers the problem is modeled as a non-linear integer-programming problem, for which explicit enumeration schemes are proposed—these, of course, are applicable only to small-scale problems. Berman and Krass [6] developed a Branch and Bound scheme with a tight upper bound to optimally solve the problem efficiently for small size problems; they also provide an efficient heuristic approach for the constant-expenditure multi-facility case. We note that the TLA procedure can be applied to the constant-expenditure case as well, providing a significant improvement on the results in [6]—these will be reported in a separate paper. 2. Spatial interaction models with concave demand functions We assume that customer demand is concentrated at discrete points N = {1, 2, . . . , n} with point i containing wi customers. All customers residing at point i are assumed to be ‘‘identical’’ with respect to their facility R. Aboolian et al. / European Journal of Operational Research 181 (2007) 598–619 601 choice and expenditure decisions—thus, each point can be thought of as containing a separate customer ‘‘segment’’ (type). Let P  N be the set of potential facility locations—some of which are already occupied by the existing facilities. We assume that each location is capable of supporting at most one facility. This assumption is made without loss of generality, since a ‘‘replicate’’ of a location in P can always be added (e.g., if the actual physical location is capable of supporting two facilities, we can simply add a second copy of this location to P). We assume N to be a metric space, equipped with distance dij, i, j 2 N (e.g., N might consist of the nodes of a network, with d being the shortest path distance, or it may be embedded in a plane, with d being the l2 norm, etc.). Let E  P be the set of existing (competitive) facilities, and E ¼ P  E be the set of locations available for new facilities. To avoid indefinite forms later, we will assume that there is always some minimal travel distance  > 0 between any customer location and any facility, in other words wj = 0 for j 2 P. We will use S  E, jSj = m to denote the set of locations for m new facilities—the objective of our problem is to find S that optimizes the objective function introduced below. As mentioned in the introduction, we will use the general framework of gravity-type spatial interaction models, which were first introduced by Huff [18] to represent customer choice decisions. These models belong to the general class of stochastic utility models—also known as the ‘‘brand share models’’ in the marketing literature (see Naret and Weverbergh [25], and Nakanishi and Cooper [24]). A basic assumption in these models is that for a customer at point i 2 N, the utility of a facility at location j 2 S can be written as a ratio of a non-decreasing function of the ‘‘attractiveness’’ Aj of the facility, and a nondecreasing function of the travel distance dij, uij ¼ F ðAj Þ : H ðd ij Þ ð1Þ A typical form of uij is [18]: uij ¼ Aaj d bij a; b > 0: ð2Þ Another popular functional form for uij involves the use of exponential functions instead of power functions in (2) (see Beaumont [4]). In some applications it might be useful to define utilities so that uij = 0 unless j is within certain maximum distance from i, since a very attractive facility far away is unlikely to affect consumer’s behavior (in marketing terms, this maximum distance represents the radius of the ‘‘consideration set’’ for point i). Note that the utilities uij are independent of the locations of facilities other than j, and thus can be estimated a priori (i.e., before locational decisions are made)—data collection and estimation techniques are fairly welldeveloped for this class of models (see Nakanishi and Cooper [24]). For our purposes, the exact functional form of uij is not important—we simply treat these quantities as input parameters to our models. However, in the computational experiments reported in Section 5 below, we use the form (2) for generating values for uij. Let MSij be the frequency (probability) of customer at i using facility at j 2 S, where we assume, for the moment, that the location set S of new facilities is fixed. The quantity MSij can be interpreted as the share of the market available at demand point i that is captured by the facility at j. A fundamental assumption of spatial interaction models is that this frequency equals to the relative utility of facility at j compared to other facilities available on the network: uij : ð3Þ MS ij ¼ P k2E[S uik ‘‘Standard’’ spatial interaction models assume that the demand (i.e., total expenditure) of each customer for the service provided by the facilities is constant and known in advance. Let gi represent the demand (expenditures) of a customer at location i 2 N. Then, the total market size available at i is given by wigi. The problem of locating m new facilities can now be written as: XX wi gi MS ij : ð4Þ max SE i2N jSj¼m j2S 602 R. Aboolian et al. / European Journal of Operational Research 181 (2007) 598–619 This problem was analyzed in detail in Berman and Krass [6]. As mentioned in the introduction, a key weakness of the model described above is the assumption that expenditures gi are constant and independent of the number and locations of the facilities—thus, this model captures the ‘‘cannibalization’’, but not the ‘‘market expansion’’ effect. To overcome this limitation, Berman and Krass [5] generalized the model above by allowing the demand of customers at location i 2 N to be non-constant; it is assumed that this demand is an increasing function of the total utility obtained by a customer at i 2 N from all of the available facilities. Specifically, for i 2 N we define the demand (expenditure) function gi ðU 0i ðSÞÞ to be a non-negative function of U 0i ðSÞ, where X uij U 0i ðSÞ ¼ j2E[S is the total utility derived by customers at i from all of the facilities in E [ S. Without loss of generality, we can assume that 0 6 gi ðU 0i ðSÞÞ 6 1, and thus the total expenditure of customers at i 2 N for a particular location set S is given by wi gi ðU 0i ðSÞÞ. The problem of locating m new facilities with variable demand function can now be re-written as: XX max ZðSÞ ¼ wi gi ðU 0i ðSÞÞMS ij : ð5Þ SE jSj¼m i2N j2S Berman and Krass [5] discussed the linear, bounded linear and step-wise forms of the expenditure function. In general, the resulting optimization models are very difficult, due to inherent non-linearities of the objective function. For the general form of the model, only small instances can be solved to optimality. In this paper we will focus on the special case of model (5), where the demand function satisfies certain concavity assumptions. Specifically, we make the following definitions. For i 2 N, let: X X Hi ¼ uij ; and U i ðSÞ ¼ uij ¼ U 0i ðSÞ  H i : ð6Þ j2E j2S Note that these expressions represent the total utility which customers at i derive from all pre-existing facilities, and the total utility which customers at i derive from all new facilities, respectively. Note also that Hi are constant in our model for i 2 N. We also define X Hi MS i ðSÞ ¼ MS ij ¼ 1  ; ð7Þ H þ U i ðSÞ i j2S which represents the share of market i captured by all facilities in S (the new facilities). For U > 0 and e ¼ U þ H i , let U   e Þ ¼ gi ð U e Þ 1  Hi V iðU : ð8Þ U þ Hi e ¼ U i ðSÞ þ H i , the quantity V i ð U e Þ represents the total demand of a customer in market i 2 N captured For U by the facilities in S. Note that the expression for the objective function can now be rewritten as follows: X e i ðSÞÞ: ZðSÞ ¼ wi V i ð U ð9Þ i2N We make the following assumption which will remain in force for the remainder of the paper: e Þ is a concave, non-decreasing function of U e 2 ½H i ; 1 for all Assumption 1. The total demand function V i ð U i 2 N. The following result shows that this assumption is automatically satisfied when the expenditure function e Þ is concave, non-decreasing and vanishes at 0. gi ð U e Þ is a twice differentiable, non-decreasing concave Theorem 1. Suppose that the expenditure function gi ð U e Þ satisfies Assumption 1. function such that gi(0) = 0 for some i 2 N. Then, V i ð U R. Aboolian et al. / European Journal of Operational Research 181 (2007) 598–619 603 Proof. To simplify the notation in this proof we suppress the index i. By the definition of V(Æ) we have:   H 0 e 0 e e Þ H P 0; V ðU Þ ¼ g ðU Þ 1  þ gð U e e Þ2 U ðU where the inequality holds since both terms are non-negative. Also,   e Þ ¼ g00 ð U e Þ 1  H þ g0 ð U e Þ 2H  gð U e Þ 2H : V 00 ð U 2 e eÞ e Þ3 U ðU ðU e Þ3 from the last two The first term is clearly non-positive since g00 ( ) 6 0. Taking out the positive factor 2H =ð U e ÞU e  gð U e Þ. However, this last quantity is clearly non-positive—this can be verified by terms, we obtain g0 ð U e ÞU e 6 0 and it vanishes at 0. This implies that V00 ( ) 6 0, which completes observing that its derivative is g00 ð U the proof. h A convenient functional form for the expenditure function that satisfies the conditions of the previous result is gi ðU 0i ðSÞÞ ¼ 1  expðki U 0i ðSÞÞ; ð10Þ where ki > 0 is the exponential decay parameter (this parameter can be interpreted as the rate at which the total demand of a customer in market i reaches its saturation point wi). We will use this functional form in the computational experiments reported later in the paper. The assumption that demand is concave (equivalent to our Assumption 1) is quite standard in much of the economic literature. It will be seen below that Assumption 1 allows us to develop efficient exact and algorithmic approaches to our problem. However we stress that this assumption does come at a price. In particular, the following result, which follows directly from Theorem 1 in [5] (which is based on a result from [1]), shows that the ‘‘market expansion’’ effect is limited when Assumption 1 holds. TheoremP 2. Let S 2 P be the set of company’s facilities that have already been located. For a facility j 2 S let e i ðSÞÞMS ij be the demand currently seen by facility j. Suppose Assumption 1 holds. Then for Dj ðSÞ ¼ i2N wi gi ð U any k 2 P  S  E, Dj ðSÞ P Dj ðS [ fkgÞ; i.e., if a new facility is located at k, the demand seen by all existing facilities cannot increase. This results indicates that, from the point of view of individual facilities, the market expansion resulting from location of any new facilities will never overcome the cannibalization of existing demand by these facilities. We stress that this does not mean that the total demand seen by all of the facilities owned by the company (both new and already located) will not increase—indeed it can, and from the ‘‘company-wide’’ point of view the location of additional facilities may well be warranted. In spite of the previous result, it should be recognized that the location model described above, even when Assumption 1 holds, does provide a lot of flexibility to the decision-maker and can be applied to many reallife situations. For example, when the distance sensitivity parameter b in (2) is small, the facilities are likely to be centrally located and in close proximity to each other (under the optimal solution)—particularly if the demand elasticity is high. In practice, this may lead the decision maker to aggregate some of the facilities into a smaller number of larger units—‘‘big box’’ retail stores provide a good example of this situation. On the other hand, if the distance sensitivity b is high, the optimal solution will tend to spread out facilities resulting in a ‘‘service network’’ type of solutions—seen in practice in the case of, for example, bank machines. The following example illustrates how the demand cannibalization and market expansion effects are captured with our model. Example 1. Consider two customer demand nodes located at the endpoints of a unit line segment. We assume that 50 units of demand are available at each node (i.e., w1 = w2 = 50). We will use the exponential form of the expenditure function (10), set the distance sensitivity b = 2 and attractiveness Aj = 1 for all j in (2), and examine model behavior under the different values of the demand elasticity parameter k (here, the higher the value of k, the less elastic the demand). 604 R. Aboolian et al. / European Journal of Operational Research 181 (2007) 598–619 We will locate two facilities sequentially in this market: first pick the optimal location for facility 1, and then find the optimal location for facility 2 assuming the location of facility 1 is already fixed. We are interested in finding location patterns and measuring cannibalization and expansion effects. We will designate the location of facilities 1 and 2 by y1, y2—here yi represents the distance from the first demand point, and 1  yi— from the second, for i = 1, 2. First note that if the expenditures are constant (or, equivalently, k is high), then we have a ‘‘classic’’ Hotelling problem (see [17]), with the optimal location for the first facility at the midpoint (y1 = 0.5), and one facility capturing 100% (i.e., all 100 units) of the available demand. When the second facility is brought in, its optimal location is to co-locate with facility 1 (y2 = 0.5), it captures 50 units of demand and the total amount of captured demand stays at 100 units. Thus there is a 50% cannibalization effect (representing the decrease in the demand seen by facility 1 when the second facility is opened) and 0% market expansion effect (since the two facilities together intercept the same amount of demand as the first facility by itself). Note that these results will hold for any distance sensitivity parameter b—since the expenditures do not react to increased distance sensitivity. It can be verified numerically that the situation described above holds for the parameter values defined above when k > 1. However, the results are quite different for low values of k (corresponding to elastic demand). For example, for k = 0.1, in the one—facility case the location y1 = 0.5 is no longer desirable since the relatively high distance sensitivity coupled with large demand elasticity mean that the customers are spending just under 33% of their potential expenditures, with the total intercepted demand of about 33 units. We can increase this to 55 units by moving the facility very close to either endpoint (y1 = 0.86 is one optimal location; y1 = 0.14 is another)—this location receives 99.3% of the available potential demand from the ‘‘nearby’’ node, and only 12.6% from the ‘‘distant’’ one, for a total intercepted demand of 56 units. Suppose y1 = 0.86 is selected. The optimal location for the second facility is y2 = 0.06—i.e., this facility is ‘‘specialized’’ to node 1. The two facilities together intercept 99.8 units of demand—the market expansion is 78% in this case. However, the demand seen by facility 1 declines to 49.3 units—representing the cannibalization effect of 12%. Thus, in the case of elastic demand, the facilities are forced to focus on local markets—instead of a ‘‘shopping centre’’ solution in the constant expenditure case (with both facilities located in close proximity to each other), the new solution resembles a ‘‘corner store’’ pattern. Table 1 illustrates these effects for a range of the values of k (columns 2 and 3 contain the facility locations, while the next two columns contain cannibalization and expansion effects). It can be seen that as k gets smaller (i.e., as the demand sensitivity increases), the market expansion effect grows, the cannibalization effect decreases, and the distance between the two facilities grows. On the other and, for larger values of k we observe the same location pattern as in the constant expenditure case. We also note that an obvious extension of the model (5) is to make the constraint on the number of new facilities to be located flexible (by changing the constraint to jSj 6 m) and changing the objective to the profit maximization (rather than revenue maximization). We now turn our attention to formulating the location model defined above as a mathematical program. Let xj be a binary variable that is equal to 1 if locating a facility at node j and 0 otherwise. With this decision variables, we can formulate the location problem defined above as the following optimization problem. Table 1 Facility locations, cannibalization and expansion effects in Example 1 k y1 y2 Cannibalization (%) Expansion (%) 0.05 0.1 0.15 0.2 0.25 0.3 0.4 0.92 0.86 0.82 0.78 0.73 0.68 0.5 0.06 0.09 0.12 0.15 0.17 0.21 0.5 6 12 18 24 30 36 40 89 78 68 57 47 38 20 R. Aboolian et al. / European Journal of Operational Research 181 (2007) 598–619 max X i2N s:t: X 0 wi g i @ X j2E uij þ X j2E xj ¼ m; 1 uij xj A 1  P P j2E uij j2E uij þ P j2E uij xj ! 605 ðPO Þ j2E xj 2 f0; 1g; j 2 E: We will refer to problem PO as the ‘‘original model’’. This model is a non-linear IP (where the non-linearity comes from both, the expenditure function gi( ) and the ratio term in the objective) that contains jEj binary variables with one Knapsack-type constraint. By relaxing the integrality condition xj 2 {0, 1} to xj 2 [0, 1] we obtain continuous relaxation P RO . ^ ¼ ðx1 ; x2 ; . . . ; xjEj Þ be a feasible solution of PO. Clearly, there is a one-to-one correspondence between Let x ^ and a feasible location set S for the new facilities. Thus we will use the notations Ui(S), a feasible solution x Z(S) and U i ð^ xÞ, Zð^ xÞ interchangeably. Note that X xÞ ¼ U i ð^ uij xj ; ð11Þ j2E is a linear function of the decision variables. From (9) and Assumption 1 it follows that each term in the objective function of PO is a composition of a one-dimensional concave non-decreasing function and a linear function. In the following section we develop an efficient a—optimal procedure for integer programs with this type of objective function. 3. TLA procedure for a class of concave integer programming models Consider the following problem: max xðxÞ ¼ n X xi ð/i ðxÞÞ ðProblem DO Þ i¼1 s:t: AðxÞ 6 b where x non-negative, integer-valued vector in Rm, where for i = 1, . . . , n; xi(y) is a twice differentiable, nondecreasing, concave function of y 2 R+, with xi(0) = 0; /i ðxÞ ¼ cTi x; i ¼ 1; . . . ; n is a linear functional with non-negative coefficient vector ci P 0, and A(x) 6 b is a system of linear constraints. Moreover, we assume  i , for each i, where /  i is a constant. Note that one possible sufficient condition under which that 0 6 /i ðxÞ 6 / this assumption holds is that the feasible region of DO is bounded. We will denote the continuous relaxation DO by DRO . Observe that each term in the objective function x(x) is a composition of a concave non-decreasing function with a linear functional. As noted at the end of Section 2 above, the competitive facility location model PO defined in the previous section is a special case of DO. We note that if there is only one linear constraint (as is the case in PO), DO is belongs to the class of nonlinear integer Knapsack problem with non-separable objective function—which are notoriously difficult to solve, as discussed in the introduction. While it may appear that the problem can be transformed to a separable non-linear Knapsack problem by a simple variable substitution (defining new decision variables yi = /i(x)), this approach does not work since the argument of the objective function is no longer integer and it may be impossible to express the constraint in terms of the new variables, which means that the transformed model no longer has the Knapsack structure (see [8] for a nice discussion). The method we develop below allows us to solve DO directly in an efficient manner and with arbitrary number of linear constraints. In our opinion, the importance of the optimization models satisfying the same properties as DO extends beyond the competitive facility location model PO defined in the previous section. For example, the ‘‘classical’’ spatial interactive model analyzed in [6] also belongs to this class. Applications outside the field of locational analysis may also exist. Our ultimate goal in this section is to develop an approximate Linear Program DRa whose optimal solution lies within a pre-specified error bound a from the solution of the relaxed non-linear problem DRO . We will then 606 R. Aboolian et al. / European Journal of Operational Research 181 (2007) 598–619 show how this approximating LP can be used to find exact and a-optimal solutions to DO within a standard Branch-and-Bound scheme. Suppose the objective function x(x) is replaced by a concave piecewise linear function xA(x), which bounds x(x) from above. More specifically, for i 2 N, let xAi ð/Þ be a concave piecewise-linear function such that  i; xA ð/Þ P xi ð/Þ P 0 for 0 6 / 6 / ð12Þ i where xi(Æ) is the ith term in x. Note that the existence of xAi ð/Þ follows from the fact that xi(/) is bounded on  i . The following result now follows. ½0; / Lemma 1. For i = 1, . . . , n let  A  xi ð/i Þ  xi ð/i Þ ai ¼ max i xi ð/i Þ 06/i 6/ ð13Þ be the maximum relative error for approximating xi(/i) by xAi ð/i Þ. Let a0 ¼ maxfai g: i2N ð14Þ Then, for any feasible solution x of DO (or DRO ), xðxÞ 6 xA ðxÞ 6 ð1 þ a0 ÞxðxÞ ð15Þ Proof. Follows directly from (12) and (14). h Suppose we want to make sure that the maximum relative error from such an approximation does not exceed some pre-specified level a > 0. Then this requirement can be met by constructing concave piecewise linear functions xai for i = 1, . . . , n that satisfy (12) and for which ai < a holds in (13). The procedure for constructing such a piecewise linear function will be developed in Section 3.1 below.  For Suppose xai ð/i Þ with the property described above has been constructed for i = 1, . . . , n, and /i 2 ½0; /. m x 2 R . Let n X xai ð/i ðxÞÞ; xa ðxÞ ¼ i¼1 and define mathematical programs DA and DRA that have xa as the objective function and the same constraint sets as DO and DRO , respectively. Note that since xa is piecewise linear and concave by construction, DA (DRA ) can be re-formulated as a linear integer program (linear program) using standard formulation techniques for problems with concave piecewise linear objectives. The following results follows immediately from the preceding lemma. Lemma 2. For any a 2 [0, 1], the optimal solution of DA is a-optimal for DO and the optimal solution of DRA is aoptimal in DRO . It follows that by applying the preceding result to our original problem PO, it will be possible to obtain aoptimal solutions by solving a linear integer problem instead of a non-linear one (provided the concave piecewise linear over-approximating functions xai can be efficiently constructed, of course). We will come back to this issue in Section 4. We also note that in constructing our piecewise linear approximator, it is desirable to minimize the number of linear segments (for a given accuracy a). This is because each linear segment requires a separate decision variable in the LP representation of a piecewise linear function. It will be shown that the TLA procedure developed below is ‘‘optimal’’ in this regard—the piecewise linear approximator it constructs has the smallest possible number of linear segments for a given relative error a. 3.1. Construction of piecewise linear over-approximator for a given relative error a In this section we develop the fundamental results allowing us to construct the piecewise linear functions xai ð/i Þ possessing the properties described above for a given a 2 (0, 1) for i = 1, . . . , n. To simplify the notation R. Aboolian et al. / European Journal of Operational Research 181 (2007) 598–619 607 in this section, we will drop the subscript i. For easy reference, we re-state the basic problem addressed in this section.  Piecewise linear over-approximator problem. Let x(/) be a concave, non-decreasing function for / 2 ½0; /. Assume that x(/) is twice-differentiable and that x(0) = 0. For a given a 2 (0, 1) we would like to construct a concave, non-decreasing, piecewise linear function xa(/) such that  ð16Þ xa ð/Þ P xð/Þ for / 2 ½0; /;  a  x ð/Þ  xð/Þ  6 a; j0 < / < / ð17Þ max / xð/Þ and xa(0) = 0. The next result, stated in Theorem 3, shows that for any piecewise linear function that satisfies (16), in order for the property (17) to be satisfied it is sufficient that it be satisfied at the ‘‘breakpoints’’. The Proof of Theorem 3 can be found in the Appendix.  such that Theorem 3. Let xA(/) be a non-decreasing concave piecewise linear concave function on ½0; /  xA ð/Þ P xð/Þ for / 2 ½0; /. Let (/) = [xA(/)  x(/)]/x(/) be the relative error (distance) between xA(Æ) and x(Æ). Then (/) is maximized at the beginning of one of the line segments of the approximating function xA(/). Theorem 3 implies that if we ensure that (/) 6 a holds at the starting points of all line segments defining our piecewise linear over-approximator, then the relative error will indeed be equal to a. In fact, it will be shown with that in order to minimize the number of line segments used it is advantageous to set (/) = a at each starting point, except for the first one, where (0) can be set to 0. This observation forms the basis of the TLA procedure. In defining our piecewise linear function xa, the following notational conventions will be used: (1) we will index each line segment by l, with l = 1, . . . , L, (2) bl will be used to represent the slope of segment l, and (3) the  Thus, segment l is given by [cl, cl+1] for starting point of segment l will be represented by cl, with cLþ1 ¼ /. l = 1, . . . , L. Intuitively, the TLA procedure operates as follows (the technical details of how different steps are performed, and the correctness of the procedure are provided in the Appendix): Step A: We start by setting xa(0) = 0 and use the point 0 as the starting point of the first segment—i.e., set c1 = 0. Note that the relative error at 0 is set to 0 by convention (since x(0) = 0). We set the slope of the first segment equal to the first derivative of x(/) at 0, i.e., set b1 = x 0 (0). Thus, the segment is tangent to the graph of x(Æ) at 0. Step B: We find the endpoint of the current segment by locating a point wa(c2) on the ray originating at 0 and with slope b1 such that the relative error (c2) = a. This provides the starting point for the next segment. Step C: To find the slope of segment l for l = 2, . . . , L, we find the slope of the ray originating at (cl, wa(cl))  that is tangent to the graph of x(Æ) on ½cl ; /. Step D: To find the endpoint of segment l for l = 2, . . . , L, we repeat step B.  The following The procedure continues until the xa(/) function has been defined for all points in ½0; /. result establishes the convergence of the TLA procedure and provides an upper bound on the number of segments in the resulting piecewise linear approximator (which is the same as the number of major iterations of the TLA). The proof can be found in the Appendix. Theorem 4 1. For any a > 0, the TLA procedures converges in finitely many steps to a piecewise linear non-decreasing cona ð/Þ  cave function xa(/) such that 1 6 xxð/Þ 6 1 þ a for all / 2 ½0; /. a 2. Let  > 0 be the length of the first linear segment in x (/). Then the total number of segments (which is the  same as the number of major iterations in the TLA) is bounded above by 1 þ d/ e. 608 R. Aboolian et al. / European Journal of Operational Research 181 (2007) 598–619 We note that the upper bound on the maximum number of ‘‘major iterations’’ in the TLA procedure derived in the preceding result is very rough. In our experience, the procedure tends to require only a small number of iterations even for small values of a. This will be illustrated through a set of computational experiments below. However, we first discuss an interesting ‘‘optimality’’ property of the TLA procedure. Recall that our intent is to use the concave piecewise linear over-approximator xa in place of x in the mathematical programming formulations. Since each segment in xa will require a separate linearization (i.e., several new decision variables will have to be added to the formulation), it is desirable that xa have as few segments as possible. The next result shows that the function constructed by the TLA procedure is, in fact, optimal in this regard. Theorem 5. For a given a > 0, the TLA procedure constructs a concave piecewise linear function, which approximates x(/) from above within relative error a with the minimum number of line segments.  Proof. Let [cl,cl+1] for l = 1, . . . , L be the lth segment in xa(/), with c1 = 0 and cLþ1 ¼ /. A Let x (/) be another concave piecewise linear function such that  xð/Þ 6 xA ð/Þ 6 ð1 þ aÞxð/Þ for / 2 ½0; / ð18Þ Assume xA(/) consists of LA segments, with segment l given by ½cAl ; cAlþ1 . To prove the result we need to show that L 6 LA. Consider l 6 L  1 and let cT 2 [cl,cl+1] be the point where xa(/) is tangent to x(/), i.e., xa(cT) = x(cT) and bl = x 0 (cT), where bl is the slope of segment l in xa(/). Consider segment l in xA(/) with endpoints cAl ; cAlþ1 and slope bAl . First we establish the following preliminary result: cAl 6 cl implies that cAlþ1 6 clþ1 : Suppose, to the contrary, that cAl 6 cl but cAlþ1 > clþ1 . Then ½cl ; clþ1   ½cAl ; cAlþ1 , and cT 2 ½cAl ; cAlþ1 Þ. By assumption, xA(cT) P x(cT) = xa(cT). We distinguish between two cases: Case 1: xA(cT) > xa(cT). By linearity of xA(/) we have xA ðcT Þ ¼ xA ðcl Þ þ bAl ðcT  cl Þ > xa ðcT Þ ¼ xa ðcl Þ þ x0 ðcT ÞðcT  cl Þ: ð19Þ a However, by construction of x we have xa ðcl Þ ¼ ð1 þ aÞxðcl Þ P xA ðcl Þ; which contradicts (19) if bAl 6 x0 ðcT Þ. It follows that bAl > x0 ðcT Þ must hold in this case. Now, xA ðclþ1 Þ ¼ xA ðcT Þ þ bAl ðclþ1  cT Þ > xa ðcT Þ þ x0 ðcT Þðclþ1  cT Þ ¼ xa ðclþ1 Þ ¼ ð1 þ aÞxðclþ1 Þ; which contradicts (18). Case 2: xA(cT) = x(cT) = xa(cT). By linearity of xA(/) on ½cAl ; cAlþ1 , it follows that bAl ¼ x0 ðcT Þ ¼ bl . Thus, xa(/) = xA(/) for all / 2 [cl, cl+1]. Recall that by the TLA procedure and Lemma 4 (see Appendix), the point cl+1 is the unique root of xðcT Þ þ x0 ð/  cT Þ ¼ ð1 þ aÞxð/Þ;  and that (1 + a)x(/) < x(cT) + x 0 (cT)(/  cT) for / > cl+1. But, since cA > clþ1 , this on ½cT ; /, lþ1 implies that xA ðcAlþ1 Þ ¼ xA ðcT Þ þ bAl ðcAlþ1  cT Þ ¼ xðcT Þ þ x0 ðcT ÞðcAlþ1  cT Þ > ð1 þ aÞxðcAlþ1 Þ; which contradicts (18). It follows that cAl 6 cl implies that cAlþ1 6 clþ1 for all l = 1, . . . , L  1. Since c1 ¼ cA1 ¼ 0 this immediately implies that cAL 6 cL . It follows that LA P L must hold, thus establishing the result. h 609 R. Aboolian et al. / European Journal of Operational Research 181 (2007) 598–619 Table 2 Number of line segments for approximation of each customer node found by TLP ki 0.1 0.1 0.1 0.1 0.5 0.5 0.5 0.5 1 1 1 1 10 10 10 10 Hi 1 2 5 10 1 2 5 10 1 2 5 10 1 2 5 10 ai = 5% ai = 2% ai = 1% ai = 0.1% 4 7 9 24 4 7 8 23 4 6 8 22 2 4 6 15 5 8 12 37 5 8 12 35 5 8 11 32 4 6 8 23 5 9 12 37 5 9 12 36 5 8 11 33 4 6 8 25 5 9 12 38 5 9 12 37 5 9 12 33 4 7 9 26 As mentioned earlier, the usefulness of the TLA procedure for the purposes of solving the optimization problem DO depends critically on having a relatively small number of linear segments in the approximating function xa(/). To evaluate the efficiency of the TLA procedure in this regard we conducted a set of computational experiments, the results of which are presented in Table 2. For x(/) we used the exponential expenditure function (10) in the competitive location model described in the previous section (we looked at the component corresponding to customer i 2 N, with i chosen arbitrarily). We examined how the number of segments in the approximating function xa(/) depends on the required relative error rate a for a = 5%, 2%, 1% and 0.1%, as well as on several other model parameters: pre-existing facilities H = Hi—as defined the exponential decay factor k = ki and the total utility derived from  was set to /  ¼ P uij  H i —this is the upper limit on Ui in (6). The upper limit of the relevant interval / j62E in (8). Table 2 shows that the most critical parameter in determining the required number of segments is the relative error rate a, with the number increasing sharply as a becomes smaller. The required number of segments also tends to increase with the exponential decay parameter k, but this effect appears to be weaker than that of the relative error rate. This effect is likely due to the fact that the exponential function becomes steeper as k increases. Finally, the required number of segments decreases with H—this is likely due to the fact H is inver sely related to /. Overall, Table 2 shows that the number of segments required to achieve a reasonable level of accuracy is not large—for example, 12 segments are sufficient to reduce the relative error to below 1%. It appears that the TLA procedure provides a practical approach for solving the competitive location problem P0 defined in (11). This approach is explained in greater detail in the following section. 4. Applying TLA to the competitive location model P0 We begin by developing a linear ‘‘a-approximating’’ model for our original model P0 given in (11). a For i 2 N, let Vi(U), given in (8), be the ith term of the objective function of P0, and P let xi ðU Þ be the a-approximator of Vi(U) constructed by the TLA procedure, with U 2 ½0; U i  and U i ¼ j62E uij . Let Li be the number of linear segments in xai ðU Þ, with the endpoints of segment l given by cil and cilþ1 for l 2 {1, . . . , Li}. Let bli be the slope of segment l, and ail ¼ cilþ1  cil be the length of this segment (projected onto the X axis). The function xai ðU Þ can be represented as follows for U 2 ½0; U i : xai ðU Þ ¼ i ðU Þ LX ail bil y il ; l¼1 where Li ðU Þ ¼ maxfl : cl 6 U g (note that Li ¼ Li ðU Þ), and ( 1 if l < Li ðU Þ l y i ¼ U cil if l ¼ Li ðU Þ: ai l 610 R. Aboolian et al. / European Journal of Operational Research 181 (2007) 598–619 We now introduce the following linear integer program P a0 : a Z ¼ max Li XX i2N X s:t: wi ail bil y il uij xj ¼ Li X ail y il for i 2 N ; l¼1 j2E X ð20Þ l¼1 xj ¼ m; j2E 0 6 y il 6 1 for i 2 N ; l ¼ 1; . . . ; Li ; xj 2 f0; 1g for j 2 E: We will denote the linear relaxation of P a0 by P aR and use Za(x) to represent the objective value of P a0 for a particular feasible location vector x (with the y variables chosen optimally for this x). It should be clear that for any feasible location vector x, X Z a ðxÞ ¼ xa ðU i ðxÞÞ; i2N P where U i ðxÞ ¼ j2E uij xj , as in (11). The following result now follows from (9), Lemma 1, Theorem 4, and the fact that any location vector x is feasible in problem P0 if and only if it is feasible in problem P a0 . Theorem 6. Let x* be the optimal location vector in problem P0 (problem PR). For some a > 0, let xa be the optimal location vector in problem P a0 (problem P aR ). Then Zðx Þ 6 Z a ðxa Þ 6 ð1 þ aÞZðx Þ and 1 Zðx Þ 6 Zðxa Þ 6 Zðx Þ: 1þa The preceding result suggests two different ways to apply the TLA procedure to the competitive location model P0: Exact procedure: We apply a ‘‘regular’’ branch-and-bound algorithm, using the relaxed approximator problem P aR to obtain upper bounds, and the objective function of the origin non-linear Integer Program P0 to obtain lower bounds. We call this approach the ‘‘exact procedure’’ because the resulting solution is guaranteed to be optimal for the original problem P0—no matter what value of a is used (here smaller values of a result in tighter upper bounds at the cost of increasing the dimensionality of the LPs that have to be solved to obtain these bounds). a-Optimal approximate procedure: Under this procedure we solve the approximator problem P a0 directly and then use the resulting optimal location vector xa as an a-optimal solution in P0. The main advantage of this approach is that the problem P a0 , which is a linear Integer Program, can be solved by any of the ‘‘off-the-shelf’’ IP solvers—it is not necessary to develop a special-purpose solution algorithm. The obvious disadvantage is that the resulting solution is only guaranteed to be a-optimal for the original model. Both of the procedures described above will be explored further in our computational experiments described in Section 6 below. However, we point out that for practical problems both procedures involve solving large-scale Integer Programs, which may not be possible if the problem dimensionality is sufficiently large. Thus there is a need for an efficient heuristic procedure capable of providing good solutions quickly. Such a procedure is developed in the following section. 5. Heuristic procedure In this section we develop a greedy heuristic for our competitive location problem. This algorithm is much faster than either the ‘‘exact’’ or the ‘‘a-optimal’’ procedures described in the previous section, but, of course, it 611 R. Aboolian et al. / European Journal of Operational Research 181 (2007) 598–619 does not provide any direct control over the size of the relative error. However it is possible to obtain a bound on the maximum relative error that can result from the application of this greedy heuristic. Heuristic algorithm (the greedy heuristic): Starting with the empty set of facilities, find the potential location that results in the greatest increase in the objective function and locate a new facility there. Proceed until all m facilities have been located. Note that the running time for this algorithm is OðN jEjmÞ. The worst case analysis of the algorithm is based on the following Theorem, whose proof can be found in Appendix. Theorem 7. The objective function Z(S) of the problem P0 (defined in (9) above) is non-decreasing and submodular. The following corollary now follows immediately by the results of [11]. Corollary 1. Let Z H by the Heuristic Algorithm. Let Z* be the optimal O be the objective function value produced ZH O value of the objective function in original problem, then Z  P 1  1e  :63. The computational results in the next section show that the average performance of the Heuristic Algorithm above is much better than the performance guaranteed by the worst case bound. In fact, the relative errors are typically well below 2%. Thus, the Heuristic Algorithm above provides an efficient way to obtain high-quality approximate solutions for our competitive location model. On the other hand, the bound established in the preceding result is tight, as shown by the following construction. Consider an instance where we are to locate m facilities, E = ; (there are no pre-existing facilities), n = 3m, with nodes 1, . . . , m being the customer nodes (i.e., having positive demand), and nodes m + 1, . . . , 3m serving as potential facility locations (i.e., jEj ¼ 2m). For convenience, we will use index j = 1, . . . , 2m to refer to the potential locations. We assume that the expenditure function gi( ) is exponential—as defined in (10) above. Let wi = ki = 1 for i = 1, . . . , m. The utility vector uij for i 2 {1, . . . , m} and j 2 {1, . . . , 2m} is given by 8 > <  lnð1  1=mÞ for 1 6 i; j 6 m; uij ¼  lnðÞ for 1 6 i 6 m; j ¼ i þ m; > : 0 otherwise; where  is chosen so that 0 6  6 (1  1/m)m. For example, for the case of m = 4, the utilities are given by: 2  lnð0:75Þ  lnð0:75Þ  lnð0:75Þ  lnð0:75Þ  lnðÞ 0 0 6 6  lnð0:75Þ  lnð0:75Þ  lnð0:75Þ  lnð0:75Þ 0  lnðÞ 0 6 6 6  lnð0:75Þ  lnð0:75Þ  lnð0:75Þ  lnð0:75Þ 0 0  lnðÞ 4  lnð0:75Þ  lnð0:75Þ  lnð0:75Þ  lnð0:75Þ 0 0 0 0 0 0  lnðÞ Let A = {1, . . . , m} and B = {m + 1, . . . , 2m}. We make the following observations. 3 7 7 7 7: 7 5 Proposition 1. Let S* represent the optimal location set, and Sh represent the location set selected by the Heuristic Algorithm above. Then 1. S* = B and z(S*) = m(1  ). 2. Sh = A and z(Sh) = m(1  (1  1/m)m). Proof 1. Consider an arbitrary solution S, consisting of k nodes from set A and m  k nodes from set B, k 6 m. It is easy to check that the objective function value of this solution is given by k k kð1  exp½lnð1  1=mÞ Þ þ ðm  kÞð1  exp½lnð1  1=mÞ þ lnðÞÞ; 612 R. Aboolian et al. / European Journal of Operational Research 181 (2007) 598–619 a quantity which is clearly decreasing in k. Thus, the optimal solution must have k = 0 and the first statement follows. 2. The proof follows by induction. By the definition of , it is clear that the first facility will be located in A by the heuristic. Suppose that the algorithm locates the first i facilities in A. Now, if the i + 1-st facility is located in B, the resulting objective function value is ðm  1Þð1  exp½lnð1  1=mÞ i1 Þ þ ð1  exp½lnð1  1=mÞ i1 þ lnðÞÞ; which increases with . By setting  to its maximum value (1  1/m)m, and substituting above we obtain ðm  1Þð1  1=mÞ i1 þ ð1  1=mÞ mþi1 i 6 mð1  1=mÞ ; where the quantity on the right can be recognized as the objective function value obtained by locating the ith facility in set A. It follows, that at each step of the algorithm, it is beneficial to locate the next facility in set A, which completes the proof. h From the preceding result it follows that Z GO 1  ð1  1=mÞm ; ¼ Z 1 which approaches 1  (1  1/m)m as  ! 0. Thus, by selecting a large value of m and a small value of , we can approach the bound in Corollary 1 arbitrarily closely. 6. Computational results We have conducted a series of computational experiments to test the efficiency of the three computational approaches developed above, namely the ‘‘exact’’ and ‘‘a-optimal’’ approaches of Section 4, as well as the Heuristic Algorithm of Section 5. Throughout this section, we assume that the expenditure function is of exponential form, as defined by (10). Each problem instance was generated as follows: • We varied the number of candidate nodes jEj between 20 and 200 (in steps of 10) and the number of new facilities m between 2 and 21. The number of pre-existing facilities jEj was drawn randomly from {1, . . . , 2m}. • We set the number of customer (demand) nodes to jEj þ jEj. Thus jN j ¼ 2ðjEj þ jEjÞ, but only the customer nodes have positive demand. Weights (amount of total demand) for each node were randomly drawn from a uniform distribution on [1, 10]. • We randomly generated the demand elasticity parameter ki for node i from a uniform distribution on [0.1, 1]. • We generated the utility coefficients uij for each customer node i and candidate location j as follows. We A used the form uij ¼ d bj , with distance sensitivity b randomly drawn from [1, 4], the shortest distances dij ij obtained by first randomly generating distances for all links in the graph and then constructing a shortest distance table (the distances for links were generated from a uniform distribution on [5, 100]), and the attractiveness parameter Aj was generated from a uniform distribution on [1, 100]. We note that this form of the utility coefficients corresponds to (2) – which is the ‘‘traditional’’ form for Huff-type models. For each combination of jEj and m values, we generated five random instances. For each problem, we applied the ‘‘exact’’ and the ‘‘a-optimal’’ approaches with a = 0.001. The time limit for all algorithms was set to 60 minutes of CPU time. The resulting optimization models were solved using CPLEX IP solver on a Dell Poweredge 6300 machine running Redhat Linux. Results of the experiments can be found in Table 3 for moderate size problems, and in Table 4 for larger problems. The following measures were recorded for each run (and then averaged for each combination of jEj and m values): the relative errors for the objective function values for the a-optimal and the heuristic solutions (columns 3 and 4), and the CPU time in seconds (columns 5–7 in Table 3 and columns 5–6 in Table 4); the 613 R. Aboolian et al. / European Journal of Operational Research 181 (2007) 598–619 Table 3 Results of the experiment for moderate size problems # Candidate sites (jEj) 20 30 40 40 50 50 100 100 100 120 120 120 160 160 200 200 # Facilities to locate (m) 4 4 4 5 4 5 2 3 4 2 3 4 2 3 2 3 Relative error CPU time (seconds) Heuristic 0.1%-Optimal Heuristic 0.1%-Optimal Exact 0.0106 0.0116 0.0000 0.0029 0.0000 0.0000 0.0012 0.0000 0.0033* 0.0000 0.0000 0.0060* 0.0000 0.0017 0.0000 0.0010* 0.0002 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0000 0.0010* 0.0000 0.0000 0.0010* 0.0000 0.0000 0.0000 0.0010* 0.00 0.01 0.01 0.02 0.02 0.02 0.02 0.03 0.03 0.02 0.03 0.03 0.02 0.03 0.03 0.03 0.93 1.63 1.52 15.40 1.11 1.93 2.57 8.28 12.35 6.92 9.85 15.77 37.59 43.59 67.91 102.38 3.26 32.78 123.45 232.09 218.12 1042.32 5.36 479.35 N/A 9.71 789.02 N/A 29.15 1452 77.42 N/A For entries marked with (*) relative error is evaluated versus the best known upper bound. relative errors were computed as (Z*  Z)/Z, where Z* is the optimal objective function value and Z is the value of the candidate solution. Entries marked ‘‘N/A’’ in the last column of Table 3 mark combinations of parameters for which the ‘‘exact model’’ failed to converge within the specified time limit—in this case the relative errors were measured versus the best available upper bounds; the best available upper bound was computed as the smaller of the following quantities: Za(xa) (i.e., the objective function of the problem P a0 corresponding to its optimal solution xa) and the best intermediate upper bound produced by the exact procedure of Section 4 before the time limit was reached. Note that since for all problems in Table 4 the ‘‘exact model’’ failed to converge, the relative errors are measured versus the best available upper bound. The following observations can be made based on experimental results: 1. The ‘‘a-optimal’’ approach appears to perform extremely well even with a = 0.1%. The computational times are quite reasonable and do not appear to grow overly fast with the number of candidate facility sites. The computation times do increase with the number of facilities to be located, however problems with up to 20 facilities can be solved in reasonable time. The deviations from the optimal solution appear to be negligible—in fact, as shown in Table 3, the optimal solution is found for over 90% of all ‘‘moderate-size’’ instances solved to optimality; the largest observed departure from optimality was 0.02%. 2. The ‘‘exact’’ approach appears to be less attractive, with solutions times much larger than for the ‘‘a-optimal’’ approach. The exact approach appears to be viable only for small-to-moderate size instances of the problem. 3. The Heuristic Algorithm performs extremely well. For the moderate-size instances (Table 3), the average deviation from the optimal value of the objective function was only 0.19%; the maximum deviation recorded was 1.3% for the instances that were solved to optimality; the optimal solution was found for 47% of these instances. Even for larger size problems (Table 4), the maximum relative error versus the upper bound did not exceed 6.2%. Moreover, there is no evidence that the size of the relative error tends to increase with the problem size. The Heuristic Algorithm is quite fast, with the average CPU time of only 0.02 seconds, increasing to 0.03 seconds for the largest instances solved. It appears that both a-optimal and the heuristic approaches are viable for dealing with large-scale instances of the problem. Of course, if the problem can be solved to optimality in a reasonable time using the a-optimal 614 R. Aboolian et al. / European Journal of Operational Research 181 (2007) 598–619 Table 4 Results of the experiment for larger size problems # Candidate sites (jEj) # Facilities to locate (m) Relative error (versus upper bound) CPU time (seconds) Heuristic 0.1%-Optimal Heuristic 0.1%-Optimal 100 100 100 100 120 120 120 120 120 120 120 120 120 120 160 160 160 160 160 160 160 160 160 160 160 160 200 200 200 200 200 200 200 200 200 200 200 200 200 200 5 10 15 20 6 6 8 8 9 9 12 12 15 15 6 6 8 8 9 9 12 12 15 15 18 18 6 6 8 8 9 9 12 12 15 15 18 18 21 21 0.0034 0.0142 0.0030 0.0010 0.0061 0.0010 0.0112 0.0260 0.0004 0.0010 0.0015 0.0015 0.0090 0.0027 0.0079 0.0010 0.0074 0.0010 0.0011 0.0010 0.0245 0.0037 0.0269 0.0010 0.0619 0.0033 0.0010 0.0059 0.0010 0.0033 0.0010 0.0202 0.0302 0.0007 0.0150 0.0038 0.0010 0.0010 0.0328 0.0012 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 0.001 3.56 3.85 4.13 4.86 4.02 4.01 3.90 3.48 4.33 5.13 5.63 5.32 6.25 5.89 6.58 8.41 6.55 6.37 7.55 7.82 8.01 8.95 9.65 9.21 10.06 9.97 15.43 15.26 12.15 16.81 15.39 17.23 17.65 18.15 21.3 19.65 21.54 20.95 23.50 22.79 542.36 875.31 1023.60 953.08 675.38 795.51 738.09 1164.86 1325.84 1369.21 1865.63 1653.28 1135.91 986.26 386.94 628.68 1212.57 996.58 1296.36 1358.60 1624.61 1298.12 2013.56 2154.09 1785.32 1568.37 1119.5 1710.23 1004.69 331.21 1652.09 1065.83 2013.84 1864.66 2216.2 2315.65 1985.32 2865.54 2132.14 2291.07 approach, it should be preferred since it provides a guarantee on the quality of the solution. Where this is not viable due to problem size, the heuristic method can be applied with a reasonable expectation of obtaining a near-optimal solution. We note that the excellent performance of the Heuristic Algorithm is not surprising—it is our experience that this performance is, in fact, to be expected for problems of this type where the sub-modularity-based bound established in Section 5 holds (see, e.g., [5]). On the other hand, we underscore that the results discussed in this section are based on a limited number of randomly generated instances, and thus may not extend to problem instances with structure significantly different from the one used above; in particular we do not know how sensitive the algorithmic performance is to the form of the utility coefficients and other problem parameters (we would expect this to be more of an issue for the Heuristic Algorithm than for the a-optimal procedure, since the latter employs a generic IP solver). R. Aboolian et al. / European Journal of Operational Research 181 (2007) 598–619 615 7. Concluding remarks In this paper, we presented a model for locating competitive service facilities with concave demand. The key advantage of the model is that it captures two key aspects of the underlying problem: the market expansion and market cannibalization effects. Previous modeling attempts in this direction had run into significant computational problems due to non-linearity and complexity of the underlying optimization models. By focusing on the concave demand functions we were able to overcome these difficulties through the development of the TLA procedure which allowed us to obtain a-optimal solutions to large-scale models using standard Integer Programming software. In addition, we developed a Heuristic Algorithm with a tight worst-case error bound, which is somewhat unexpected in view of the non-linearity of the underlying model. Our computational experiments show that this heuristic has excellent performance both, with respect to solution quality, and the computational times. In fact, in view of the performance of the heuristic it may be asked whether the a-optimal procedure is really necessary. In our opinion, the answer is ‘‘yes’’. While the performance of the heuristic in our computational experiments is certainly encouraging, it must be remembered that, in the worst case, the heuristic solution may be up to 37% worse than the optimal one. Given the strategic nature of the facility location problems, and the large capital commitments involved in the locational decisions, it has been our experience that real-world decision makers are unwilling to accept this level of potential deviation from optimality, especially when the uncertainty regarding the solution quality can be resolved with a few hours (or even days) of computer time. The a-optimal procedure has a distinct advantage in this regard since by setting the a level sufficiently low one can guarantee that the potential ‘‘error’’ in the solution is below the ‘‘noise’’ level introduced by other parameters in the model. In future we plan to concentrate on incorporating facility design decisions into our model, i.e., focusing on problems where the decision is not just where to locate the new facilities, but what kind of facilities in terms of their design (and thus their attractiveness level Aj) they should be, where ‘‘design’’ can include many aspects such as facility layout, size, customer service, etc. Appendix Proof of Theorem 3. Let the lth segment defining xA(Æ) have endpoints cl, and cl+1 and the slope bl. To prove the result it suffices to show that  A  xA ð/Þ x ðcl Þ xA ðclþ1 Þ 2 max ; : /2½cl ;clþ1  xð/Þ xðcl Þ xðclþ1 Þ This holds trivially if x(/) is constant on [cl, cl+1]. For the remainder of the proof, assume that x(/) is increasing on this interval. Þxðcl Þ For / 2 [cl, cl+1] let Sð/Þ ¼ xðcl Þ þ ð/  cl Þ xðcðclþ1 be a straight line connecting points (cl, x(cl)) and lþ1 cl Þ Þxðcl Þ be the slope of this line. Since x(/) is assumed to be concave and increas(cl+1, x(cl+1)), and let sl ¼ xðcclþ1 lþ1 cl ing on [cl, cl+1], it follows that sl > 0 and Sð/Þ ¼ xð/Þ for / 2 fcl ; clþ1 g and Sð/Þ < xð/Þ for / 2 ðcl ; clþ1 Þ: ð21Þ Thus, xA ð/Þ xA ð/Þ P Sð/Þ xð/Þ for / 2 ½cl ; clþ1 ; with equality holding for / = cl and / = cl+1. A l Þxðcl Þ Let r ¼ bsll . By definition of the relative error (/), we have ðcl Þ ¼ x ðcxðc . Now, for / 2 [cl,cl+1], lÞ xA ð/Þ ¼ xA ðcl Þ þ bl ð/  cl Þ ¼ ð1 þ ðcl ÞÞxðcl Þ þ rsl ð/  cl Þ ¼ ð1 þ ðcl ÞÞ½xðcl Þ þ sl ð/  cl Þ þ ½r  ð1 þ ðcl ÞÞsl ð/  cl Þ: ð22Þ 616 R. Aboolian et al. / European Journal of Operational Research 181 (2007) 598–619 It follows that xA ð/Þ sl ð/  cl Þ ¼ ð1 þ ðcl ÞÞ þ ½r  ð1 þ ðcl ÞÞ : Sð/Þ xðcl Þ þ sl ð/  cl Þ A A ð/Þ ð/Þ is non-decreasing is r P 1 + (cl) and is decreasing otherwise. In either case, xSð/Þ is We conclude that xSð/Þ maximized at one of the endpoints of [cl, cl+1]. The result now follows from (21) and (22). Detailed algorithm and technical results for the TLA procedure Steps B and C of the informal procedure described in Section 3.1 above are obviously the crucial ones in the procedure described above. They are based on the two results presented below.  and d 2 R such that x(c) 6 d. For a point cT 2 ½c; /,  let bL ðcT Þ ¼ xðcT Þd be the Lemma 3. Consider c 2 ½0; / cT c slope of the line from (c, d) to (cT, x(cT)), and let L(/) = d + bL(cT)(/  c), / 2 [c, cT] be the equation of that line. Then  such that x 0 (cT) = bL(cT) (i.e. the slopes of L(/) and x(/) coincide at cT) if and only if 1. There exists cT 2 ½c; /  /   cÞ þ d:  P x0 ð/Þð ð23Þ xð/Þ  2. If condition (23) holds, then cT is unique, and can be found by solving the following equation for / 2 ½c; /: x0 ð/Þð/  cÞ þ d  xð/Þ ¼ 0: ð24Þ  consider the function T(/) = x 0 (/)(/  c) + d  x(/). Note that T 0 (/) = x00 (/)(/  c). Proof. For / 2 ½c; /,  (in fact, it is strictly Since, by assumption, x(/) is concave, x00 (/) < 0, and thus T(/) is non-increasing on ½c; /  decreasing on ½c; /). Recall that d P x(c). If d = x(c), then T(c) = 0 and thus cT = c provides the unique solution to T(/) = 0.  P c, it follows that T ð/Þ  6 0, which shows that Finally, since T(/) is non-increasing, T(c) = 0, and / condition (23) holds as well.  It follows that T(/) = 0 has a Assume now that d > x(c). Then T(c) > 0, and T(/) is decreasing on ½c; /. solution if and only if condition (23) holds. Let cT be the root of T(/) = 0. Then bL(cT) = x(cT) by definition of bL, thus completing the proof. h Note that if the solution to (24) exists, then it can be found by a standard line search procedure—e.g., bisection. The convergence of the line search procedure is guaranteed by the monotonicity of the function T(Æ) defined in the proof of Lemma 3 above.  For / 2 ½c; /,  let L(/) = x(c) + x 0 (c)(/  c) be the line tangent to x(/) at point Lemma 4. Consider c 2 ½0; /.  such that c. Then, for a > 0 there exist cE 2 ½c; / LðcE Þ  xðcE Þð1 þ aÞ ¼ 0 if and only if  P xð/Þð1  Lð/Þ þ aÞ: ð25Þ ð26Þ Moreover, if a solution to (25) exists, then it is unique.  This function is concave (by the concavity of x(/)), and Proof. Let T(/) = x(/)(1 + a)  L(/) for / 2 ½c; /.  If (26) holds, then T ð/Þ  6 0 and such a root exists. The positive at c. Thus, it has at most one root on ½c; /. converse also holds. h We note that if (26) holds, the unique root of (25) can be found by a standard line search technique (e.g., bisection)—the convergence of the latter is guaranteed by the uniqueness of the root. We can present the details of TLA Procedure for constructing the approximating function xa. The TLA procedure Step 1: Set l = 1, c = c1 = cT = 0, b1 = x 0 (0).  L = l and STOP. Step 2: If condition (26) fails, set clþ1 ¼ /, R. Aboolian et al. / European Journal of Operational Research 181 (2007) 598–619 617 Else, use the bisection method to find the unique root of (25) and set cl+1 equal to that root. If  set L = l and STOP. Else, set c = cl+1 and proceed to Step 3a. clþ1 ¼ /, Step 3a: Set l = l + 1.  If If condition (23) holds, proceed to Step 3b. Else, set L = l and cL ¼ /.  xðcl Þð1 þ aÞ 6 xð/Þ; ð27Þ  lÞ then set bl ¼ xð/Þð1þaÞxðc and STOP.  /cl Else, set bl = 0 and STOP. Step 3b: Use the bisection method to find the unique root of (24), set cT equal to that root and bl = x 0 (cT).  set cl+1 = cT, L = l and STOP. Step 4: If cT ¼ /, Else, repeat Step 2. Comments Note that cT is used as a temporary variable in the procedure above—it is the point where segment l is tangent to x(/). We note that from the two preceding results it follows that if a stopping point has not been reached (i.e., if the next iteration is initiated in Step 4), the relative error for the current segment will be a at its right endpoint. Thus, the relative error at the left endpoint is 0 for the first segment and a for every subsequent segment. Since we know from Theorem 3 that the maximum relative error for each segment is reached at one of the endpoints, the overall relative error of a is assured as long as the relative error at the right endpoint of the last segment does not exceed a. This is easy to check by analyzing the possible stopping conditions of the procedure above. The procedure has several stopping points. In Step 2 it stops if either condition (26) is not satisfied (which—  and thus the current segment can be extended by Lemma 4—implies that the relative error is less than a at /,  In either case,  or if the right endpoint of the current segment (computed by Lemma 4) coincides with /. to /), the relative error for the last segment cannot exceed a. In Step 3a, the stopping point is reached if condition (23) fails, indicating that the tangent line to x(/) can In this case, it is natural to simply extend the current segment to the point not be constructed on ½cl ; /.  xð/ÞÞ,  thus ensuring that the relative error for the final segment is 0. However, this extension may lead ð/; to the negative slope for the last segment if condition (27) holds—which would violate the requirement that xa(/) must be non-decreasing. Thus, if condition (27) holds, we simply set the slope of the last segment to 0.  6 xa ð/Þ  ¼ xa ðcl Þ ¼ ð1 þ aÞxðcl Þ 6 ð1 þ aÞxð/Þ  (where we have used the fact that Note that in this case xð/Þ x(/) is non-decreasing and that the relative error is a at the previous endpoint). This shows that the relative error is no more than a for the final segment in this case as well. Finally, if the stopping point is reached in Step 4, then the tangent point for the last segment has coincided  and thus the relative error at /  is 0. with /, The discussion above shows that the relative error for the piecewise linear function xa(/) constructed by the TLA procedure does not exceed a (in fact, the relative error is equal to a at the right endpoint of every segment, with the possible exception of the last segment). Moreover, by construction, the slope of every segment exceeds that of the following segment, since the slope of each segment equals to the slope of x(/) for some / belonging to the segment, and x(/) is concave. This ensures that xa(/) is concave and non-decreasing. Proof of Theorem 4. The properties of xa(/) constructed by the TLA procedure have already been established in the discussion above. It remains to show that the procedure converges in finitely many steps. Consider segment l in xa(/) for l 2 {1, . . . , L  1}. It consists of two sub-segments: [cl, cT] and [cT, cl+1]. Note that xa(cT) = x(cT) and xa(cl+1) = (1 + a)x(cl+1). Thus, clþ1  cl > clþ1  cT ¼ ð1 þ aÞxðclþ1 Þ  xðcT Þ axðclþ1 Þ P 0 ; x0 ðcT Þ x ðcT Þ where the last inequality follows since x(/) is increasing. Since x(/) is increasing and concave, the quantity on the right-hand side above is increasing with l. Thus, the lengths of the segments in xa(/) are bounded below by an increasing sequence (in particular, the length of the sub-segment [cT, cl+1] is strictly increasing with l). Let  > 0 be the length of segment 1. Since for this 618 R. Aboolian et al. / European Journal of Operational Research 181 (2007) 598–619 segment, cT = c1 = 0,  is the length of the [cT, c2] sub-segment, and thus the length of every other segment is  L ¼ 1 þ d/ e segments, and the TLA procedure stops after at bounded below by . Thus, xa(/) has at most b most b L major iterations, which concludes the proof. P Proof of Theorem 7. To prove that Z O ðSÞ ¼ i2N wi V i ðU i ðSÞÞ the set function in the original problem is nondecreasing we note the reader to the fact that Vi(Ui(S)) is a non-decreasing function in Ui(S) and since from definition Ui(S) is a non-decreasing P function in Ui(S), then Vi(Ui(S)) is a non-decreasing function in S. Now since wi P 0 for i 2 N, therefore i2N wi V i ðU i ðSÞÞ is also a non-decreasing function in S. P P Define the following terms: ti ¼ j2T uij , si ¼ j2S uij , where S  T  E. To prove that ZO(S) is submodular we have to show that ZO(T [ {k})  ZO(T) 6 ZO(S [ {k})  ZO(S) where S  T  E and k 62 T. Since Vi(U) is concave, for 0 6 u1 6 u2 and r P 0 we have V i ðu1 þ rÞ  V i ðu1 Þ P V i ðu2 þ rÞ  V i ðu2 Þ for i 2 N : ð28Þ Since uik P 0 and si 6 ti for i 2 N, from (28) we obtain V i ðsi þ uik Þ  V i ðsi Þ P V i ðti þ uik Þ  V i ðti Þ for i 2 N : ð29Þ From (29) it follows that X X X X wi V i ðti Þ: wi V i ðti þ uik Þ  wi V i ðsi Þ P wi V i ðsi þ uik Þ  i2N i2N i2N ð30Þ i2N Thus, ZO(S) is sub-modular, which concludes the proof. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] R. Aboolian, Competitive facility location and design problem, Ph.D. thesis, University of Toronto, Toronto, Canada, 2002. D. Achabal, W.L. Gorr, V. Mahajan, MULTILOCC, a multiple store location decision model, Journal of Retailing 58 (1982) 5–25. M.S. Bazaraa, H.D. Sherali, C.M. Shetty, Nonlinear Programming: Theory and Algorithms, John Wiley and Sons Inc., 1993. R.J. Beaumont, Spatial interaction models and the location-allocation problem, Journal of Regional Science 20 (1) (1980). O. Berman, D. Krass, Locating multiple competitive facilities: Spatial interaction models with variable expenditures, Annals of Operations Research 111 (2002) 197–225. O. Berman, D. Krass, Flow intercepting spatial interaction model: A new approach to optimal location of competitive facilities, Location Science 6 (1998) 41–65. R.D. Buzzell, Mathematical Models and Marketing Management, Boston, Division of Research, Graduate School of Business Administration, Harvard University, 1964. K.M. Bretthauer, B. Shetty, The nonlinear Knapsack problem—algorithms and applications, European Journal of Operational Research 138 (2002) 459–472. A. Caprara, D. Pisinger, P. Toth, Exact Solution of the quadratic knapsack problem, INFORMS Journal on Computing 11 (1999) 125–1371. E.W. Cheney, Introduction to Approximation Theory, McGraw-Hill, New York, 1966. M.L. Cornuejols, G.L. Nemhauser, L.A. Wolsey, The uncapacitated facility location problem, in: R.L. Francis, P. Mirchandani (Eds.), Discrete Location Theory, Wiley Interscience, New York, 1990. T. Drezner, Competitive facility location in the plane, in: Z. Drezner (Ed.), Facility Location, Springer-Verlag, 1995, pp. 291–298. T. Drezner, Optimal continuous location of a retail facility, facility attractiveness, and market share: An interactive model, Journal of Retailing 70 (1) (1994) 49–64. J.P. Dussault, J.A. Ferland, B. Lemaire, Convex quadratic programming with one constraint and bounded variables, Mathematical Programming 36 (1986) 90–104. A. Ghosh, S.L. McLafferty, C.S. Craig, Multifacility retail networks, in: Z. Drezner (Ed.), Facility Location, Springer-Verlag, 1995, pp. 301–330. D.S. Hochbaum, A nonlinear knapsack problem, Operations Research Letters 17 (1995) 103–110. H. Hotelling, Stability in competition, Economic Journal 39 (1929) 41–57. D.L. Huff, Defining and estimating a trade area, Journal of Marketing 28 (1964) 34–38. T.D. Klastorin, On a discrete nonlinear nonseparable knapsack problem, Operations Research Letters 9 (1990) 233–237. J.D.C. Little, L.M. Lodish, A media planning calculus, Operations Research 18 (1969) B25–B40. L.M. Lodish, D.B. Montgomery, F.E. Webster, A dynamic sales call policy model, working paper 329-68, Sloan School of Management, M.I.T. Cambridge, 1968. R.R. Meyer, Two-segment separable programming, Management Science 25 (4) (1979) 385–395. R. Aboolian et al. / European Journal of Operational Research 181 (2007) 598–619 619 [23] R.R. Meyer, Computational aspects of two-segment separable programming, Mathematical Programming 26 (1983) 21–39. [24] M. Nakanishi, L.G. Cooper, Parameter estimates for multiplicative competitive interaction models-least square approach, Journal of Marketing Research 11 (1974) 303–311. [25] Ph. Naret, M. Weverbergh, On the predictive power of market share attraction models, Journal of Marketing Research 18 (1981) 146– 153. [26] G. Nurenberger, Approximation by Spline Functions, Springer-Verlag, Berlin, 1989. [27] P.M. Pardalos, H. Ye, C. Han, Algorithms for the solution of quadratic Knapsack problem, Linear Algebra and its Applications 152 (1991) 69–91. [28] F. Plastria, The generalized big square small square method for planar single facility Location, European Journal of Operations Research 62 (1992) 163–174.