A greedy variable neighborhood search heuristic for the maximal covering location problem with fuzzy coverage radii

Knowledge-Based Systems 41 (2013) 68–76

Contents lists available at SciVerse ScienceDirect

Knowledge-Based Systems journal homepage: www.elsevier.com/locate/knosys

A greedy variable neighborhood search heuristic for the maximal covering location problem with fuzzy coverage radii Soheil Davari a, Mohammad Hossein Fazel Zarandi a,⇑, I. Burhan Turksen b,c a

Department of Industrial Engineering, Amirkabir University of Technology, 424 Hafez Avenue, Tehran, Iran Department of Mechanical and Industrial Engineering, University of Toronto, Toronto, Ontario, Canada M5S 2H8 c TOBB Economics and Technology University, Sogutozu Caddesi No: 43, Ankara, Turkey b

a r t i c l e

i n f o

Article history: Received 2 July 2012 Received in revised form 6 November 2012 Accepted 25 December 2012 Available online 3 January 2013 Keywords: Maximal covering location problem (MCLP) Credibility theory Variable neighborhood search Fuzzy simulation Facility location

a b s t r a c t The maximal covering location problem (MCLP) seeks location of facilities on a network, so as to maximize the total demand within a pre-defined distance or travel time of facilities (which is called coverage radius). Most of the real-world applications of MCLP comprise many demand nodes to be covered. Moreover, uncertainty is ubiquitous in most of the real-world covering location problems, which are solved for a long-term horizon. Therefore, this paper studies a large-scale MCLP on the plane with fuzzy coverage radii under the Hurwicz criterion. In order to solve the problem, a combination of variable neighborhood search (VNS) and fuzzy simulation is offered. Test problems with up to 2500 nodes and different settings show that VNS is competitive, since it is able to find solutions with gaps all below 1.5% in much less time compared to exact algorithms. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction and problem description The problem of locating facilities is not new to the operations research community. The challenge of where to best site facilities has inspired a rich, colorful and ever growing body of literature. Facility location has a decisive role in success of supply chains with numerous applications in locating retail stores, hospitals, production facilities, emergency units, etc. One of the traditional location problems, which has been well studied since the very beginning of location science, is the covering location problem. In a covering location problem, one seeks a solution to cover a subset of customers considering one or more objectives. The covering location problem is often categorized as location set covering problem (LSCP) and maximal covering location problem (MCLP). While the LSCP tries to cover all the demand using the least possible number of facilities, the covered population is maximized in MCLP using a fixed number of facilities. It is to be noted that in both of these problems, a node is covered if there is at least one facility within a defined distance to it. Such a distance is called coverage radius and should be determined carefully. The study of LSCP and MCLP root back to 1970s and initially presented by Toregas et al. [1], and Church and ReVelle [2] respectively.

⇑ Corresponding author. Tel.: +98 21 64 54 53 78; fax: +98 21 66 41 30 25. E-mail address: [email protected] (M.H. Fazel Zarandi). URL: http://cil.aut.ac.ir (M.H. Fazel Zarandi). 0950-7051/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.knosys.2012.12.012

An omnipresent feature of most of the location decisions is the uncertainty associated with their parameters such as demands and travel times. A possible approach of modeling location problems under uncertainty is using stochastic models which are futile when there is not enough historical data to be used in the model or when gathering reliable historical data is not reasonable, especially for larger datasets. Another well-known approach is modeling problems using fuzzy logic. Moreno Pérez et al. [3] claimed that in real applications, the facility locations can be full of linguistic vagueness. Such vagueness can be appropriately modeled using networks with the fuzzy values which describe nodes, weights or importance of nodes, lengths of paths, etc. Like that, the knowledge of experts could be elicited as fuzzy membership functions, such as representing the process duration as ‘‘something between 20 and 30 min’’ or the demand as ‘‘about 100 units per week’’. This paper considers a problem where coverage radii are not precisely known and treated as fuzzy variables. Some applications are locating warning sirens, locating lighting posts to illuminate an area, or even radio stations. In these cases, the coverage depends on some elements which could not be modeled easily for long-term horizons. Take the example of locating wireless internet facilities to provide service for the people in an area. The coverage of demand nodes depends on network traffic, weather conditions, topography of the area, and some other factors which vary by time. The same happens in locating lighting facilities in highways, where weather conditions influence the coverage status to a great extent. These two instances are major applications of the problem discussed in this paper.

69

S. Davari et al. / Knowledge-Based Systems 41 (2013) 68–76

In the literature, there exist many approaches to solve the crisp version of MCLP. These include exact, heuristic, and metaheuristic methods which will be discussed in the next section. It is known that solution of large instances of MCLP is unmanageable using exact methods. Therefore, heuristics and metaheuristics have been employed in order to solve MCLPs of larger sizes. ReVelle et al. [4] presented using heuristic concentration to solve MCLP in a crisp environment. They solved MCLPs with 900 nodes. In this paper, we will solve problems with 1000, 1500, and 2500 nodes and will show the performance of variable neighborhood search (VNS) on a set of test problems. The main contributions of this work are as follows. First, the concept of fuzzy coverage radii is presented and the chanceconstrained model of fuzzy MCLP (FMCLP) is offered. Moreover, a hybrid algorithm of fuzzy simulation and VNS is used to solve the fuzzy version of the problem. Besides, large-scale MCLPs are solved using the proposed algorithm and some analysis is given. The rest of the paper is organized as follows: First, a concise literature review of covering problems and related issues are presented in Section 2. Then, fuzzy variables and basics of credibility theory are discussed. Section 4 is dedicated to FMCLP and the chance-constrained version of MCLP. The proposed solution algorithm is presented in Section 5 and numerical examples appear in Section 6. Moreover, results are analyzed and discussions are given in this section. Finally, to bring the paper to a close, conclusions and outlooks for potential future research are given in Section 7.

2. Literature review The literature of covering models is too diverse to be completely studied in this paper. Since this paper does not deal with a comprehensive review of the literature in the area of MCLP, for further information, one may refer to valuable reviews such as ReVelle and Eiselt [5] and Berman et al. [6]. In what follows, some major contributions to the literature of MCLP are given. As already stated, MCLP was first introduced by Church and ReVelle [2]. Later, various extensions to the basic model have been proposed. The hierarchical MCLP is studied by Moore and ReVelle [7]. Plastria and Vanhaverbeke [8] addressed MCLP in competitive environments. Gradual covering problems are studied in publications such as Berman et al. [9], Berman and Krass [10], and Berman et al. [11]. Younies and Wesolowsky [12] introduced a zero-one mixed integer programming formulation for a maximal covering location problem where points are covered by inclined parallelograms in a plane. The computational hurdle of many location models limited research to deterministic cases which are usually very hard to solve when the problem size grows. However, recently the ability to solve combinatorial problems has been improved dramatically, and some scholars proposed various approximate procedures to deal with MCLP. There is some evidence that local search metaheuristic have been more popular in dealing with MCLPs. Another stream of research has been devoted to real-world case studies. To solve a covering problem, Murawski and Church [13] proposed a model to assume that the established facilities are fixed and the accessibility of demand nodes are to be improved. Their model is called maximal covering network improvement problem which was formulated as an integer-linear programming problem and a case study in Ghana is surveyed. Curtin et al. [14] studied the maximal covering location problem in order to determine the distribution of police patrol areas in Dallas. Recently, Basar et al. [15] studied the multi-period MCLP to locate emergency medical services in Istanbul, Turkey. The first wave of published location models were deterministic and, thus, do not account for uncertainties, which may exist in

Table 1 Some problems solved in fuzzy environment using credibility theory. Paper

Problem

Solution procedure

Peng and Liu [19] Zhao and Liu [20]

Parallel machine scheduling Standby redundancy optimization Vehicle routing problem Quadratic assignment problem Portfolio selection Location-routing problem

Genetic algorithm Genetic algorithm

Zheng and Liu [21] Liu and Li [22] Huang [23] Fazel Zarandi et al. [24] Ke and Liu [25] Wen and Iwamura [26] Davari et al. [18] Fazel Zarandi et al. [27]

Project scheduling Location-allocation problem Maximal covering location problem Location-routing problem

Genetic algorithm Genetic algorithm Genetic algorithm Simulated annealing Genetic algorithm Genetic algorithm Simulated annealing Simulated annealing

location problems, such as covering location variants. Some examples of uncertainty in covering location problems are the uncertainty associated with coverage radius of a facility or the demand of a node to be covered. The uncertain MCLP has been addressed by some scholars such as Batanovic et al. [16] who suggested maximal covering location problems in networks with uncertainty. They studied problems with equal importance of demand nodes, relative deterministic weights of demand nodes, and linguistic terms for weights of demands. In a publication by Berman and Wang [17], they surveyed a case in which the demands are random variables with unknown probability distributions. Davari et al. [18] is another recent attempt in modeling covering problems in the fuzzy sense. They modeled a closed-loop supply chain network and used fuzzy goal programming to solve it. Employing hybrid algorithms to solve problems with fuzzy parameters is not new to the literature. Table 1 presents some of the recent publications in the literature concerning various mathematical programming problems which have been solved using hybrid algorithms of metaheuristics and fuzzy simulation. Our review of the MCLP literature clearly shows that uncertain covering problems have not gained much attention so far. To the best of our knowledge, the only publication dealing with uncertain covering location problem using fuzzy variables is Davari et al. [18]. Moreover, although heuristics and metaheuristics have been a major tool in dealing with MCLPs, the performance of VNS needs to be examined further. Based on these findings, this paper contributes to the literature by filling these two gaps by proposal of a fuzzy MCLP and also presenting an effective solution procedure to solve it. 3. Fuzzy variables Fuzziness is a possible type of uncertainty which was introduced by Zadeh [28]. In order to measure fuzzy events, the concept of possibility measure was presented in 1978. Later, Liu and Liu [29] presented the credibility measure which is a self-dual measure. This section presents some basic concepts, axioms, and definitions in credibility theory to be used in this paper. Definition 1 30. A fuzzy variable is defined as a measurable function from the credibility space (H, P, Cr) to the set of real numbers.

Definition 2 30. Let n be a fuzzy variable defined on the credibility space (H, P, Cr). Then, its membership function is derived from the credibility measure as:

70

S. Davari et al. / Knowledge-Based Systems 41 (2013) 68–76

lðxÞ ¼ ð2Crfn ¼ xgÞ ^ 1 x 2 R

ð1Þ

where ^ is the minimum operator. Definition 3 (30). (Credibility inversion theorem) Let n be a fuzzy variable with membership function l. Then for any set B of real numbers, we have:

1 suplðxÞ þ 1  suplðxÞ Crfn 2 Bg ¼ 2 x2B x2Bc

! ð2Þ

which is a self-dual measure. (Possibility and necessity measures lack the self-duality property.) Now let us consider an example of a triangular fuzzy variable n = (r1, r2, r3). From the definitions of possibility, necessity and credibility, it is easy to obtain:

Posfn P rg ¼

Necfn P rg ¼

Crfn P rg ¼

8 > <1

r 3 r

r 3 r2 > : 0

8 > <1 > :

if r 2 6 r 6 r 3

if r 6 r 1

0

if r P r 2

2ðr 2 r 1 Þ

ð3Þ

if r P r3

if r1 6 r 6 r 2

r3 r > > 2ðr3 r2 Þ > > : 0

-nsup ðaÞ þ ð1  -Þninf ðaÞ

ð8Þ

In Hurwicz criterion, the parameter x 2 ð0; 1, which reflects the degree of the decision maker’s optimism, must be determined by the decision maker a-priori. The proper value of x varies from person to person. By varying the parameter x, the Hurwicz criterion becomes various criteria, e.g., when x = 1, the criterion is the optimistic criterion; when x = 0, it degenerates to a pessimistic criterion. 4. MCLP with fuzzy coverage radii

if r 6 r2

r 2 r r 2 r 1

8 1 > > > 2r2 r1 r > <

ninf(a) > ginf(a), where ninf ðaÞ and ginf ðaÞ are the a-pessimistic values of n and g, respectively. In this paper, a combination of a-optimistic and a-pessimistic values will be used to rank solutions. Such a combination is performed using the Hurwicz criterion [31] which could be stated as follows:

ð4Þ

if r 6 r 1 if r 1 6 r 6 r 2 if r 2 6 r 6 r 3

ð5Þ

if r P r 3

Definition 4 (30). Let n be a fuzzy variable, and a 2 ð0; 1. Then, (6) and (7) are respectively called a-optimistic and a-pessimistic values to n.

nsup ðaÞ ¼ supfrjCrfn P rg P ag

ð6Þ

ninf ðaÞ ¼ inffrjCrfn 6 rg P ag

ð7Þ

This means that the random variable n will reach upwards of the aoptimistic value at least a of time, and will be below the a-pessimistic value at least a of time. It should be noted that a is a parameter defined by decision maker(s) and directly affects the solution(s).

Before addressing MCLP with coverage radii, the model assumptions are given. It should be noted that the classical formulation of MCLP which was presented in Revelle et al. [4] has been used in this paper. In this formulation, let I be the set of demand nodes, and J be the set of eligible facilities. The problem is studied on a plane. It is assumed that the cost to establish facilities is identical for all facilities. Thus, the fixed cost of establishing facilities is excluded from the model. For each potential facility, there is a coverage radius which determines the coverage status of nodes. It is to be highlighted that coverage radii are independent fuzzy variables. The distance between nodes is a symmetric matrix. Each demand node has a specific demand to be met and proximity of demand nodes to a facility is desirable. The ultimate goal of the problem is to maximize the total demand covered by establishing a definite number of facilities. In the sequel, the crisp version of MCLP is presented as follows: Problem parameters hi p dij b Ni

Demand of node i (The weight of ith node) Number of facilities to be located Distance/time between nodes i and j Coverage radius {j|dij 6 b} = the nodes j that are within the coverage radius of node i

Decision variables Theorem 1 (30). Let n be a fuzzy variable. Then, we have:

 xj ¼

(a) ninf (a) is an increasing and left-continuous function of a. (b) nsup(a) is a decreasing and left-continuous function of a. Lemma 1 (30). Suppose that fuzzy variables n and g are independent. Then, for any a 2 ð0; 1, we have: (a) (b) (c) (d)

(n + g)sup(a) = nsup(a) + gsup(a). (n + g)inf(a) = ninf(a) + ginf(a). if c P 0, then (cn)sup(a) = cnsup(a) and (cn)inf(a) = cninf(a). if c 6 0, then (cn)sup(a) = c ninf(a) and (cn)inf(a) = cnsup(a).

In order to rank fuzzy variables, Liu [30] proposed four criteria including optimistic and pessimistic value criteria. Optimistic value criterion [30]. We say n > g if and only if, for some predetermined confidence level a 2 ð0; 1, we have nsup(a) > gsup(a), where nsup(a) and gsup(a) are the a-optimistic values of n and g, respectively. Pessimistic value criterion [30]. We say n > g if and only if, for some predetermined confidence level a 2 ð0; 1, we have

 yi ¼

1 if a facility is established in node j 0 otherwise 1 if node i is covered by at least one facility 0

otherwise

Now the problem of maximal covering location problem (MCLP) is stated as follows:

Max z ¼

X hi yi

ð9Þ

i2I

yi 6

X

xj

i2I

ð10Þ

j2Ni

X xj ¼ p

ð11Þ

j2J

0 6 yi 6 1 i 2 I

ð12Þ

xj 2 f0; 1g j 2 J

ð13Þ

S. Davari et al. / Knowledge-Based Systems 41 (2013) 68–76

Objective function (9) states that the objective is to maximize the total covered demand. Constraint (10) guarantees that the demand of a node is covered if there is at least one facility located with a distance less than the coverage radius to it. Constraint (11) holds that there are p facilities to be located. Constraint (12) states that yi does not need to be an integer in the model. In other words, when x variables are integers, the y variables always take an integer value. Finally, constraint (13) imposes binary restriction on decision variables. Since the coverage state of a node is determined using the coverage radii of facilities, it is obvious that coverage radii have an implicit role in fitness of a solution. Clearly, since coverage radii are considered to be fuzzy variables, there is a need to rank solutions based on a rational ranking criteria of fuzzy variables. Hence, the Hurwicz criterion is applied which has been already presented. Eqs. (14) and (15) hold, following the independence of fuzzy variables.

(

) ( ) X sup rjCr hi yi P r P a ¼ ¼

X

X hi y i

i

!

i

ðaÞ sup

ðhi yi Þsup ðaÞ

ð14Þ

i

(

) ( ) X hi yi 6 r P a ¼ inf rjCr i

X hi yi

! ðaÞ

i

¼

inf

X ðhi yi Þinf ðaÞ

ð15Þ

i

Now, the chance-constrained version of MCLP is formulated as (16):

max fx  X xj yi 6

X X ðhi yi Þsup ðaÞ þ ð1  xÞ  ðhi yi Þinf ðaÞg i

i

i2I

j2Ni

X xj ¼ p

ð16Þ

j

0 6 yi 6 1 i 2 I xj 2 f0; 1g j 2 J

x 2 ½0; 1 5. Solution procedure of MCLP and FMCLP The proposed procedure in this paper can be considered as a juxtaposition of VNS and fuzzy simulation. While VNS is used to explore the search space as a tractable approach, fuzzy simulation performs the simulation of each solution to give an approximate value for its fitness. To put it in simpler terms, a hybrid intelligent algorithm has been designed integrating fuzzy simulation and VNS to solve the fuzzy model of this paper. As is shown in (16), the fitness of each solution should be estimated as a weighted average of pessimistic and optimistic values. To estimate the optimistic and pessimistic values of a solution, there is a need to use fuzzy simulation as will be shown in Figs. 1 and 2. Then, VNS is integrated

Step 1. Randomly generate θk from the credibility space (Θ, P, Cr), write υk=(2Cr{θk})∧1

71

Step 1. Randomly generate θk from the credibility space (Θ, P, Cr), write υk=(2Cr{θk})∧1 and produce ξk = ξ(θk),k=1,2,…,N, respectively. Equivalently, randomly generate ξk and write υk=μ(ξk) for k=1,2,…,N, where μ is the membership function of ξ. Step 2. Find the minimal value r such that L(r)≥α holds where L(r) is found as: 1 L(r ) = (max{ν k | f ( x, ξ k ) ≤ r} + min{1 −ν k | f ( x, ξ k ) > r}) 1≤ k ≤ n 2 1≤k ≤n

Step 3. Return r. Fig. 2. Fuzzy simulation to estimate U2(x).

with fuzzy simulation to produce a hybrid intelligent algorithm to solve the problem. Since a greedy approach is carried out in initialization and also the local search procedure of the proposed approach, we call our procedure as a greedy VNS. 5.1. Review of solution procedure modules In recent years, stochastic-based search algorithms have been widely used to find optimum solutions to academic and practical problems [32]. A clear indication for the popularity of these methods is the large number of publications in the literature which have used one or some of these methods for various problems. Variable neighborhood search (VNS) is a recently developed metaheuristic which looks for improvements in solution quality by systematically changing the neighborhood structure. It has shown its ability to solve many problems such as: scheduling [33], vehicle routing problem [34], and generalized assignment problem [35]. The main steps of the basic VNS are cited below:  Step 1. Select the set of neighborhood structures Nk (k = 1, 2, . . . , kmax) that will be used in the search. Find an initial solution x and choose a stopping criterion or a set of stopping criteria.  Step 2. Repeat the following steps until the stopping criterion is met: (1) Set k = 1. (2) Repeat the following steps if k 6 kmax:  Shaking Randomly generate a solution x’ from the kth neighborhood of x.  Local Search Apply some local search method with x’ as the initial solution, and denote the so obtained local optimum as x’’. It is worth mentioning that in this paper, to find the local optimum, the local search is run for a number of pre-defined iterations and the best found solution is regarded as the local optimum.  Move or not If x00 is better than the incumbent, update x = x00 , and set k = 1; otherwise set k = k + 1. As already stated, the Hurwicz criterion is applied in this paper to rank fuzzy variables. Liu [30] proposed an algorithm to compute the uncertain functions of the form U 1 : x ! maxff jCrff ðx; nÞ P f g P ag and U 2 : x ! minff jCrff ðx; nÞ 6 f g P ag which is presented in Figs. 1 and 2. Take note that contrary to deterministic problems, the fitness value of each solution in VNS should be defined using the procedures in Figs. 1 and 2.

and produce ξk = ξ(θk),k=1,2,…,N, respectively. Equivalently, randomly generate ξ k and write υk=μ(ξk) for k=1,2,…,N, where μ is the membership function of ξ.

5.2. Initialization of solutions

Step 2. Find the maximal value r such that L(r)≥α holds where L(r) is found as: 1 L(r ) = (max{ν k | f ( x, ξ k ) ≥ r} + min{1 −ν k | f ( x, ξ k ) < r}) 1≤ k ≤ n 2 1≤k ≤n

Step 3. Return r. Fig. 1. Fuzzy simulation to estimate U1(x).

Obviously, metaheuristics are quite sensitive to their initial solutions. Therefore, initialization of good solutions plays a pivotal role in reaching better solutions. To be resourceful in finding good initial solutions of VNS, a greedy heuristic will be proposed which shows promising results.

72

S. Davari et al. / Knowledge-Based Systems 41 (2013) 68–76

In order to get a clear insight of the problem and its characteristics, first a couple of test problems were solved to optimality. Then, we came to notice that in optimal solutions, there is a low level of intersection between covered regions. Therefore, we take advantage of such a finding and propose to generate initial solutions with low degree of intersection. To do so, the total coverage of a node is found adding up the demands of nodes with a distance less than the coverage radius to it. Then, the first facility to be established is found using the roulette wheel selection. In other words, the node with the highest amount of coverage has a higher chance of being selected. All the covered nodes by the established facility are removed from the space and the coverage of each node is re-calculated using the remaining nodes. This procedure is followed until p facilities are established. Fig. 3 shows that (a) covers 15 nodes. Moving one of the circles to the right, some other nodes are added to the list of covered nodes, while all the covered nodes in (a) remain covered. Such a change increases the number of covered nodes from 15 to 21. This shows the role of coverage intersection on the quality of solutions. An additional module of the proposed algorithm is the re-start module. This module is activated whenever a specific region of the search space is exhaustively explored without reaching a better solution for a number of consecutive runs. Such an approach is used in order to overcome the possibility of getting stuck in local optima. The re-start mechanism is to pursue the search from a new solution which is generated using the initialization procedure. Our preliminary experiments clearly showed the effectiveness of this module in improving the solution fitness. 5.3. Encoding scheme An effective encoding scheme has a considerable impact on the performance of VNS. A possibility to encode the solutions of our problem is to use binary representation to show establishment of facilities in candidate nodes. Considering n nodes available and p facilities to be established, a solution is shown using a vector of binary elements with n bits where exactly p bits take a one value. The procedures of encoding and decoding a candidate solution are illustrated by applying them to an example. Consider a problem with eight candidate nodes to establish a facility in which four facilities are to be located. A feasible solution looks like [1 0 1 0 1 0 1 0] which shows that four facilities are to be located at 1st, 3rd, 5th, and 7th nodes. One special feature of the presented MCLP is that any generated solution is feasible, as long as p facilities are located. This special characteristic of model is a great advantage which facilitates searching the solution space. Since the sizes of problems targeted in this paper are very large, there is a need to devise a method to

Fig. 4. The performance of shake(5) on a sample solution with n = 20, p = 8.

calculate the fitness of a solution in a fast and easy way which is derived from ReVelle et al. [4] and is shown below.

z¼

X

hi  Min½1;

i2I

X xj 

ð17Þ

j2Ni

Such a procedure reduces the computational load to a great extent and has shown its effectiveness in various test problems. 5.4. Neighborhood search structures The set of neighborhoods used for shaking is at the heart of the VNS. Each neighborhood should strike a proper balance between perturbing the incumbent solution and retaining the good parts of the incumbent solution [36]. In this paper, a shaking phase is carried out whenever k random open facilities are closed and k close facilities are randomly selected and opened simultaneously. In other words, a shake (k) move is associated with generating solutions in which the set of open facilities differ from the current solution in k indices. Since there are exactly p open facilities in a sample solution, the value of k in the shaking phase should be below p. Based on some preliminary runs, kmax is determined to equal 10. To select the minimum value of k, it should be noted that the shaking phase should provide a fair amount of diversification in solutions. Thus, the minimum value of k is equal to 5 and there are six possible shaking types to be applied in this paper. Fig. 4 depicts a sample shaking move using k = 5 where the changed bits are highlighted. To design an effective and efficient local search procedure, it has been taken into account that an effective local search for MCLP should be of a greedy nature. Therefore, a greedy two-opt move is selected as the local search module of the proposed algorithm. This move operates as follows: First a list of open facilities is constructed and the total coverage of each is computed. Moreover, the potential coverage for each of the non-open facilities is found. Then, a node from the non-open facilities replaces an open facility. It is to be noted that selection of these two nodes is carried out using the roulette wheel selection (RWS) mechanism. This increases the chance of selection for closed facilities with higher and open facilities with lower coverage levels. Undoubtedly, the proposed local search module is greedy and stochastic. 6. Numerical examples 6.1. Test problems

Fig. 3. The effect of coverage intersection on the solution fitness.

To generate the demands of the test problems, the same procedure as ReVelle et al. [4] has been followed, where demand of each node takes a value following the random distribution between 0 and 100. In addition, three different types of demand distribution on the plane have been employed in generating test problems. These three sets differ in the distribution of demand nodes on the plane, where the spread of nodes follow gamma, normal, and uniform distributions. The reason to select these three types is the fact that while the gamma distribution is a skewed distribution, the other two are symmetrical. Moreover, normal and uniform distributions are different considering the fairness of distribution of nodes on the plane. Thus, these three kinds of

S. Davari et al. / Knowledge-Based Systems 41 (2013) 68–76

73

Fig. 5. Three kinds of demand distribution on the plane.

distributions are able to represent three various types of demand locations. Fig. 5 depicts these three types of distribution using the data of this paper. For each of these demand types, instances with n = 1000, 1500, 2500, b = 2, 2.5, 3, and p = 20, 30, 40 are solved. Therefore, to examine the performance of the proposed VNS, 81 instances were solved. It should be noted that while ReVelle et al. [4] solved problems with 900 nodes, we considered problems of larger sizes up to 2500 nodes. The crisp problems were solved to optimality using CPLEX. Then, the fitness and runtime of running the proposed VNS on the same problems were compared with the solutions obtained from CPLEX. In the fuzzy sense, the same set of problems is solved. The only difference between the crisp and fuzzy versions of the problem is that in the fuzzy version of the problem, coverage radii are assumed to be fuzzy variables. To put it another way, the coverage radius of each facility is assumed to be a fuzzy variable with a triangular membership function. The parameters of these variables are (2, 2.25, 3) for each facility. Such an assumption has been already justified in the introduction of this paper. 6.2. Parameter setting As already stated, six types of shaking are employed in this paper. Moreover, in the local search stage, the greedy RWS-based

procedure is applied. Our experiments showed that running the algorithm for more than 100 s is not effective in nearly all the test instances. Experiments showed that in about one percent of runs, there has been some improvements in the solution quality, which is not considerable. Therefore, reaching 100 s of runtime has been selected as the stopping criterion of the algorithm. Furthermore, to estimate critical values, 10,000 solutions were generated. 6.3. Effect of using search structures An analysis has been carried out to examine the role of using the proposed shaking types on the performance of the proposed procedure. To do so, six replicates of each of the scenarios have been solved by excluding each operator from the proposed VNS, while other parameters of the procedure do not change. Fig. 6 shows the gaps of these experiments graphically using a box-plot diagram, where Si is the result of excluding the ith move from the solution procedure. It is to be noted that the S7 is the one in which no operator is excluded. It is interesting to see that S7 is the best scenario in terms of solution quality. Moreover, since the runtimes of the algorithm are identical for all of the scenarios, there is not any difference between the scenarios regarding their runtime. Hence, the multi-operator setting is both effective and efficient and employing S7 is justified, based on the numerical experiments.

Fig. 6. Box-plot diagram of gaps associated with each scenario.

74

S. Davari et al. / Knowledge-Based Systems 41 (2013) 68–76

6.4. Results, validation, and discussions

Table 4 The results of VNS on the problems with uniform distribution of nodes.

Test Problems have been solved using two approaches. First, each problem was solved using the CPLEX solver. Then, solutions were compared versus the results of the proposed VNS. Tables 2– 4 summarize the results of comparing the performance of the pro-

Table 2 The results of VNS on the problems with gamma distribution of nodes. n

1000 1000 1000 1000 1000 1000 1000 1000 1000 1500 1500 1500 1500 1500 1500 1500 1500 1500 2500 2500 2500 2500 2500 2500 2500 2500 2500

p

20 20 20 30 30 30 40 40 40 20 20 20 30 30 30 40 40 40 20 20 20 30 30 30 40 40 40

b

2 2.5 3 2 2.5 3 2 2.5 3 2 2.5 3 2 2.5 3 2 2.5 3 2 2.5 3 2 2.5 3 2 2.5 3

CPLEX

VNS

Gap (%)

Fitness

Time

Fitness

Time

37415.54 42328.46 47292.33 42151.55 45714.90 48546.69 46903.71 49800.75 50144.09 54895.12 64050.10 69838.21 63762.22 71786.50 72818.90 65927.73 70341.32 75238.37 90343.92 102972.79 115236.07 104533.25 115898.82 123446.83 112734.18 123136.64 125826.42

3 21 21 4 17 19 15 13 6 27 156 78 158 37 252 106 24 103 238 1382 1390 1133 1488 1164 1578 3496 1873

37263.91 41846.17 47052.96 41751.98 45484.19 48204.11 46438.67 49360.91 50027.83 54152.28 63831.76 69336.65 63174.57 71776.08 72145.28 65560.92 69836.68 74535.75 90063.45 102025.95 114797.35 103955.35 115488.81 122456.14 111914.87 121333.60 125262.33

100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100

0.41 1.14 0.51 0.95 0.50 0.71 0.99 0.88 0.23 1.35 0.34 0.72 0.92 0.01 0.93 0.56 0.72 0.93 0.31 0.92 0.38 0.55 0.35 0.80 0.73 1.46 0.45

Table 3 The results of VNS on the problems with normal distribution of nodes. n

p

b

1000 1000 1000 1000 1000 1000 1000 1000 1000 1500 1500 1500 1500 1500 1500 1500 1500 1500 2500 2500 2500 2500 2500 2500 2500 2500 2500

20 20 20 30 30 30 40 40 40 20 20 20 30 30 30 40 40 40 20 20 20 30 30 30 40 40 40

2 2.5 3 2 2.5 3 2 2.5 3 2 2.5 3 2 2.5 3 2 2.5 3 2 2.5 3 2 2.5 3 2 2.5 3

CPLEX

VNS

Gap (%)

Fitness

Time

Fitness

Time

39632.05 43598.80 48875.17 44434.79 45028.41 49859.39 48095.96 49991.92 50114.09 58606.98 68827.43 69661.38 68660.78 73380.80 75141.34 69067.58 73186.45 75238.37 94470.48 112968.10 118450.97 111221.62 122796.74 125014.85 111063.77 117125.25 125286.42

6 3 29 42 8 20 24 18 1 208 143 332 309 445 174 12 247 3 2820 807 2371 4246 7428 1492 10,784 2211 2206

39290.44 43259.94 48513.40 44152.77 44868.60 49571.17 47673.64 49293.15 50080.68 58024.86 68518.19 69284.10 68276.75 72911.82 74192.09 68546.19 72789.47 74630.81 93920.11 112716.17 117707.31 110642.97 121883.18 124813.67 110307.75 116387.66 125192.60

100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100

0.86 0.78 0.74 0.63 0.35 0.58 0.88 1.40 0.07 0.99 0.45 0.54 0.56 0.64 1.26 0.75 0.54 0.81 0.58 0.22 0.63 0.52 0.74 0.16 0.68 0.63 0.07

n

p

b

1000 1000 1000 1000 1000 1000 1000 1000 1000 1500 1500 1500 1500 1500 1500 1500 1500 1500 2500 2500 2500 2500 2500 2500 2500 2500 2500

20 20 20 30 30 30 40 40 40 20 20 20 30 30 30 40 40 40 20 20 20 30 30 30 40 40 40

2 2.5 3 2 2.5 3 2 2.5 3 2 2.5 3 2 2.5 3 2 2.5 3 2 2.5 3 2 2.5 3 2 2.5 3

CPLEX

VNS

Gap (%)

Fitness

Time

Fitness

Time

23045.63 30541.76 37483.48 31489.22 40387.27 46919.42 38381.36 46954.32 50134.54 32742.47 41738.17 54065.55 45230.65 57964.91 70062.71 55662.08 68394.82 69913.02 50486.76 68366.69 88943.59 70593.83 94177.60 112384.98 87809.07 108169.31 124312.74

1 1 3 1 2 11 1 4 18 3 3 4 3 8 38 3 26 118 9 14 15 9 13 1566 10 653 2461

22955.63 30339.70 37208.69 31125.87 40270.70 46828.02 38084.19 46668.89 50000.85 32278.07 41511.42 53681.49 44926.89 57617.41 69441.02 55203.91 67478.76 69641.03 50327.35 67556.89 88471.02 69969.38 93269.87 111619.15 86741.00 107883.91 123148.43

100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100 100

0.39 0.66 0.73 1.15 0.29 0.19 0.77 0.61 0.27 1.42 0.54 0.71 0.67 0.60 0.89 0.82 1.34 0.39 0.32 1.18 0.53 0.88 0.96 0.68 1.22 0.26 0.94

posed VNS against exact results obtained from CPLEX. Results indicate that although the proposed algorithm is unable to solve most of the cases to optimality, it works great to solve problems of various sizes and the gaps are all below 1.5%. Additionally, in 71 out of 81 instances solved, the gap is below 1% which equals more than 87% of instances. Since the proposed VNS is able to reach near-optimal solutions in much less time compared to CPLEX, it is reasonable to conclude that it is competitive with CPLEX. The reason behind this statement is the tradeoff between the solution quality and the runtimes and the fact that solution quality is not the only measure of interest, especially in simulation-embedded algorithms. Another characteristic of the proposed VNS is its strength in escaping local optima. Fig. 7 obviously shows the evolution of the solution fitness throughout the algorithm for a sample problem. It is clear that although the improvement in solution quality is almost gradual, it reaches a near-optimal solution in 100 s. It may be argued that since CPLEX reaches optimal solutions in all the instances, there is no need to use a heuristic or metaheuristic. Although this is true for cases of this paper, it should be noted that for larger instances, CPLEX fails to reach optimal solutions in much longer periods of time. Moreover, solution quality is not the only measure of interest. Since a fuzzy simulation will be embedded in the algorithm of crisp version, there is a need to use a tractable and effective solution method for the crisp version. This justifies using metaheuristics to solve MCLP and FMCLP. Tables 2–4 show the average performance of the proposed VNS in ten runs for each test problem and testify that as the problem size grows, generally, it takes more and more time for CPLEX to reach optimal solution. This is due to the combinatorial nature of MCLP. Despite this fact, the runtime of our proposed approach remains quite the same. The efficiency of the proposed algorithm has been evidenced by yielded results. It means that the proposed algorithm works well in a crisp sense. We believe that the performance of the algorithm is quite promising. Thus, fuzzy simulation is embedded in the proposed VNS in order to solve FMCLP.

S. Davari et al. / Knowledge-Based Systems 41 (2013) 68–76

75

Fig. 7. The evolution of the solution using the proposed VNS for n = 2500, p = 30, b = 3 and normal distribution of nodes.

Fig. 8. Four scenarios with different levels of a for n = 1500, p = 20 and gamma distribution of nodes.

Fig. 8 illustrates the solution of a MCLP with fuzzy radii for four different scenarios with x = 1 and a = 0.8, 0.85, 0.9, and 0.95. Results attest that the four scenarios lead to different results.

7. Conclusion and future research areas To sum up, it is to be highlighted that this paper studies a fuzzy maximal covering location problem. A hybrid algorithm of fuzzy

simulation and VNS has been proposed to solve the fuzzy version of MCLP. Moreover, a heuristic was developed to find good initial solutions of VNS. We have shown that while our proposed approach is superior to the exact method in terms of runtime, there are negligible gaps compared to the optimal solutions. Future works could be directed to various issues. A possible future study could be to compare using various heuristics/metaheuristics on this problem. Another avenue for future research could be to assess the performance of VNS for other variants of covering

76

S. Davari et al. / Knowledge-Based Systems 41 (2013) 68–76

location problem such as set covering problem or considering other parameters of the problem as fuzzy variables. Furthermore, FMCLP could be enriched by adding some other assumptions such as gradual covering. In such a problem, our heuristic should be modified to some extent. The authors of this paper are currently working on the modification of the proposed heuristic for the gradual covering location problem. References [1] C. Toregas et al., The location of emergency service facilities, Operations Research 19 (1971) 1363–1373. [2] R. Church, C.S. ReVelle, The maximal covering location problem, Papers of the Regional Science Association 32 (1974) 101–118. [3] J.A. Moreno Pérez, J.M. Moreno Vega, J.L. Verdegay, Fuzzy location problems on networks, Fuzzy Sets and Systems 142 (3) (2004) 393–405. [4] C. ReVelle, M. Scholssberg, J. Williams, Solving the maximal covering location problem with heuristic concentration, Computers & Operations Research 35 (2) (2008) 427–435. [5] C.S. ReVelle, H.A. Eiselt, Location analysis: a synthesis and survey, European Journal of Operational Research 165 (1) (2005) 1–19. [6] O. Berman, Z. Drezner, D. Krass, Generalized coverage: new developments in covering location models, Computers & Operations Research 37 (10) (2010) 1675–1687. [7] G.C. Moore, C.S. ReVelle, The hierarchical service location problem, Management Science 28 (1982) 775–780. [8] F. Plastria, L. Vanhaverbeke, Discrete models for competitive location with foresight, Computers & Operations Research 35 (3) (2008) 683–700. [9] O. Berman, D. Krass, Z. Drezner, The gradual covering decay location problem on a network, European Journal of Operational Research 151 (3) (2003) 474–480. [10] O. Berman, D. Krass, The generalized maximal covering location problem, Computers & Operations Research 29 (6) (2002) 563–581. [11] O. Berman et al., The ordered gradual covering location problem on a network, Discrete Applied Mathematics 157 (18) (2009) 3689–3707. [12] H. Younies, G.O. Wesolowsky, A mixed integer formulation for maximal covering by inclined parallelograms, European Journal of Operational Research 159 (1) (2004) 83–94. [13] L. Murawski, R.L. Church, Improving accessibility to rural health services: the maximal covering network improvement problem, Socio-Economic Planning Sciences 43 (2) (2009) 102–110. [14] K.M. Curtin, K. Hayslett, F. Qiu, Determining optimal police patrol areas with maximal covering and backup covering location models, Networks and Spatial Economics 10 (2007) 125–145. [15] A. Basar, B. Catay, T. Unluyurt, Amulti-period double coverage approach for locating the emergency medical service stations in Istanbul, Journal of the Operational Research Society 62 (2011) 627–637. [16] V. Batanovic, D. Petrovic, R. Petrovic, Fuzzy logic based algorithms for maximum covering location problems, Information Sciences 179 (1–2) (2009) 120–129.

[17] O. Berman, J. Wang, The minmax regret gradual covering location problem on a network with incomplete information of demand weights, European Journal of Operational Research 208 (3) (2011) 233–238. [18] S. Davari, M.H. Fazel Zarandi, A. Hemmati, Maximal Covering Location Problem (MCLP) with fuzzy travel times, Expert Systems with Applications 38 (12) (2011) 14535–14541. [19] J. Peng, B. Liu, Parallel machine scheduling models with fuzzy processing times, Information Sciences 166 (1–4) (2004) 49–66. [20] R. Zhao, B. Liu, Standby redundancy optimization problems with fuzzy lifetimes, Computers & Industrial Engineering 49 (2) (2005) 318–338. [21] Y. Zheng, B. Liu, Fuzzy vehicle routing model with credibility measure and its hybrid intelligent algorithm, Applied Mathematics and Computation 176 (2) (2006) 673–683. [22] L. Liu, Y. Li, The fuzzy quadratic assignment problem with penalty: new models and genetic algorithm, Applied Mathematics and Computation 174 (2) (2006) 1229–1244. [23] X. Huang, Fuzzy chance-constrained portfolio selection, Applied Mathematics and Computation 177 (2) (2006) 500–507. [24] M.H. Fazel Zarandi, A. Hemmati, S. Davari, The multi-depot capacitated location-routing problem with fuzzy travel times, Expert Systems with Applications 38 (8) (2011) 10075–10084. [25] H. Ke, B. Liu, Fuzzy project scheduling problem and its hybrid intelligent algorithm, Applied Mathematical Modelling 34 (2) (2010) 301–308. [26] M. Wen, K. Iwamura, Fuzzy facility location-allocation problem under the Hurwicz criterion, European Journal of Operational Research 184 (2) (2008) 627–635. [27] M.H. Fazel Zarandi, et al., Capacitated location-routing problem with time windows under uncertainty. Knowledge-Based Systems, 2013. [28] L.A. Zadeh, Fuzzy sets, Information and Control 8 (1965) 338–353. [29] B. Liu, Y.K. Liu, Expected value of fuzzy variable and fuzzy expected value models, IEEE Transactions on Fuzzy Systems 10 (2002) 445–450. [30] B. Liu, Theory and Practice of Uncertain Programming, vol. 3rd, 2009 . [31] L. Hurwicz, Optimality criteria for decision making under ignorance, in: Cowles Commission Discussion Paper, vol. 370, 1951. [32] T.-Y. Chen, T.-M. Chi, On the improvements of the particle swarm optimization algorithm, Advances in Engineering Software 41 (2) (2010) 229–239. [33] M. Yazdani, M. Amiri, M. Zandieh, Flexible job-shop scheduling with parallel variable neighborhood search algorithm, Expert Systems with Applications 37 (1) (2010) 678–687. [34] K. Fleszar, I.H. Osman, K.S. Hindi, A variable neighbourhood search algorithm for the open vehicle routing problem, European Journal of Operational Research 195 (3) (2009) 803–809. [35] S. Mitrovic´-Minic´, A.P. Punnen, Local search intensified: very large-scale variable neighborhood search for the multi-resource generalized assignment problem, Discrete Optimization 6 (4) (2009) 370–377. [36] V.C. Hemmelmayr, K.F. Doerner, R.F. Hartl, A variable neighborhood search heuristic for periodic routing problems, European Journal of Operational Research 195 (3) (2009) 791–802.