The capacitated maximal covering location problem with heterogeneous facilities and vehicles and different setup costs: An effective heuristic approach

Article history: Received June 16 2020 Received in Revised Format August 12 2020 Accepted September 5 2020 Available online September, 5 2020 In this research, a maximal covering location problem (MCLP) with real-world constraints such as multiple types of facilities and vehicles with different setup costs is taken into account. An original mixed integer linear programming (MILP) model is constructed in order to find the optimal solution. Since the problem at hand is shown to be NP-hard, a constructive heuristic method and a meta-heuristic approach based on genetic algorithm (GA) are developed to solve the problem. To find the most effective solution technique, a set of problems of different sizes is randomly generated and solved by the proposed solution methods. Computational results demonstrate that the heuristic method is capable of producing optimal or near-optimal solutions in a rational execution time. © 2021 by the authors; licensee Growing Science, Canada


Introduction
Finding facility location is a long time frame strategic decision for companies. Therefore, it is considered as a decisive constituent in companies' strategic planning. Covering problem (CP) is one of most interesting and most appealing issues among the facility location problems due to its successful applications in different settings (Berman et al., 2016;Zhang et al., 2017). The CP is composed of two main classifications: set-covering problem (SCP) and maximal covering location problem (MCLP) (Schilling et al., 1993). The objective of the SCP, first introduced by Toregas et al. (1971), is to find the minimum number of facilities so that all demand points are satisfied within a standard time or distance. In a real-world application, due to some inadequate assigned resources such as budget, the SCP is not able to cover all demand points (Zanjirani Farahani et al., 2012). Thus, Church & Velle (1974) proposed an MCLP with the objective of maximizing the total covered demand by positioning a determined number of facilities. MCLP has many applications in various domains such as determining the location of emergency warning sirens (Current & O'Kelly, 1992), bank branches (Alexandris & Giannikos, 2010), the design of police patrol areas, the location of aeromedical and ground ambulances (Erdemir et al., 2010), Wi-Fi equipment (Lee & Murray, 2010), and the design of police patrol areas (Curtin et al., 2010) alongside RFID network design asset tracking in healthcare (Oztekin et al., 2010). Among numerous types of MCLP models that have been presented so far, a rudimentary assumption is that the facilities being located have unlimited capacities (Yin & Mu, 2012). In spite of this fact that the above assumption has many applications in real-world situations, there definitely exist situations in which this assumption severely restricts the application of the covering models (Current & Storbeck, 1988). Current & Storbeck (1988) first introduced a capacitated version of MCLP. Then, they demonstrated the theoretical links between the model and the capacitated plant location problem, the capacitated p-median problem and, the generalized assignment problem. Pirkul & Schilling (1989) developed a mathematical model for the capacitated MCLP (CMCLP) with workload capacities on facilities and the allocation of multiple levels of backup or prioritized service for all demand points. They also proposed a solution technique based on Lagrangian relaxation to solve the problem. Haghani (1996) constructed two mathematical models with two different objective functions for the CMCLP. The objective function of his first model was only set to maximize the total covered demand under capacity constraints, whereas the objective function of his second proposed model was set to maximize the weighted coverage and minimize the average distance from the uncovered demands to the located facilities. He also developed two heuristic approaches based on the Greedy adding algorithm and Lagrangian relaxation to solve the problems. Griffin et al. (2008) approached the CMCLP under budget and capacity constraints where siddeferent costs could be fixed. Besides, they provided a mathematical model to maximize the weighted demand coverage of the needy population. Shariff et al. (2010) examined three different set of data for the CMCLP. They first applied CPLEX to solve the optimization problem at hand. As the CPLEX provided results that contravene the capacity constraints when the constraints were tight, a different meta-heuristic solution method based on Genetic Algorithm (GA) was then applied to solve the optimization problem. Yin and Mu (2012) addressed a modular CMCLP and developed a mathematical model which allows several possible capacity levels for the facility at each potential site. Allocations of the demands beyond the service covering standard were also incorporated into the mathematical model to optimally locate emergency vehicles. Salari (2014) presented a generalization of the MCLP. In this generalization, any given facility had the ability to offer different service levels in accordance with the allocated budget to the facility. There, a given facility was basically able to offer different service levels according to the budget allocated. ElKady and Abdelsalam (2016) investigted a CMCLP and developed a modified Particle Swarm Optimization (PSO) algorithm to solve the problem. Recently, Bagherinejad and Shoeib (2018) combined the CMCLP and dynamic MCLP with one another, in which a dynamic capacity constraint was considered for the facilities. They developed two meta-heuristic algorithms based on GA and Bee Algorithm to solve the proposed problem.
A fleet of vehicles used in a real-world industrial application such as goods distribution company is rarely homogeneous. This is due to the fact that the lifetime of a vehicle in the fleet is often long, during which many technological and market changes happen. In addition, companies usually desire to maintain a wide variety of vehicle types in terms of their capacity and flexibility in order to cope with operational constraints (Hoff et al., 2010). In this study, a heterogeneous fleet of vehicles with different carrying capacities is taken into account. On the other hand, in practical application, the available budget to open different facilities and the available space at each potential facility site is limited. Hence, we consider budget and space constraints as well in our proposed research problem. As the MCLP has been proved to be NP-hard (Megiddo et al., 1983), the CMCLP which is an expansion to the MCLP is also NP-hard. Due to this NP-Hardness, exact techniques are not appropriate in solving large-size optimization problems in reasonable computational times. Therefore, a heuristic approach is proposed to solve the problem at hand. It is worth mentioning that to the best of the authors' knowledge, there is no research study in the body of the literature that considers heterogeneous vehicles alongside budget and space constraints at the same time for the CMCLP. As a result, to better evaluate the results obtained using the heuristic method, a GA is developed to solve the problem as well.
In summary, we approach a CMCLP with the following assumptions in this study: a) A heterogeneous vehicles fleet with different shipping capacities is considered. b) Limited budget to open different facilities is also considered. c) Confined space is taken into account at any potential facility site. d) The capacity of each facility is different.
The objective is to determine the numbers and the types of required vehicles in addition to locating new facilities so that the covered demand is maximized.
The rest of this paper is organized as follows. In Section 2, a mathematical model is developed for the proposed research problem. In Section 3, the heuristic procedure and GA are introduced to heuristically solve the problem. The specifications of the test problems used to compare the performance of the proposed algorithms along with parameter settings are examined in Section 4. In Section 5 numerical results and comparison analysis involving the proposed algorithms are presented. Conclusions and future research areas are discussed in Section 6.

The mathematical model
In this section, an original MILP model is developed for the proposed research problem. The Indices, parameters, decision variables, and the whole mathematical model are as follows.

Indices:
Index of a demand node The total available space at site ; The maximum available budget B The coverage distance or the coverage radius R if the distance from candidate site to an existing demand node is not greater 1

Decision variables:
f demand node is covered by a facility located at site , The mathematical formulation of the problem is The objective, as presented by Eq.
(1), is to maximize the total covered demand. Constraint (2) is introduced into the model to ensure the budget allocated to locating facilities and vehicles does not exceed the total available budget. Constraint (3) guarantees that if there are no facility type at site or demand node is not covered by site , the corresponding decision variable is forced to be equal to 0. Constraint (4) states that each demand node is covered by at most one facility. According to Constraint (5), at most one facility type can be located at each site . Constraint (6) guarantees that the total capacity of the vehicles is greater than or equal to the total covered demand. Constraint (7) ensures that the total capacity of the vehicles cannot be greater than the total capacity of the facility type k at site j. Constraint (8) indicates that the sum of the occupied space by the facilities and the vehicles does not exceed the total available space ( ) at each site . The decision variables are described in Eqs. (9)-(11).

The solution approaches
As described in Section 1, the CMCLP belongs to the class of NP-hard problems. Hence, exact approaches are computationally inefficient, notably for large-sized problems. Consequently, heuristic and meta-heuristic solution procedures are required to come up with rough optimal solutions when the size of the problem increases.

The heuristic method
In this section, the steps involved in the proposed heuristic approach to solve the problem are discussed. The peculiarities of the proposed heuristic are as follows: : is the maximum number of iterations : is the iteration counter : is the budget at th iteration, and : is the value of the objective function. Z

In addition,
For it=1 to IT, set j=1, k=1, and BestFunction=0, where BestFunction is the best value of the objective function, start taking the following steps: . if > , go to Step 3. Otherwise, go to Step 11.
Step 4. Otherwise, set Z 0, where Z is the value of the objective function of the facility type at site and go to Step 9.
, then calculate Z and by the formulae mentioned in Step 5 and go to Step 7. Otherwise, go to Step 6.
At the end of this step, the remaining demand nodes ( ) are the covered demands.
Step 8 . Otherwise, go to Step 9. Step 10 . if the stopping criterion (the maximum number of iterations) is met, the process is terminated and the value of the final objective function ( ) is stored. Otherwise, continue t Z Step 11 he process. The flowchart of the proposed heuristic is shown in Fig. 1.   Fig. 1. The flowchart of the heuristic approach

Genetic algorithm (GA)
The GA, initially introduced by Holland (1975), is a population-based stochastic algorithm that evolves through a number of initial solutions called chromosomes. At each iteration, a series of basic operators suvh as selection, crossover and mutation are performed on the population with the aim of performing both exploration and exploitation (Michalewicz, 1996;Hansheng & Lishan, 1999). In this research, due to the remarkable performance of GA in solving the MCLP (Jia et al, 2007;Zarandi et al., 2011;Shariff et al., 2012;Atta et al., 2018), it has been considered for solving the research problem under examination. The following sections provide the details on the adopted algorithm.

Chromosome representation, initial population, and fitness evaluation
A chromosome is defined as an encoded solution while the decision variables within a solution (chromosome) are genes (Talbi, 2009). A solution representation in the proposed GA includes ( ) n z t m    components and three parts that have been combined. The first part represents an n m  matrix of binary elements, in which each row and column displays a demand node and a facility site, respectively. If demand node i is covered by a facility at the site j , then the corresponding element of this matrix takes a value of one; otherwise, it takes a value of zero. The second part consists of a z m  matrix of binary elements, in which each row and column represents a facility type and a facility site, respectively. If facility type k is located at the facility site j , then the corresponding element in this matrix is set to one. Otherwise, it is set to zero. The last part shows a t m  matrix with the elements randomly generated between 0 and N . Each row and column of this matrix represents a vehicle type and a facility site, respectively. In addition, each element of this matrix displays the number of type-k vehicles assigned to a potential site j . Fig. 2 illustrates the general structure of the proposed chromosome. To initialize, all elements (genes) of a chromosome take a value of zero. Then, the following steps are taken to assign values to the elements.
Step 1. Calculate the total demand at each site j by Step 2. Select one j and one k on a random basis, then, set 1. jk y  Step 3 x randomly takes a value between 0 and 1.
Step 6. Select one  in a random manner and then randomly select a value as j n  in the range between 0 to N , where N is determined as follows:

The length of a chromosome
Due to the possibility of generating an infeasible solution, a penalty function is considered in the reformed objective function formulated in Eq. (12).
where M is a big constant number, m is the mean violations happened for each constraint set, and l g is defined as a degree of infeasibility of constraint set l which is set as the difference between the left-hand side and the right-hand side of the constraint set . l

Selection
The objective of the selection procedure is to reproduce more copies of individuals with higher fitness values (Pham & Karaboga, 2000). The selection procedure has a considerable impact on driving the population towards better solutions. The selection strategy determines which individuals are selected for reproduction (Talbi, 2009). In this study, the tournament selection method is used where a competition takes place among q randomly selected members to choose the one with the highest fitness value.

Crossover
This operator is used to produce two new individuals (offspring) from two existing individuals (parents) chosen from the current population by the selection operator (Pham & Karaboga, 2000). In the proposed GA, one cut-point on each side of the chromosome (a parent) is randomly selected and two new chromosomes are produced by swapping the segments of the parents. Selecting the crossover cut-points is based on a uniform random distribution. An example is illustrated in Fig. 3 to clarify how the crossover operator functions. The part of the population on which the crossover operator is executed in an iteration depends on the crossover probability ( c p ). Fixed part (yellow numbers)

Mutation
The goal of the mutation operation is to maintain the population diversity, thus avoiding any premature convergence. It makes a new offspring up from a single parent and assists the algorithm to avoid to be trapped in local optima (Beldar & Costa, 2018). During the mutation operation, a certain percentage of the genes of a chromosome ( u m ) is selected and the mutation operator is performed. It is worth mentioning that the developed chromosome is formed of two binary sections and one discrete section. Therefore, two different strategies are required for doing mutation. The Bit flip mutation (Sivanandam & Deepa, 2008) is utilized for the first two sections where it involves changing 0 to 1 and 1 to 0. Fig. 4 illustrates an example of the mutation-flipping for better clarification. A chromosome is considered as a parent at random. For a 0 in the selected chromosome, the corresponding bit is flipped to 1 and a child chromosome is produced. A related point to consider is that if ij a of the selected gene (element) from the first section of the chromosome is equal to zero, then the selected gene remains unchanged. On the other hand, in the discrete section, the value of the selected gene is replaced with a new value in the range from 0 to N . An example is presented in Fig.  5 to illuminate the point. The portion of the population in which mutation procedure is applied depends on the mutation probability ( m P ).

Creating new population
All the newly generated chromosomes yielded by the crossover and the mutation operators are merged with the current population. Then, the best P chromosomes are selected in terms of the objective function value, where P is equal to the size of the population.

Parent
New offspring Parent New offspring

Termination criterion
While various types of stopping criteria are discussed in Sivanandam & Deepa (2008), in this research, the aforementioned process stops if there is no improvement to the population's best fitness for a specified number of iterations (four iterations).

Computational experiments
In order to compare the performances of the proposed solution approaches, a set of test problems are randomly generated. The whole set of the test problems are solved by both heuristic and GA, based on which the performances of these algorithms are compared. The test problems specifications and the parameter tuning approaches are discussed in the subsequent sub-sections.

Test problems
The procedure adopted for data generation is derived from Karasakal & Karasakal (2004). The demand nodes are generated using the discrete uniform distribution in the ranges of [0, 50], and [0, 100]. Similarly, the potential facility sites are produced generated using the discrete uniform distribution in the ranges of [5, 45], and [10, 90]. It is assumed that the facility sites can fully cover the demand nodes within the range of 20. The Euclidean distance is utilized to calculate the distances between the demand nodes and the potential facility sites. Data generation for the remaining variables is provided in Table 1. In this study, the total number of investigated test problems is 30.

Parameters calibration
Both the effectiveness and the efficiency of meta-heuristic algorithms can be enhanced by providing appropriate values to their control parameters (Beldar & Costa, 2018). Therefore, in this section, the behavior of the proposed algorithm is evaluated while changing the control parameters' values. The Taguchi method is used here to adjust the parameters because it can examine a large number of decision variables with a few experiments (Phadke, 1989). In the Taguchi method, the factors are categorized into two groups of noise and controllable factors. The noise factors are those on which one does not have direct control over. Taking into account the fact that it is often impossible to remove the noise factors, the Taguchi method seeks to minimize the noise effect and to determine the optimal level of controllable factors on the basis of the sustainability concept (Al-Aomar, 2006). The Taguchi method both calculates the relative importance of each factor in terms of their major impact on the response variable and determines the optimal values of the factors. The main reason the Taguchi's method is considered as a sustainable design is that it transmits the values of the response variable within a ratio called signal-to-noise (S/N). The term "signal" indicates the desirable value (mean/average response variable) and the term "noise" indicates the undesirable value (standard deviation). Hence, the S/N ratio represents the amount of variation in the response variable. Taguchi categorizes all objective functions into three groups: of "Smaller-the-better," "Larger-the-better," and "Nominal-the-best". For the proposed model, the larger-the-better is utilized and the corresponding S/N ratio is calculated as follows (Phadke, 1989): where i y is the observed response value and n is the number of replications. For GA, the most influencing parameters that need to be adjusted are c P , m P and u m . Each control parameter is assumed to vary at three levels, conforming to low, medium and high levels, as depicted in Table 2. From the standard table of orthogonal arrays, L9 is selected as the most appropriate orthogonal array. Each test problem is iterated seven times for each scenario of the orthogonal array. Table 3 shows the optimal values of the control parameters obtained by the Taguchi-based calibration method. Other control parameters, including the population size, the maximum number of iterations, and the tournament size are obtained using a trial-and-error approach involving an extended number of instances. Thus, the values of the maximum number of iterations, the tournament size and the number of population, are 10, 30 and at least 100, respectively.

Numerical results and comparison analysis
The proposed solution approaches were coded in MATLAB 2012a. The whole set of experiments was performed on a PC with 3 GHz CPU and 1 GB RAM. In order to validate the outcomes obtained by the two proposed solution methods, the Relative Percentage Deviation (RPD) performance indicator as shown in Eq. (14) is used.
Optimal Solution Algorithm Solution Optimal Solution RPD 100    Note in Eq. (14) that the global optimal solution of only 23 designed test problems can be attained using the ILOG CPLEX (version 12.2). Thus, the RPD values based on the remaining 7 test problems are obtained using their local optima achieved in three hours of CPLEX execution. The elapsed time, the objective function value and the RPD value of each test problem are provided in Table 4. As it can be comprehended from the results shown in Table 4, in the last row, the total averages (Average) underline the efficacy of the heuristic approach over GA in terms of both the solution quality and the required computational time. In addition, when the size of the problem increases, the weakness of GA becomes more evident. MINITAB 16 commercial package was utilized to complete the statistical analysis. The Box plot diagram at a 95 percent confidence level depicted in Fig. 6 verifies that there is a significant statistical difference between the solution approaches and that the heuristic methodology outperforms the GA.