Dynamic capacitated maximal covering location problem by considering dynamic capacity

Article history: Received January 15 2017 Received in Revised Format April 1 2017 Accepted May 28 2017 Available online May 29 2017 Capacitated maximal covering location problems (MCLP) have considered capacity constraint of facilities but these models have been studied in only one direction. In this paper, capacitated MCLP and dynamic MCLP are integrated with each other and dynamic capacity constraint is considered for facilities. Since MCLP is NP-hard and commercial software packages are unable to solve such problems in a rational time, Genetic algorithm (GA) and bee algorithm are proposed to solve this problem. In order to achieve better performance, these algorithms are tuned by Taguchi method. Sample problems are generated randomly. Results show that GA provides better solutions than bee algorithm in a shorter amount of time. © 2018 Growing Science Ltd. All rights reserved


Introduction
Facility location problem is a special class of optimization problems whose primary goal is to locate a limited number of facilities that satisfy particular constraints (Máximo et al., 2017).Facility location problems have been studied widely during recent years due to their extensive application in real situations (Correia & Captivo, 2006).Location problems can be defined according to two factors; space (planning area) and time (time of location).Space and time issues have been taken into account in static facility location problems and dynamic facility location problems respectively (Boloori Arabani & Zanjirani Farahani, 2012).Boloori Arabani and Zanjirani Farahani (2012) classified different types of static and dynamic location problems that have been studied by the literature review.They studied multi-period facility location problem as a type of dynamic location problems.Static location problem considers only one period.If a time horizon is considered for more than one period, the location problem becomes dynamic (Canel et al., 2001).By considering a time horizon with more than one period, determining the appropriate time for facility location, specifying the best locations and better prediction of favorable and unfavorable fluctuations of demand in time horizon can be achieved; whereas single period models do not have these characteristics (Miller et al., 2006) .
Dynamic models can be classified into two categories: explicitly dynamic models and implicitly dynamic models.In explicitly dynamic models, in order to respond to changes in parameters over time, facilities are closed and opened in pre-specified times and locations.In implicitly dynamic models, all facilities are to be open in the beginning of time horizon and remain open throughout the time horizon.These models are considered to be dynamic because they try to consider changes in parameters such as demand changes over time (Current et al., 1998).
Due to largely capital outlaid, facility location problems are frequently long-term in nature.Facilities such as schools, hospitals and dams operate for decades (Current et al., 1998).While, some facilities such as buses may move around in order to meet the demands of the population (Datta, 2012).Decision makers should select locations and consider time of relocation of facilities which are able to response demand fluctuations over the time (Daskin et al., 1992).Therefore, in order to control probable fluctuations in the future as well as parameters fluctuations a dynamic model seems to be necessary (Boloori Arabani & Zanjirani Farahani, 2012).Theoretically, opening/closing of facilities could impose no cost (Hormozi & Khumawala, 1996).One of the objectives in facility location problem is to minimize the total cost for assigning facilities to satisfy the demand nodes (Jahantigh & Malmir, 2016).
One of the traditional location problems is covering location problem.Covering location problem seeks for a solution to cover a subset of customers by considering one or more objective (Davari et al., 2013).Although covering models are not new, due to their application in real cases especially for emergency service facilities, they have always been attractive topics for researchers (Fallah et al., 2009).Static models can be transformed to their equivalent dynamic models.In these models instead of single period, T periods are considered.Therefore, maximum covering location problem could be studeid as a multiperiod and dynamic problem (Boloori Arabani & Zanjirani Farahani, 2012).
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. Section 3 is dedicated to the definition of the problem.The proposed solution algorithm is presented in Section 4 and numerical examples and parameter setting appear in Section 5.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 6. Schilling et al. (1993) classified covering location problems in two categories named maximal covering location problem (MCLP) and set covering location problem (SCLP).In covering problems, a demand is said to be covered if at least one facility is located within a predefined distance of it.This predefined distance is often called coverage radius.The objective of SCLP is to cover all demand with the minimum number of facilities.Covering location problem was introduced for the first time by Hakimi (1965).The objective of that model was to determine the minimum number of polices that was necessary to cover demand nodes on a network of highways.SCLP was introduced formally by Toregas et al. (1971) and was extended slightly by Berlin and Liebman (1974).MCLP was introduced for the first time by Church and Revelle (1974).The objective of the MCLP is to locate a fixed number of facilities in such a way that the total covered demand is maximized.

Literature review
The main assumption of covering location problems is that the facilities are uncapacitated (Salari, 2013).But, practically this assumption is not always valid (Pirkul & Schilling, 1991) and usually limits the applications of covering models (Current & Storbeck, 1988).Most service facilities are capacitated (Murray & Gerrard, 1998;Liao & Approach, 2008).Therefore some covering models considered a capacity constraint for facilities.Although incorporating capacity constraint in formulation of location problems is not difficult, it increases computational complexity.Therefore, most research efforts focus on improvement of solving method (Pirkul & Schilling, 1989).Current and Storbeck (1988) incorporated capacity limitation to MCLP and LSCP.There is a theoretical link between these models and capacitated plant location problem, the capacitated P-median location problem and generalized assignment problem.Small and moderately sized problems can be solved with existing solution methods.Theoretical links give insight into developing new heuristics for large sized capacitated covering problems.Pirkul and Schilling (1991) developed the model by considering workload limits on the facilities and quality of service delivered to the uncovered demand zones.Facility's workload limits the demand amount which a facility can serve.In this model, all demands are allocated to facilities regardless of whether they are in covering radius or not.The quality of service is modeled as the total distance from uncovered demand zones to the nearest facility.This model is solved by an approach based on the Lagrangian relaxation (Pirkul & Schilling, 1991).Haghani (1996) proposed a capacitated maximal covering location problem in which weighted covered demand is maximized and average distance from the uncovered demands to the located facilities is minimized.They solved the problem with two heuristic methods.The first one was based on greedy adding technique and the second one was based on Lagrangian relaxation.Correia and Captivo (2003) considered modular capacitated location problems.In this model, instead of considering only one fixed capacity level for each facility, they considered a discrete and limited set including available capacity levels.Capacity level of facility is selected from this set.This model can be applied in schools, warehouses and other public facilities.Griffin et al. (2008) proposed capacitated maximal covering location problem by considering three capacity levels for each facility.In their model, there is no composing relationship (such as that between the number of ambulances and emergency stations) between facilities' capacity levels.Yin and Mu (2012) proposed modular capacitated maximal covering location problem (MCMCLP) in two situations.In these models, it is assumed that each facility has a capacity which is related to number of vehicles assigned to that facility.Vehicles have a fixed capacity but the capacity of each facility is equal to the total capacity of vehicles assigned to that facility.In the first model, the number of vehicles is predefined but in the second model, number of vehicles as well as number of facilities is predefined.Yin and Mu (2012) stated that this is a static model and disregard dynamic factors such as daily population movement.On the other hand, since MCMCLP is NP-hard, proposing a heuristic for this problem is important (Yin & Mu, 2012).Although the papers surveyed above, considered capacitated MCLP in only one period, the concept of dynamic (multi-period) covering location problem is not new in the literature (Fazel Zarandi et al., 2013).Schilling (1980) is among the first researchers who considered dynamic maximal covering location problem.Also, other researchers have taken into account this problem.Fazel Zarandi et al. (2013) considered large scale dynamic MCLP.Dell'Olmo et al. (2014) proposed a multi period MCLP for the optimal location of intersection safety cameras on an urban traffic network.Vatsa and Jayaswal (2016) present a new formulation for a multi-period MCLP with server uncertainty, motivated by its relevance with respect to primary health centers.
In this paper, by integrating modular capacitated maximal covering location problem and multi-period maximal covering location problem, a developed model is proposed.

Problem definition
In the proposed model, a time horizon consisted of T periods is considered.The objective of this model is to find optimal location of q facilities in a time horizon in such a way that with locating at most pt vehicles in period t, the maximum covering is achieved in the whole time horizon.This model can be applied in location facilities such as ambulance bases and vehicles such as ambulances.In this model, it is assumed that each vehicle has a fixed capacity (Yin & Mu, 2012) equal to maximum amount of demand that it can serve in each period.Capacity of each facility in each period is related to the number of vehicles stationed in that facility (Yin & Mu, 2012).So, facilities' capacity is changing periodically and we call it as dynamic capacity.For example, if capacity of ambulance is C and is the number of ambulances located in ambulance base in location j in period t, the capacity of that ambulance base will be CZ .
Another assumption is that potential locations are identical in all periods.In each potential location, only one facility can be located and if a facility were located in location j in period t, facility would serve in this location until the end of the time horizon.In other words, the cost/penalty of closing facilities is so high that prevents closing facilities (for instance buildings such as hospital).Since the importance of costs in public sectors is inconsequential compared to provided services (Fazel Zarandi et al., 2013), it is assumed that opening and closing of vehicles and relocation of them has no cost.Therefore, vehicles are closed at the end of each period and are relocated again in the next period (if they were available).
Since some vehicle may become unavailable in each period because of being out of use, etc. 0 or some new facilities are added ( 0), the number of vehicles are not considered to be identical in all periods.In especial situation, the number of facilities is identical in each period ⋯ .In the proposed model, each facility in each period is as a potential location for stationing of vehicles.If there is no facility in a period, there would be no potential location for stationing of vehicles.Therefore, to maintain the feasibility of the problem, the constraint on the number of vehicles in each period is considered as the maximum number of vehicles.The maximum number of vehicles is given in the beginning of each period and no limitation on the number of vehicles which can be stationed in a facility is considered.In this model, a constraint on the number of facilities is considered in the whole time horizon.If a constraint on the minimum number of new facilities which can be located in period t is not imposed, the minimum number of new facilities in period t will be zero ( 0 , if a constraint on the maximum number of new facilities which can be located in period t is not imposed, the maximum number of new facilities in period t will be q ). q is the total number of facilities in the time horizon.In some periods, a constraint on the minimum or maximum number of new facilities located might be imposed in each period.In this situation, the decision maker determines the minimum number of facilities in each period in such a way that sum of minimum number of facilities in the time horizon would not exceed the total number of facilities and considers the maximum number of facilities in each period more than the minimum number of facilities in each period.It is assumed that the minimum and maximum number of new facilities in each period is certain and predefined.It is assumed that at the end of each time horizon, all facilities and vehicles are closed.So, in the beginning of each time horizon, no facilities are located in potential locations 0 .Hereby, the proposed dynamic capacitated MCLP is presented.First, problem parameters and variables are defined.

Sets and parameters
i, I: The index and set of demand nodes j, J: The index and set of eligible facility locations t, T: The index and set of time periods : The population or demand at node i in period t d: The Euclidean distance from demand node i to facility at j S: The distance (or time) standard within which coverage is desired N={j|d ≤ S}: the set of nodes that are within a distance less than S from node i : Maximum number of vehicles in period t q: The number of facilities to be located throughout the time horizon : Minimum number of new facility in period t (∑ ) : Maximum number of new facility in period t ( ) c: capacity of each vehicle Variables : A binary variable which equals one if a facility is sited at location j in period t. : A binary variable which equals one if node i in period t is covered by one or more facilities stationed within a distance of S. : An integer variable which indicates the number of vehicles which are located in period t and site j.
Then, the proposed model will be as follows: (1) (2) (3) ∀ , (10) The objective function (1) maximizes the overall covered demand.Constraint (2) shows that q facilities are to be established in T periods.Constraint (3) specifies the maximum number of vehicles to be located in each period.Constraint (4) assures that if minimum or maximum number of new facilities in period 1 be predefined, the number of new facilities in period 1 (t=1) will be in the predetermined interval.
Otherwise, it will be between zero and q (the constraint on the minimum and maximum number of new facilities in t=1).Constraint (5) specifies that if minimum or maximum number of new facilities for t>1 be predefined, the number of new facilities will be in the predetermined interval.Otherwise, it will be between zero and the number of available facilities (number of facilities is not located until period t).Constraint (6) specifies that the demand point i in period t is covered, if it does not exceeds the total capacity of facilities which can cover this demand point.Constraint (7) ensures that the total covered demand in each period cannot exceed total capacity of facilities in that period.Constraint (8) specifies that if a facility is located in period t, it will remain open until the end of the time horizon.Constraint (9) specifies that if a facility is located in period t in site j, the number of vehicles assigned to this facility cannot exceed pt (the maximum number vehicles in period t).In other words, in each period vehicles can be stationed in site j if a facility is located in that site and this cannot be more than the number of vehicles predefined for each period.Constraint (10) specifies that decision variables and are binary.Constraints (11) restrict the integer decision variable zjt, which ranges from zero to pt.

Genetic algorithm (GA)
4.1.1.Review of GA GA was first proposed by Holland (1975) as an evolutionary algorithm.It is based on Darwinian evolution: good traits survive and mix to form new while the bad traits are eliminated from the population (Zanjirani Farahani & Hekmatfar, 2009).Beasley and Chu (1996) seem to be the first to apply GA for covering model.The Simple GA is as follows: 1. Generate an initial population mostly in a random way.The rest of this section is devoted to elaboration of the proposed GA.

Encoding scheme
An appropriate chromosome representation must be defined for a GA.Encoding very depends on the problem.Most previously adopted representations, such as the bit string, are linear or one-dimensional.Some real problems are naturally suitable for two-dimensional representation.In this paper, GA with multi chromosome representation is applied and a separate chromosome is considered for each variable (Matthias et al., 2013).Therefore, three chromosomes are defined which are two-dimensional.The first chromosome is binary matrix which has I rows and T columns.Each bit in this matrix represents the status (covered/uncovered) of the node i in period t.The second chromosome is a binary matrix which has J rows and T column and each bit in this matrix represents the status (located/unlocated) of facility at location j in period t.In other words, one value in the ith bit means that there is a facility at location j in period t.The third chromosome is a matrix consisted of integer number which has J rows and T column and each bit in this matrix represents the number of vehicles at location j in period t.Each bit of third chromosome is an integer number between zero and pt in such a way that total number of vehicles in each period is not more than pt.Therefore, an upper and lower bound are considered for the number of vehicles in each period.At first, an initial population is generated randomly for each chromosome.

Selection
Selection is the stage of GA in which chromosomes are selected from the population to be parented for the next generation.There are many methods in selecting the best chromosomes.In this paper, the Roulette Wheel Selection (RWS) has been used as the selection method.Each individuals of the population is allocated a section of an imaginary roulette wheel.The size of each sections is directly proportional to their fitness values, such that the fittest individual has the biggest slice of the wheel and the weakest individual has the smallest.The wheel is then spun and the individual associated with the winning section is selected.

Crossover
Crossover and mutation are two basic genetic operators used to make new off springs.Type and implementation of operators depend on encoding and also on a problem.In this paper, crossover in first and second chromosomes is performed using one of these two methods: one-point crossover and twopoint crossover.But in the third chromosome, the arithmetic crossover is applied.Fig. 1 shows different types of crossover.

One-point crossover
One-point crossover randomly select a crossover point and then copy everything before this point from the first parent and then everything after the crossover point copy from the second parent.Here chromosomes are two-dimensional and crossover point can be a bit or column.In one-point crossover (column), third column in each parent is selected randomly and exchanged.In one-point crossover (bit), two bits in each parent selected randomly and exchanged (Fig. 1).

Two-point crossover
Two-point crossover randomly selects two crossover points within a chromosome then interchanges the two parent chromosomes between these points to produce two new offspring.In Fig. 1, in the case of two point crossover first and forth columns in each parent are selected and they are swapped.

Arithmetic crossover
Arithmetic crossover is a linear combination of two chromosomes as follows: . 1 (20) where C is a parent, is an offspring, and α is random matrix (between 0 and 1).This operation is performed for each bit (Köksoy & Yalcinoz, 2008).In Fig. 1, random matrix is given.The number of bits are selected randomly.According to the Eq. ( 20) and Eq. ( 21) arithmetic crossover is operated.

Mutation
Mutation is the genetic operator that randomly changes one or more of the chromosome's gene.The purpose of the mutation operator is to prevent the genetic population from converging to a local minimum and to introduce to the population of new possible solutions.The mutation is carried out according to the mutation probability.Mutation rate is usually set to a very low level.However, different references have found that a higher mutation rate is necessary when the GA has converged (Jaramillo et al., 2002).In this paper, crossover in first and second chromosomes is performed using one of these methods: swap, binary and reversion mutations.In third chromosome, the mutation method is only arithmetic mutation.Fig. 2 shows different types of mutation.

Binary mutation
In Binary mutation, some bits are selected randomly and their values are changed from zero to one and one to zero.In Fig. 2, in case of binary mutation, colored bits are selected randomly and binary mutated.

Swap mutation
In Swap mutation, two bits are selected randomly and are exchanged with each other.In Fig. 2, in the case of swap mutation, colored bits are selected randomly and swapped.

Reversion mutation
In reversion mutation, two bits are selected randomly and bits between them are reversed.In Fig. 2, in the case of reversion mutation, colored bits are selected randomly and reversed.

Arithmetic mutation
The definition of arithmetic mutation is as follows: where is upper bound, C is a parent, is an offspring, and α is random matrix (between 0 and ).At each iteration, offspring as many as the population size are created.In order to produce new generation, these offspring replace less fit individuals in the existing population.In Fig. 2, in the case of arithmetic mutation, upper bound and random matrix is given ( 10).Certain number of bits are selected randomly.According to the Eq. ( 22) arithmetic mutation is operated.

Termination criteria
The algorithm will iterate until the maximum number of iterations is attained.

Review of Bee algorithm
Bee algorithm is one of the swarm-based algorithms which imitate the food foraging behavior of honeybee swarms (Tsai, 2014).The basic version of the algorithm performs a kind of local search combined with random search and it can be used for combinatorial optimization and functional optimization (Pham et al., 2006).The foraging process begins in a colony by scout bees being sent to search for promising flower patches.At first, these bees spread out in flower patches randomly.Each selected food source represents a feasible solution of the problem.The nectar amount in the food source represents the quality of the problems' solution.When the foraging process is finished, based on the nectar amount, food sources are ranked into three categories: elite, good (average) and bad (unselected) sites.Each scout bee performs a dancing known as the waggle dance above a certain quality threshold deposit their nectar.This way, it transfers information of that region (in comparison to the hive), its distance from the hive and its quality rating to other bees.This information helps the colony to send follower bees to flower patches.Most follower bees go to region where are more promising to find nectar.
In other words, more follower bees are assigned to elite sites than those of good sites.Assigning follower bees is as generating a neighbor for a solution.In each iteration of the algorithm for each elite site, a certain number of neighbors are generated.In this paper, the method for generating neighbors is similar to mutation methods in GA.The fitness of these neighbors is calculated and for each site the best neighbor is recognized.The quality of the best neighbor is compared to that of current bee.If the quality of the best neighbor bee is better than that of current bee, it will be replaced.Otherwise, current bee stays in its site with no movement.This process is performed on the good site.The difference is that the number of neighbors generated for good site is less than those of elite sites.Bad sites are abandoned and scout bees assigned to bad sites are assigned to other sites randomly.Therefore, the new population is generated.This process continues utile the criteria is met.In this algorithm, chromosomes defined in GA are considered as bee and Mutation operator is performed to generate neighbors.Bee algorithm is as follow: 1.The bees of initial population are randomly generated (n: number of initial bees).
2. The fitness are calculated.
3. Certain number of best bees (m) are selected to finding neighborhood.Among selected bees (m), a certain number of them (e) are considered as elite bees.The rest of them (m-e) are good bees.4. By mutation method, the neighborhood are generated for elite and good bees and their fitness are calculated.The number of neighborhood are generated for each elite bees are more than good bees.5.Among neighborhoods of each bees, the best bee is selected.These bees are transferred to the next generation.6. Good bees (n-m) are used to random search and their fitness are calculated.These bees also transmitted to the next generation.7. If a stop condition is met, the algorithm stops.Otherwise, go to step 3.

Termination criteria
The algorithm will iterate until the maximum number of iterations is attained.

Test problems
To generate test problems, a similar approach to Revelle et al. (2008) was employed.In this approach, the locations of nodes were randomly generated using a uniform distribution between 0 and 30 for both x and y coordinates.The distances between the nodes are then defined as their Euclidean distance.Populations on the demand nodes for each time period were randomly generated using a uniform distribution between 0 and 100.Fazel Zarandi et al. (2013) by considering 5 periods use this approach for generating dynamic sample problems.This paper by considering time scale uses this method for generating 30 sample problems.The minimum number of new facilities in each period is randomly generated using a uniform distribution between 0 and q in such a way that sum of minimum number of new facilities in the whole planning horizon be less than or equal to q ∑ .The maximum number of new facilities in each period is generated by random numbers larger than or equal to the minimum number of new facilities in each period .Lingo 8.0 was used to solve these problems and results were compared against those obtained using the GA and bees algorithm.

Parameter setting
Heuristic and metaheuristic algorithms are sensitive to their parameters; A small change can affect the quality of the solution.So, tuning algorithms are necessary to find a better solutions (Pasandideh et al., 2015).The most widely used method to tune the algorithms is a full factorial design (Chan et al., 2015).This method does not seem effective when the number of parameters significantly increases, since it requires arduous task to conduct the experiment.A family of matrices is used to reduce the number of experiments.In Taguchi method, we utilize the orthogonal arrays to investigate a large number of decision variables with a small number of experiments (Raju et al., 2014).
In this method, factors are classified into two groups: Controllable factors (signal) and uncontrollable factors (noise).Also, objective functions are categorized into three groups: "the smaller the better", "the larger the better" and "the nominal value is expected".The objective functions of the proposed model is "the larger the better".S/N ratio (the large the better) is calculated by Eq. ( 23).
where = observed response value and n= number of replications.According to Taguchi's design of experiments, for 4 parameters and 3 levels L9 Taguchi orthogonal array was selected (Table 1 and Table  2).For calibration of each algorithm, 6 sample problems are iterated for 5 times in each scenario.Since sample problems' dimension is not identical, so the differences between their objective functions are large and using raw data cause to wrong analyses.In other words, the dimension of the problem should be excluded form data.So after conversion of raw data to relative deviation index (RDI), the S/N ratio is calculated.Relative deviation index is calculated by Eq. ( 24): 0 1. .6 , 1. .5 , 1. .9 where ijk OF is the objective function value related to iteration j in sample problem i in scenario k. li and ui are minimum and maximum values for ith sample problem respectively.The S/N rate for scenario k can be calculated (using relative deviation index of objectives function) by Eq. ( 25).
In Taguchi method, S/N rate is considered as the first criterion.There could be no meaningful difference between different S/N levels, so, another criterion named RDIk is introduced for scenario k which is calculated by Eq. ( 23). is considered as the smaller the better.
It is time consuming to set all effective parameters in bee algorithm.Therefore, we set the most important effective parameters and other parameters have been determined by try and error.According to S/N (Fig. 3 and Fig. 5) and RDI (Fig. 4 and Fig. 6) the selected levels are colorful.

Results and discussions
In this paper, Lingo software package has been applied to find the exact solutions for some sample problems.Lingo uses branch and bound method to solve the problems.Objective bound specifies the theoretical bound of objective function.This bound specifies how much a solver can improve the objective function.The best solution cannot exceed the objective bound.Colored rows specify problems in which Lingo cannot find the optimal solution in one hour.In such cases instead of optimal value, the objective bound and the best solution found in one hour is reported (Niroomand, 2008).In such cases, metaheuristic/heuristic algorithms might find better solution than what lingo finds in one hour.In such situation, the gap will be negative (Mehdizadeh & Afrabandpei, 2012).Gap is calculated as follow: 100 According to computational results of 30 sample problems, it could be concluded that Lingo can find the optimal solution for only one thirds of the problems.In half of the sample problems, GA finds a solution  0).In other problems except one problem, GA can achieve a gap less than 1.9%.The run time of GA in the largest problems is less than 2.5 minutes.Bee algorithm finds a solution better than or equal to what lingo finds in one hour in less than half of the sample problems ( 0).In other problems except two problems, Bee algorithm can achieve a gap less than 2.9%.The run time of Bee algorithm is less than 4 minutes in the largest problems.Totally, it can be stated as following: Therefore, GA can find better solution in a shorter time.The average run time and objective function value in 5 iterations is reported in Table 3.

Conclusion and future research areas
In this paper, capacitated MCLP and dynamic MCLP were integrated to each other and dynamic capacity constraint was considered for facilities.Therefore, the MCLP has been extended to the capacitated dynamic MCLP.The developed model was solved by GA and bee algorithm and the results were compared to the exact solutions of Lingo.We have shown that while GA and bee algorithm are superior to the exact method in terms of runtime, there are negligible errors compared to the optimal solutions.GA found better solutions in a shorter amount of time than the bee algorithm.Although GA shows great performance to solve this model, one may assess the performance of other methods in finding solutions to this problem.Another opportunity for research is to add a constraint on the number of vehicles which can be located in each facility.A possible future study could be to integrate this model with gradual covering location problem.Another future research is to consider cost for each vehicle.Cost can be dynamic and changes in each period.Objective function could be maximization covered demands while the cost of vehicles is minimized.Some parameters can be fuzzy such as covering radius.Covering radius can be dynamic, too.
2. Select individuals for reproduction.3. Perform genetic operations and generate a new generation.4. Insert offspring into population and form the new population.5.If the predefined stopping criteria are met, stop the algorithm, otherwise, return to step 2.

Table 1 Table 2
Parameter and their levels in GA Parameters and their levels in Bee algorithms

Table 3
Computational results