Discretization based heuristics for the capacitated multi-facility Weber problem with convex polyhedral barriers

Article History: Received 04 September 2016 Accepted 13 September 2017 Available 23 October 2017 The Capacitated Multi-facility Weber Problem (CMWP) tries to determine the location of I capacitated facilities in the plane and to satisfy demand of J customers so as to minimize the total transportation cost. The CMWP assumes that the facilities can be located anywhere on the plane and customers are directly connected to them. This study considers an extension of the CMWP where there exist convex polyhedral barriers blocking passage and locating facilities inside. As a result, the distances between facilities and customers have to be measured by taking into account the polyhedral barriers. The CMWP with convex polyhedral barriers (CMWP-B) is a non-convex problem that is difficult to solve. We propose specially tailored discretization based heuristic procedures. Since CMWP-B is novel to the literature, a new set of test problems is randomly generated. Then, the performance of the suggested methods are tested on the test instances. Our results imply that the suggested heuristics yield quite accurate and efficient solutions for the CMWP-B.


Introduction
The Capacitated Multi-facility Weber Problem (CMWP) tries to locate I facilities with capacity restrictions in the plane and to meet the demand of J customers while minimizing the total cost of transportation.The objective function of the CMWP is shown to be neither convex nor concave [1].When capacity restrictions of the facilities are removed, the CMWP becomes the so called Multi-facility Weber Problem (MWP).Both of the MWP and CMWP are NP-hard as shown by Sherali and Nordai [2], and Meggido and Supowit [3], respectively.Moreover, they both generalize the Weber Problem (WP) that aims to find the optimal location of a single facility.The objective function of the WP is convex, and thus, easy to solve.However, this is not true for its generalizations of the MWP and CMWP as additional allocation decisions have to be made for them.One can partition the customer set into non-intersecting subsets each of which is served by an uncapacitated facility for the MWP.This indicates that each customer is served from exactly one facility at the optimal solution of the MWP.On the other hand, customer demands may need to be met from multiple facilities for the CMWP when capacity restrictions are imposed.Nevertheless, all of the WP, MWP and CMWP assume that each facility can be freely located in the plane.Such an assumption can be misleading in practice since there may exist physical areas that obstruct to travel and to open facilities inside.For example, lakes, mountains, glaciers and forests are natural barrier areas which can prohibit both travelling and facility location inside.Indeed, traveling (or passage) is frequently blocked by existing work-shop (office room) areas in a manufacturing plant (in an office building).This work focuses on an extension of the CMWP where the existence of convex polyhedral barriers *Corresponding Author are considered as obstacles to pass inside and to locate facilities.Namely, the CMWP with convex polyhedral barriers (CMWP-B) is addressed.We suggest several discretization based heuristic approaches to efficiently solve the CMWP-B.Mostly, the transportation cost between a facility and customer is a function of the distance and the amount of shipment between them.In particular, the distance is usually modelled with the Euclidean, squared Euclidean, rectilinear and ℓ r (with 1 ≤ r < ∞) norms in location-allocation type problems such as the CMWP.Alpaydın et al. [4] and Brimberg et al. [5] provide excellent surveys that include techniques and functions to model distances.In this work, the distances between facilities and customers are measured by the so called "barrier distance" which uses the Euclidean distance, i.e. ℓ 2 , as its underlying distance measure.Observe that the Euclidean distance may not be feasible when polyhedral barriers exist between facilities and customers.Therefore, "barrier distance" is determined as the shortest path, which does not pass through the barrier regions where the edge costs are calculated with the Euclidean distance, between a facility and a customer.Clearly, the "barrier distances" require extra efforts to calculate the distances which makes the CMWP-B a more realistic and difficult problem to solve than the unrestricted CMWP.Nonetheless, the "barrier distance" still preserves the metric properties [6].This work has the following contributions.For all we know, the CMWP-B is novel to the literature and we suggest a mathematical programming formulation for it.Next, we offer three discretization based (DB) heuristics which reduce the continuous location space into a discrete space using a discretization strategy.The first DB heuristic employs the solution of an approximating mixed-integer linear programming (MILP) formulation suggested for the CMWP-B.The second DB heuristic applies a Lagrangean Relaxation (LR) scheme on the suggested MILP formulation.The third DB heuristic performs a Tabu Search (TS) using neighborhoods defined over the discretized location space.All upper bounds are calculated with alternate location-allocation type heuristics for the DB heuristics.We generate new test instances for the CMWP-B where they inherit the data of the standard CMWP instances.For that purpose, convex polyhedral barriers are randomly constructed such that the feasibility of the CMWP-B instances is maintained.Lastly, we perform extensive computational experiments on randomly generated test instances derived from the standard CMWP instances.
The remainder of this study is organized as follows.In Section 2 a brief review of the relevant literature is presented.This is followed by the mathematical programming formulation of the CMWP-B in Section 3. Section 4 presents details of the suggested DB heuristics.Section 5 is where we discuss our computational findings.Lastly, the paper is concluded and future research directions are given in Section 6.

Literature Review
In this section, we present a short review of the relevant literature.Both the MWP and CMWP have attracted the attention of researchers starting with the seminal works by Cooper [1,7].There exist several exact solution procedures suggested for the MWP ( [8][9][10][11]) as well as heuristic approaches for the MWP ( [7,[12][13][14][15][16][17][18]).Furthermore, we can cite the works by Cooper [1], Sherali et al. [19] and Akyüz et al. [20] as examples of exact solution methods developed for the CMWP.On the other hand, we can mention Cooper's [7] Alternate Location-Allocation (ALA) type heuristics [21,22], Discrete Approximation (DA) heuristics [23,24] and metaheuristics [25] developed for the CMWP.Akyüz et al. [26,27] also offer heuristic procedures for a multi-commodity extension of the CMWP.A survey on the location-allocation type problems can be found in the work by Brimberg et al. [28].As an extension of the single facility WP, the WP-B arises when there exist barriers which prohibit to travel inside and to locate a facility.Katz and Cooper [29] introduce the WP-B having a single circular barrier region.The WP-B with rectilinear distances is considered in the works by Larson and Sadiq [30] and Batta et al. [31].Aneja and Parlar [32] design an algorithm for the WP-B as well as for a variant of the WP-B where travelling is permitted within the barriers.Butt and Cavalier [33] address the WP-B where barriers are convex polygons and propose an algorithm which yields local optimum solution.Klamroth [34] derives a reduction result for the WP-B having convex polyhedral barriers.Then, an exact and a heuristic procedure is developed for the WP-B.McGarvey and Cavalier [35] develop a branch and bound approach that partitions the continuous location space to optimally solve the WP-B in the presence of polyhedral barriers.Bischoff and Klamroth [36] deal with the WP-B where the barriers are convex polyhedral sets.Their method benefits from the reduction result by Klamroth [34] and decomposes the WP-Bs into multiple subproblems (i.e.WPs) each can be solved over a convex and bounded region.The authors [36] apply a genetic algorithm to decrease the number of subproblems to solve.
The MWP with convex polyhedral barriers (MWP-B) is first considered by Krau [10] which applies a column generation approach based on partitioning of the customer set.Bischoff et al. [37] suggest two ALA type heuristics that have similar allocation phases for the MWP-B.Their first ALA type heuristic alternately solves multiple WP-Bs and set partitioning problems resulting in an inefficient algorithm.The second ALA type heuristic suggested by Bischoff et al. [37] yields better results.In the location phase, the reduction result of [34] is extended to the multifacility case.Unlike their first ALA type heuristic, this results in solving multiple WPs over convex restricted regions and increase the efficiency of the location phase for the second ALA type heuristic.In what follows, we give a formal definition of the CMWP-B.

Capacitated Multi-Facility Weber Problem with Convex Polyhedral Barriers
Let I, J and P denote the number of facilities, the number of customers and the number of polyhedral barriers, respectively.The coordinates of customer j is shown by a j = (a j1 , a j2 ) T and its demand is denoted by q j .The parameter s i stands for the capacity of facility i.There are two decisions to be made for the CMWP-B: location decisions and allocation decisions.Then, the unknown location of facility i is represented as x i = (x i1 , x i2 ) T .The unknown amount of flow between facility i and customer j is denoted by f ij .c ij is the cost of transportation per unit flow between customer j and facility i per unit distance.A polyhedral barrier p is indicated with set B p and the union of barriers are stated with the set B, i.e.B = P p=1 B p .Now, the feasible region to locate facilities is given with the set X that can be defined as X = E 2 \ B where E 2 is the two dimensional Euclidean space.The Euclidean norm, i.e.
, is used as the underlying distance measure to calculate the "barrier distance" between facility i and customer j.Note that barrier distance is represented with d B (x i , a j ) between facility i and customer j.The details about its calculation will be discussed later on.Now a mathematical programming formulation of the CMWP-B is given as follows. CMWP-B: The objective function ( 1) is to minimize the sum of total transportation cost between facilities and customers.Notice that constraints (2)-( 4) are the constraints of the well-known Transportation Problem (TP).Constraints (2) state that the total flow sent from a facility i is equal to its capacity s i .Constraints (3) imply that the demand q j of customer j is exactly met.Constraints (4) ensure the nonnegativity of flows.Constraints (5) guarantee that facilities are placed within the feasible region X .It is possible to further restrict the set of feasible region X .Wendell and Hurter [38] show that an optimal solution of the WP can be found within the convex hull of customers.The result by Wendell and Hurter [38] can be directly generalized to the WP-B as well as to the CMWP-B in the case where barrier regions stay within the convex hull of customer locations.In this case, it is enough to consider the region within convex hull of customers that is outside of the polyhedral barriers.However, unlike the CMWP, the resulting feasible region to locate facilities is non-convex for the CMWP-B.Moreover, distances between facilities and customers are measured with barrier distances d B (x i , a j ) for the CMWP-B.Observe that, the CMWP employs the Euclidean distance between each facility and customer.However, barrier distance is calculated as the shortest path distance between any two points when there exist barriers.The CMWP-B can be reduced to the CMWP when all barriers are removed.The optimal objective value of the CMWP is a lower bound on the optimal objective value of the CMWP-B since barrier distance between any two points is larger than or equal to their Euclidean distances.That is to say, the CMWP-B is more difficult to solve than the CMWP which is NP-hard.This motivates the implementation of efficient heuristic methods for the CMWP-B.The discussion above only considers the barrier distances among the nodes existing in the visibility graph.For an arbitrary point x within the set of feasible locations, i.e. x ∈ X , the calculation of barrier distances is slightly different.One may consider re-constructing the visibility graph from scratch such that x is also added into the node set.Fortunately, it is sufficient to determine only the set of nodes visible from x, say V x , in the visibility graph.Then, the barrier distances from the point x to any customer location a is measured using the formula: where h stands for the customer locations and/or extreme points of barriers that are visible from the point x.SP D(h, a) represents the shortest path distance between points h ∈ V x and a in the visibility graph.The barrier distance function d B (x, a) defines a metric on X satisfying the following properties of a metric: positivity, definiteness, symmetry and triangle inequality.Klamroth [6] is an excellent reference to resort for more details on the single facility location problems with barriers and their properties.In the left-hand side of Figure 1, an example with four customers and single tetragon barrier is given.In the right-hand side of Figure 1, the corresponding visibility graph is illustrated.

Discretization Based (DB) Heuristics
The CMWP-B reduces to solving I WP-B's when the amount of shipments are known.This results in a two dimensional minimum location problem in continuous space.When all barriers are removed, the resulting problem is a classical WP and using the results by Wendell and Hurter [38] and Hansen et al. [41], an optimal solution can be found within the convex hull of customer locations.Clearly, the single WP is easy to solve because it is a convex programming problem.That is to say, it can be solved by the Weiszfeld's [42] algorithm or one of its generalizations ( [43,44]).
On the other hand, the objective function of the CMWP-B is non-convex and the problem itself is not easy to solve on continuous space.This sparks a discretization strategy to solve the CMWP-B using a discretized location space.It is conceivable that the CMWP-B can be transformed into a MILP problem formulation when facility locations are chosen from a set of candidate locations.Indeed, the resulting MILP problem formulation may yield the optimal solution of the CMWP-B when the set of candidate locations to place facilities includes the optimal facility locations.Obviously, it is not possible to determine the optimal facility locations in advance when constructing the set of candidate facility locations.However, solving a MILP problem considering a set of candidate facility locations produces an approximate solution for the CMWP-B.For that purpose, a systematic way of choosing candidate facility locations is of high importance and may provide good solutions for the CMWP-B.In fact, such a strategy is formerly offered by Hansen et al. [14] and Aras et al. [23] for the MWP and CMWP, respectively.They use customer locations as the set of candidate facility locations and solve approximating MILPs.Hansen et al. [14] solve an approximating p-median problem and obtain highly accurate solutions for the MWP.Similar results also hold for the CMWP from the work by Aras et al. [23].The discretization strategy that we pursue in this study is analogous to the previous ones [14,23].In what follows, we first present the approximating MILP problem formulation and the first DB (DB-I) heuristic.This part also elaborates how we tackle with the difficulties imposed by the barrier distances for the single WP-B.Then, the efficiency of DB-I heuristic is improved by a LR scheme.This heuristic is denominated as DB-II.Lastly, a Tabu Search (TS) algorithm is employed to further improve the efficiency of the DA-I heuristic.This heuristic is called as DB-III heuristic.

Approximating MILP Problem Formulation: DB-I Heuristic
Let k = 1, . . ., K denote candidate facility locations with known coordinates given as v k = (v k1 , v k2 ) T .The decision variables w ijk represent the amount of flow between facility i located at candidate point k and customer j.Binary variables u ik take a value of 1 if and only if facility i is opened at candidate point k and 0 otherwise.c ijk is the corresponding transportation cost for flow w ijk .Specifically, it is calculated as ).Now, an approximating MILP problem formulation of the CMWP-B can be stated as follows. DA: s.t.
Here, constraints ( 8), ( 9) and ( 11) are analogous to the TP constraints (2)-( 4) for the approximating MILP problem formulation of the CMWP-B.Constraints (10) ensure that every facility is opened at exactly one candidate facility location.
Before giving the details of the DB-I heuristic, we first introduce an improvement heuristic that is employed within the DB heuristics in the following.This improvement heuristic is an adaptation of the famous Alternate Location-Allocation (ALA) heuristic that is offered for the MWP by Cooper [7].It consists of alternately solving two subproblems: allocation and location subproblems.These subproblems respectively arise when the facility locations and allocations between facility and customers are fixed.ALA heuristic repeats alternating until the objective function value does not change from one iteration to another.This also indicates that facility locations and allocation values do not significantly change and the objective value becomes stable.It can be observed that initialized from given facility locations, the CMWP-B reduces to solving a TP, which is the "allocation" subproblem, in order to determine allocations, i.e. flow between facilities and customers.TP can be solved straightforwardly using linear programming solvers.On the other hand, with a given feasible shipment plan, the CMWP-B can be decomposed into I single facility WP-Bs so called the "location" subproblems.Solving the resulting location subproblems, i.e. the WP-Bs, is not trivial and our approach for solving the location subproblems is summarized in the following.Recall that the WP-B is a nonconvex problem and several solution approaches can be used to solve ( [32, 33, 35]).Instead of an exact solution procedure, which can be prohibitive to use within an efficient heuristic approach, we prefer to make use of the Weiszfeld's [42] procedure that solves the unrestricted WP optimally.Weiszfeld's algorithm employs the gradient direction to update facility locations from one iteration to the other and eventually converges to the optimal facility locations.This procedure works well on a convex problem like the WP.However, there is no guarantee of optimality and it cannot be directly used for the non-convex WP-B.Nevertheless, observe that, when Weiszfeld's algorithm ends up with a solution x * that is within the feasible location space outside of the polyhedral barriers, i.e. x * ∈ X , then it is also optimal for the WP-B.Otherwise, when Weiszfeld's algorithm terminates with a solution x * that is within a polyhedral barrier p, then, the facility location can be approximated by choosing the closest point x on the border of the barrier p.In other words, our solution approach finds a solution for the location subproblems that yields approximate (or hopefully optimal) facility locations.This procedure is illustrated with Figure 2.This improvement heuristic solves the location and allocation subproblems as described until the objective value remains the same from one iteration to another.This heuristic is named as ALA with barriers (ALAB) heuristic.ALAB heuristic is frequently resorted within our DB heuristics.
The DB-I heuristic works as follows.First, the DA formulation is solved using the set of candidate facility locations consisting of customer locations.The DA formulation yields the facility locations which are used to initialize the improvement heuristic.Second, the ALAB heuristic is used to enhance the solution initialized with the facility locations obtained in the first phase.Then, the best solution found is reported as the outcome of the DB-I heuristic.Note that applying the ALAB heuristic in the second phase of the DB-I heuristic usually improves the initial solution obtained from the DA formulation.

A Lagrangean Relaxation (LR) Scheme: DB-II Heuristic
The DA formulation can be intractable for large problems and thus computationally very expensive to solve exactly.Therefore, it may be better to use an approximate solution approach.To this end, a LR scheme and subgradient optimization is employed to compute good feasible solutions for the DA.The demand constraints (9) are relaxed with Lagrangean multipliers µ j to obtain the Lagrangean subproblem of DA, namely the LDA. LDA(µ): s.t. ( 8), ( 10), ( 11), ( 12), Observe that constraints (14) are introduced as basic upper bounds on the flow quantities to the LDA formulation.Clearly, these constraints are redundant for the DA formulation.However, they may significantly improve the optimal objective value, say Z * LDA (µ), of the Lagrangean subproblem LDA.The last term of the objective function ( 13) is constant and LDA(µ) can be further decomposed over the facilities.That is to say, solving the following subproblems for each facility i = 1, . . ., I is equivalent to solving the Lagrangean subproblem LDA(µ) LDA i (µ): s.t.
w ijk ≤ min{s i , q j } j = 1, . . ., J; k = 1, . . ., K, (18) where the unit costs c ijk are determined as c ijk = (c ijk − µ j ) with a given Lagrange multiplier vector µ.The resulting subproblem LDA i (µ) can be solved by an inspection procedure for each candidate facility location k.Note that, when facility i is placed at a candidate location k, then LDA i (µ) further reduces to the following problem LDA ik (µ): LDA ik (µ) is a continuous bounded knapsack problem that can be optimally solved in polynomial time [45].This approach requires sorting of the cost coefficients for each candidate location k.Hence, we need to run a sorting procedure K times.The least cost candidate point k * is chosen as the optimal location of facility i, and thus for the subproblem LDA i (µ).Now, the optimal value of LDA i (µ) is calculated as It is conceivable that the values of binary variables are set as u ik * = 1 for the least cost candidate location k * and u ik = 0 for k = 1, . . ., K and k = k * .Once all subproblems LDA i (µ) are solved, for a given Lagrange multiplier vector µ, the optimal value of the Lagrangean subproblem LDA(µ) is determined as Z * LDA (µ) = I i=1 Z * LDA i (µ) + J j=1 µ j q j .Notice that Z * LDA (µ) constitutes a lower bound on the optimal value of the DA for any Lagrange multiplier vector µ.The best lower bound can be obtained by solving the Lagrangean dual problem max µ {Z * LDA (µ)} that is accomplished using subgradient algorithm by Held et al. [46].For the sake of brevity, we do not give details of the subgradient algorithm.The subgradient algorithm calculates upper bounds during its run, and hence, feasible solutions for the CMWP-B.For that purpose, the ALAB heuristic is used.The values of binary variables u ik are used to determine the locations of each facility for each solution of the resulting Lagrangean subproblem LDA(µ) with given multiplier values µ.Lastly, the DB-II heuristic reports the best feasible solution obtained from the subgradient algorithm.

A Tabu Search Algorithm: DB-III Heuristic
Inspired with the promising results obtained by using the customer locations as the set of candidate facility location, the DB-III heuristic consists of applying a TS algorithm on the set of candidate facility locations.A simple neighborhood search structure can be defined by exchanging the location of a facility i ′ at a candidate location k ′ with another one in order to move from one feasible solution to another one.This yields a feasible solution for the CMWP-B.Further, feasible solutions can be improved by applying a ALAB heuristic.
The suggested TS algorithm basically exchanges the location of a facility with another candidate location at each iteration.The selected candidate location is declared as tabu for that facility and its status can not be revoked during the tabu tenure.A greedy strategy of selecting the closest candidate location is applied to move from the current feasible solution to a neighbor solution.
To increase the intensity of the neighbor search the closest N neighbor feasible solutions can be checked and the best feasible solution is picked as the new feasible solution.The TS algorithm searches over the candidate location set and reports the best solution found.For further details on TS we refer to the work by Glover and Laguna [47].
The TS algorithm uses a tabu list T that keeps record of candidate facility locations which are declared as tabu for a duration (tabu tenure) say α iterations for each facility i.Each facility i is associated with a set of candidate locations which are declared as tabu.A tabu declared candidate location k for a facility i can not be in the current solution until the tabu tenure record, represented as T (i, k), decreases to zero.Let t denote the current iteration number and θ stand for the maximum number of iterations completed by the TS, namely, the iteration limit of the TS.At each iteration, a facility i * is chosen and a non-tabu candidate location k * with T (i * , k * ) = 0 is declared as tabu for i * .The facility i * is determined with respect to their facility index in the order from smallest facility index to largest facility index.In other words, if facility i is selected at iteration t then facility i+1 is the next facility to be selected at iteration t + 1.Note that, when i * = I for the selected facility, the facility to be declared as tabu is the one with index number (1) for facility i * and Z best T abu is updated.In addition, the current tabu solution is updated so that facility i * is located on the candidate location k ′ (1) .Similarly, such an approach is followed until N th closest candidate location whenever an improvement is achieved.If there is no improvement for the best upper bound Z best T abu after the exchange of N th closest candidate location k ′ (N ) , then facility i * is located at the candidate location which gives the lowest objective value among these N neighbor candidate locations.Selection of the candidate location k * of facility i * to move to another neighbor solution requires further attention.The tabu status of candidate locations terminates after a while and it is likely that the same candidate location is selected repetitively.To avoid such a case, a tabu frequency list (FL) keeps record of the total number of times a candidate point k is declared tabu for a facility i within the TS algorithm.F L(i, k) stands for the FL, and initially, F L(i, k) = 0 for all i = 1, . . ., I; k = 1, . . ., K. When a candidate location k for facility i * is tested within the neighborhood search described at an iteration, F L(i * , k) is increased by one.The candidate location k which has the highest F L(i * , k) value is excluded from the neighborhood search in the next iteration.Then the set of neighbor candidate locations does not contain a tabu declared candidate location or the most frequent candidate location for i * .Therefore, this strategy favors diversification of the solutions during the TS.The upper bound Z t U B obtained at iteration t is associated with a solution vector (k 1 , k 2 , . . ., k I ) which represents the candidate locations of each facility in the solution.Here, k 1 Algorithm 1: Tabu Search Algorithm Step 1. (Initialization): Find initial upper bound Z 0 U B and its associated solution vector of candidate locations (k 1 , k 2 , . . ., k I ) 0 .Set Z best T abu = Z 0 U B and iteration counter t = 1.Set F L(i, k) = 0 for i = 1, . . ., I; k = 1, . . ., K. Set T (i, k i ) = α and F L(i, k i ) = 1 for each facility i = 1, . . ., I.
Step 2. For facility i = 1, . . ., I, (i) determine the neighborhood set K i as 1) and find its associated objective value as Step 3. If Z best T abu does not improve, find k * such that k * = argmin Step 4. Set t = t + 1 and decrease each tabu tenure value T (i, k i ) > 0 by one.If t = θ or Z best T abu does not improve for 30 consecutive iterations STOP and report Z best T abu , otherwise go to Step 2.
is the candidate location of the first facility, k 2 is the candidate location of the second facility and so forth, in the solution at iteration t.A formal outline of the suggested TS algorithm is given in Algorithm 1.

Computational Experiments
In this section, first the test bed used in this work is given.Second, the results obtained with our DB heuristics are presented for the CMWP-B.A Dell Precision T5810 workstation with Intel(R) Xeon(R) E5-1650v3 processor of 3.50 GHz and 64 GB RAM operating within Microsoft Windows 7 Pro 64-bit environment is employed as our computing platform.The callable library of Gurobi 5.6.3 with default settings is used to solve MILP and LP formulations presented and all codes are written in C ++ programming language.

Test Bed
We have performed our computational experiments on randomly generated test instances that are produced using standard CMWP test instances from the literature.Standard test instances are taken from the works by Sherali et al. [19] and Boyacı [48].Recall that the CMWP instances do not have barriers forbidding location and travel.However, the CMWP-B can inherit facility capacities, customers' demand and location information from the CMWP instances.We have generated convex polyhedral barriers for each CMWP instance by considering the following two issues.First, we ensure that barrier regions do not contain any customer location inside.Second, interior areas of barriers should not be coinciding with others.Then, these barriers are integrated with the data of the CMWP instances resulting in CMWP-B test instances.
The optimal objective value of the CMWP constitutes a lower bound on the CMWP-B.We have employed the best known objective values of standard CMWP test instances as benchmark values to make a comparison among the performance of the suggested DB heuristics.These benchmark values representing the best known objective values of the CMWP instances are reported in Table 1.The first column indicates the names of the instances.Original instance numbers are used as stated in the work by Sherali et al. [19].We add a prefix "S" before the corresponding instance number.For example, the CMWP instance number 10 in the work by Sherali et al. [19] is shown as S10.11 CMWP test instances from Sherali et al. [19] are considered.Analogously, a total of 40 CMWP test instances from Boyacı [48], which have nonunique cost values, are represented with a prefix "B" followed by an instance number.The second column stands for the size of the instance by giving the number of facilities and the number of customers in parenthesis, respectively.The instances which have the same number of facilities and customers are distinguished by adding a suffix letter starting from "a" to "e" in accordance with the original denomination presented in the reference works.The last column includes the best known objective values that are taken from the reference works ( [19,27,48].As a remark, these benchmark values are not necessarily optimal for the CMWP as the latter two studies present heuristic outcomes.
To construct barriers for each standard CMWP test instance given, we have randomly generated different number of extreme points and convex polygons.P polygons are chosen from the set P ∈ {1, 3, 5, 10} for the CMWP-B test instances.When P = 1, that is a single polygon, the total number of extreme points, denoted as B, of the polygon is determined within the set H = {2, 3, 4, 5, 6, 7, 8, 9}.For P > 1, when there are more than one polygon, B is calculated by multiplying the elements of H with the number of polygons P .For example, for P = 3, B is chosen such that B ∈ H = {6, 9, 12, 15, 18, 21, 24, 27}.This makes a total of 4 × 8 = 32 combinations of CMWP-B test instances which are generated for each source instance of the CMWP.Further, two different strategies are followed to generate instances that are grouped as regular and random instances.Regular instances contains polygons with the same number of extreme points.For example, an instance is called regular when all barriers are triangles.On the other hand, random instances contains polygons each of which do not necessarily have the same number of extreme points as long as their total does not exceed B. Shortly, each of regular and random instance groups for the CMWP-B contain 4×8 = 32 test instances.This makes 11 × 32 × 2 = 704 and 40 × 32 × 2 = 1920 CMWP-B test instances that are constructed using existing CMWP test instances Sherali et al. [19] and Boyacı [48] instances, respectively.

Computational Results
We have employed the best known values of the CMWP as our benchmark values as given in Table 1 to compare the performance of our DB heuristics.The percent deviations of the heuristics are determined using the following formula.81055.87b a Results from [27] b Results from [48] where Z U B is the upper bound value found by the suggested heuristic for the CMWP-B.Z best CM W P is the best known objective value for standard CMWP test instances as shown in Table 1.Clearly, the formula ( 25) yields an upper bound on the performance of the heuristics.The instances are divided into two groups as small and large instances depending on their number of facilities and customers.For example, large instances contain I > 10 facilities to locate.In addition, instances with J ≥ 100 customers are also considered as large instances.The upper bounds are determined with the ALAB heuristic as described in Section 4.1 for our DB heuristics.The DB-I heuristic employs the ALAB heuristic as an improvement step by using the facility locations obtained after the proposed MILP formulation is solved.For the DB-II heuristic, an initial upper bound is needed to increase the efficiency of the subgradient algorithm.Initially, facilities are assumed to have unlimited capacity and they can meet all customer demand.Then, facilities are assigned to the candidate locations from which they serve the customers with least total cost.Once the corresponding candidate locations are fixed as the starting facility locations, the ALAB heuristic is run initialized from those facility locations.Lastly, the resulting upper bound is used as the initial upper bound of the subgradient algorithm for the DB-II heuristic.Besides, the ALAB heuristic is systematically run (i.e.once in ten iterations) to update best upper bound found within the subgradient algorithm.Finally, the ALAB heuristic is applied on the best solution found by the TS algorithm within the DB-III heuristic.Table 2 gives the outcomes of the DB heuristics on both random and regular CMWP-B test instances generated from the instances by Sherali et al. [19].The first three columns explain instance properties.The first column gives the name of the standard CMWP instance as mentioned.(I, J) stands for the number of facilities and customers in the second column.The third column states the total number of barrier regions P in the instances."% Dev." and "CPU" denote the percent deviation calculated by the formula (25) and the CPU time of the corresponding DB heuristic in seconds, respectively.To be precise, the cells under these columns are the average of 8 test instances having different number of total extreme points from the set H depending on the number of barriers P .We should emphasize that all instances in Table 2 are small instances.The last row indicates the average of each column.Starting with column 4, six columns are consecutively dedicated to the outcomes of the DB heuristics for random and regular test instances.Best percent deviations at each row are marked in bold characters.On these instances, we observed that, DB-I heuristic outperforms both DB-II and DB-III heuristics in terms of accuracy.Indeed, the DB-I heuristic yields outstanding accuracy for both random and regular test instances.Using LR scheme as in the DB-II heuristic increases the efficiency and yields outcomes in almost half time of the DB-I heuristic.Additionally, the DB-II heuristic produces similar accuracy to the DB-I heuristic in 12 out of 44 cases of random instances.This value is 11 out of 44 cases when regular instances are considered.
The tabu tenure parameter α for the DB-III heuristic is set to α = max{I, 20} in the light of our preliminary experiments.A maximum number of tabu iterations θ is set to 300.Another stopping condition of 30 non-improving tabu iterations is also imposed to avoid from unnecessary computations for the DB-III heuristic.A neighborhood search width denoted with N is empirically determined as N = min{⌈J/3⌉, 30}.Here, ⌈a⌉ is the smallest integer value that is larger than or equal to a.These TS settings are employed in all of our computational experiments.Unfortunately, the performance of the DB-III heuristic is not promising on these small instances.Nevertheless, running time of the DB-III heuristic is shorter than the DB-I heuristic on the average for small instances.
Although the instances based on Sherali et al. [19] are small instances, their sizes are of limited variety.Therefore, Table 3 shows additional results obtained for random and regular CMWP-B instances (small instances) based on Boyacı [48] instances.Diversity of these instances is higher than the ones by Sherali et al. [19].Table 3 can be read similar to Table 2 since they share the same structure.The results of Boyacı [48] instances are similar to Sherali et al. [19] instances and strengthen the success of the DB-I heuristic for small instances.The DB-I heuristic is the winner for both random and regular CMWP-B instances over all test instances in terms of accuracy.In particular, the DB-I heuristic yields percent deviations with a difference of 10% more than that of the closest approach, the DB-II heuristic, on the average for small instances.For these small instances, it becomes evident that the efficiency of the DB-I heuristic decreases drastically, i.e. 110.6 and 108.96 seconds on the average for random and regular instances, respectively.The performance of the DB-II heuristic using LR scheme outperforms the DB-I heuristic in CPU times.Moreover, the DB-II heuristics yields the best results in 18 out of 88 cases for the random instances.Analogously, the DB-II heuristic finds the best results in 17 out of 88 cases for the regular instances.The DB-II heuristic is slightly more efficient than the DB-III heuristic with CPU times of 0.75 (0.74) and 0.98 (0.98) seconds on the average, respectively, for random (regular) instances.The accuracy of the DB-III heuristic is poor for small instances.The efficiency of the DB-I heuristic deteriorates especially for the instances with more than 30 customers (or similarly 30 candidate facility locations).The DB-II and DB-III heuristics show different behavior than the DB-I heuristic does since they are not significantly affected from number of customers.Observe that increasing the number of facilities deteriorates the CPU time of all DB heuristics.Table 4 presents the results of the large CMWP-B instances based on Boyacı [48] instances.The

Conclusion
In this work, we have focused on the capacitated multi-facility Weber problem with polyhedral barriers (CMWP-B).We have suggested a mathematical formulation for the CMWP-B and its discrete equivalent that is a MILP problem.
Discretization based heuristic procedures for the CMWP-B have been proposed.We have carried out extensive computational experiments on randomly generated test instances that are based on standard CMWP instances.The proposed methods compute upper bounds for the optimal objective value of the CMWP-B.
The first discretization based heuristic solves the mixed integer linear programming formulation using the customer locations as the set of candidate facility locations.Its efficiency is improved by using a Lagrangean relaxation scheme and subgradient algorithm that resulted in the second discretization based heuristic.Lastly, the third discretization based heuristic employs a tabu search algorithm using a neighborhood search over customer locations.Among the discretization based heuristics, the first one yields the highest accuracy.However, its performance is limited with relatively small instances.As a remedy, we have applied a Lagrangean relaxation scheme within the second discretization based heuristic which constitutes a compromise between accuracy and efficiency.The performance of the third discretization based heuristic is usually poor, however, its performance is promising for large instances.Implementing exact solution procedures can be a good direction of research for the CMWP-B in the future.Efforts to provide an effective branch and bound algorithm is of high importance for this type of location-allocation type problems.Last but not least, a probabilistic extension of the CMWP-B, where the customer locations and/or their demand quantities are stochastic, is a worthwhile open research area.

Figure 1 .
Figure 1.An illustrative example with 4 customers, a single tetragon barrier and the corresponding visibility graph.

Figure 2 .
Figure 2. Two possible cases for the output of a Weiszfeld-like procedure.
1 at the next tabu iteration.Once tabu facility i * is set at iteration t, the candidate point k * is searched over the closest neighbors of the current candidate location, say k ′ , of facility i * .Clearly, all facilities other than i ′ , and so forth.The upper bound obtained at tabu iteration t and the best upper bound found so far are respectively represented with Z t U B and Z best T abu .The upper bound Z t U B is determined as the lowest upper bound value obtained by exchanging k ′ with the candidate locations of the set {k ′ (1) , k ′ (2) , . . ., k ′ (N ) } one by one from the closest (k ′ (1) ) to farthest (k ′ (N ) ) neighbor candidate location.Here, N is the number of neighbor solutions checked (or the width of the neighborhood) and k ′ (N ) is the N th closest candidate location to k ′ .This strategy is followed as long as Z t U B does not improve the best upper bound Z best T abu .When an improvement is obtained, the neighborhood search is stopped.The best upper bound Z best T abu is updated whenever the improvement occurs and the TS proceeds to the next iteration t + 1 with the next facility.For example, when the feasible solution obtained by exchanging k ′ with the closest candidate location k ′ (1) yields a better upper bound than the best upper bound Z best T abu , then Z t U B is calculated with candidate location k ′ * maintain their facility locations determined at previous iteration t − 1 and the new feasible solution only changes the location of facility i * at iteration t.Let k ′(1) be the closest candidate location to k ′ , k ′(2) denote the second closest candidate location to k

Table 1 .
Benchmark values for the standard CMWP test instances from the literature.

Table 2 .
[19]performance of the DB heuristics over random and regular CMWP-B instances based on Sherali et al.[19]instances.

Table 3 .
[48]performance of the DB heuristics over small sized random and regular CMWP-B instances based on Boyacı[48]instances.