A CONGESTED CAPACITATED LOCATION PROBLEM WITH CONTINUOUS NETWORK DEMAND

. This paper presents a multi-objective mixed-integer non-linear programming model for a congested multiple-server discrete facility location problem with uniformly distributed demands along the network edges. Regarding the capacity of each facility and the maximum waiting time threshold, the developed model aims to determine the number and locations of established facilities along with their corresponding number of assigned servers such that the traveling distance, the waiting time, the total cost, and the number of lost sales (uncovered customers) are minimized simultaneously. Also, this paper proposes modified versions of some of the existing heuristics and metaheuristic algorithms currently used to solve NP-hard location problems. Here, the memetic algorithm along with its modified version called the stochastic memetic algorithm, as well as the modified add and modified drop heuristics are used as the solution methods. Computational results and comparisons demonstrate that although the results obtained from the developed stochastic memetic algorithm are slightly better, the applied memetic algorithm could be considered as the most efficient approach in finding reasonable solutions with less required CPU times


Introduction
As one of the elements of paramount importance in the field of logistics, the facility location problem (FLP) has been scrutinized by a myriad number of researchers, especially in the context of operations research applications.According to the existing literature reviews [1][2][3], FLP could be categorized into various classes based on the nature of applied parameters, mobility of servers, type of variables, objective function(s), and applications.
Based on the uncertainty of the model parameters such as traveling times, demands, and costs, FLP could be classified into stochastic or deterministic problems [4,5].As instance, [6] studied a warehouse location problem subject to data contamination to address uncertain environment and human errors in disaster response.Li et al. [7] proposed a bi-level stochastic optimization model for selecting the shelter locations in the case of disaster evacuation based on dynamic traffic assignment.An et al. [8] studied a stochastic location problem integrated with en-route traffic congestion, in-facility queuing delay, and facility disruptions for disaster response planning.
According to servers, FLP could be subdivided into two categories: mobile and immobile servers.The mobile servers are predominantly used in emergency location problems.In this case, the servers travel to the location of customers [9][10][11].In contrast, customers are supposed to visit the servers in the case of immobile servers [12,13].
FLPs could be studied on discrete or continuous space problems.In discrete space problems, the facilities are to be located on a finite set of pre-defined candidate locations [14,15], but in the continuous space problems, any point in the plane could be opted for establishing facilities [16,17].
Based on the objective function, FLPs are categorized into center, covering, and median problems.Center problems minimize the maximum traveling distance or service times [18].Median problems minimize the total traveling or waiting times [19].Covering problems maximize the number of covered demands [20].These problems are classified into set covering and maximal covering location problems [21].Set covering aims to satisfy a specified level of coverage, while maximal covering tries to increase the possible coverage rate.
In contrast with the incapacitated problem with no bound for the services, capacitated location problems refer to those FLPs in which the required resources are limited and all the assigned demands may not be satisfied [22,23].The unsatisfied demands, also known as outliers, are defined as those portions of the demand that their satisfaction is not economically justified [24].The outliers will be referred to as "lost sales" in the rest of this paper.On the premise that the number of assigned customers cannot exceed the capacity of each facility, Torrent et al. [12] studied a capacitated location-allocation problem with lost sales to maximize the number of covered demands.The idea of lost sales has been rarely considered in the capacitated median problems.While respecting the capacity restrictions, the majority of related studies have assumed that there would be the possibility of covering all the demands using the available resources [7,10,25].Harris et al. [26] developed a bi-objective capacitated inventory-location model to minimize the environmental effects of established facilities and the financial costs simultaneously.The capacity restrictions reflect the maximum number of items that can be stored and the maximum number of customers that can be served.Referring to the extended enterprise as a set of entities with different capacity levels, such as internal units and external partners, Abbal et al. [27] studied a bi-level facility location problem considering customers-depots and depots-plants relationships at different levels.
The idea of congestion has been used rampantly in FLPs.In congested location problems, customers may wait in line to receive the service.As it is reviewed by Boffey et al. [28], congested location problems can be classified according to criteria such as the number of servers in each facility, rules for assigning customers to the facilities, and the probabilistic distribution of arrival and service times.Aboolian et al. [29] presented one of the first multiple server congested center location problems using the M/M/C queuing system to minimize the maximum traveling time plus the expected waiting time.Considering the same queuing framework, Sayarshad and Chow [30] proposed a non-linear programming model for the relocation problem of on-demand mobility services.Their model minimizes the travel time between customers and servers, the relocation costs to get the vehicle to the new locations, and the delay costs.
Incorporating capacity restrictions into the node-based congested location problems has attracted great interest in recent years.Applying the queue length and the number of servers restrictions, Tavakkoli et al. [31] proposed a multi-objective capacitated location model without lost sales.Shavarani et al. [32] studied a multi-level capacitated fuzzy facility location problem to determine the locations of launching and refueling stations in an aerial delivery system in which the Unmanned Aerial Vehicles (UAVs) are used as the delivery means.While considering a general distribution for the service time, they assumed that there are enough resources to cover all the demands.The idea of lost sales has also been considered in some of the studied capacitated congested location problems.Using an exponential service time in a single-server congested system with queue length restrictions (M/M/1/K), Berman et al. [33] studied a location problem to determine the minimum number and locations of required facilities such that the number of demands lost due to coverage radius and congestion restriction does not exceed a pre-defined threshold.Using the same queuing framework, Rahmati et al. [34] Table 1.A comparison of the present study with the most related works.

Ref.
Type Queuing Capacity Lost* Servers Solution Method Li et al. [7] node space-resource immobile Lagrangian relaxation Aboolian et al. [29] node M/M/C immobile GA, Descent Sinha [25] node service capacity immobile Add-Drop-Interchange An et al. [8] node M/M/1 service capacity immobile Lagrangian relaxation Pasandideh et al. [36] node M/M/1 service capacity immobile SA, GA, RS Abbal et al. [27] node service capacity immobile Branch and Cut Ilkhanizadeh et al. [10] node flying distance mobile Exact (CPLEX) Torrent et al. [12] node service capacity immobile Clustering-based SA Harris et al. [26] node service capacity immobile SEAMO2 Rahmati et al. [34] node M/M/1/K queue length immobile NSGAII, NRGA Etebari [35] node M/M/1/K service & queue immobile GA, PSO Berman et al. [33] node M/M/1/K queue length immobile Add-Drop, Tabu Search Tavakkoli et al. [31] node M/M/C/K queue length immobile MOVDO, NSGAII, NRGA Sayarshad and Chow [30] node M/M/C immobile Lagrangian decomposition Zamani et al. [13] node M/M/1 service capacity immobile Ant Lion, GA Shavarani et al. [32] node M/G/M endurance mobile MA Shavarani et al. [11] edge flying distance mobile GA, MA Jafari & Arkat [37] edge M/M/1 immobile SA, GA, MA Golabi et al. [38] edge flying distance combined GA, MA, SA Shavarani et al. [39] edge endurance mobile NSGAII, NSGAIII Golabi et al. [40] edge M/M/C total servers immobile GA, MA, SA This study edge M/M/C servers at each facility immobile MA, SMA, ADD, DROP considered a bi-objective FLP to minimize the establishing costs along with the summation of traveling and waiting times, simultaneously.The lost demand was only applied in their calculation of effective waiting time.The same queuing system was applied by Etebari [35] to maximize the total profit gained from covered demands.On the premise that customers leave the system in the case of facing a disrupted facility, Zamani et al. [13] applied an M/M/1 queuing framework with service capacity to maximize the profit gained from fully covered demands.Table 1 illustrates a summary of the aforementioned congested and/or capacitated location problems, as the main focus of the current study.Adverse to the conventional location problems in which a network of nodes is studied, flow capturing location problems (FCLP) and flow refueling location problems (FRLP) consider a network of edges each directed toward a specific flow [41].Inspired by FCLPs, the edge-based location problem has been recently introduced as an alternative to conventional node-based FLPs.It is believed that the idea of distributed demands along the network edges reflects a more realistic image of the problem.The edge-based location studies are considered a new branch of facility location problems.Although only a few studies focus on this variant, all the definitions, and categorizations used for node-based studies also work for edge-based facility location problems.One of the first works focusing on edge-based FLPs was studied by Jafari and Arkat [37].They proposed a novel model for a single-server congested location problem using an M/M/1 queuing system.Golabi et al. [38] developed a mixed-integer non-linear programming model for a stochastic edge-based median problem considering both mobile and immobile servers depending on the nature of the demand.Their objective function minimizes the traveling time of both customers and immobile servers over a set of feasible scenarios.The idea of distributed demands was also used in a study by Shavarani et al. [11] wherein regarding the flying range restrictions, a bi-level edge-based mobile location problem is modeled to determine the location of launching and refueling stations in a drone-based delivery system.Later, Shavarani et al. [39] considered a capacitated bi-objective mobile facility location problem with uniformly distributed demands in UAV-supported delivery operations to minimize the facility costs and lost sales simultaneously.
To the best of the authors' knowledge, the only study considering the congested capacitated edge-based location problem is done by Golabi et al. [40], where the authors introduced an expanded version of [37] by using a multiple-server environment and considering the maximum available servers restriction.Assuming a pre-known number of opened centers, the location of facilities and their corresponding number of servers were the decision variables.Despite using a restriction on the total number of assigned servers, the authors did not consider any limitation for the maximum number of servers at each location.They also assumed that the available servers have the potential of covering all the customers and thus no lost demand was considered.Using the M/M/C queuing framework, the objective function was to minimize the summation of customers' traveling and waiting times.
To reflect a more realistic image of the problem and account for the available space limitation, the current study extends [40] by applying a restriction on the number of servers assigned to each facility and incorporating the idea of lost sales.The classical constraints for the assignment of customers and the same formulation introduced in [40] for decomposing the network edges are applied in this study.Besides, new constraints are introduced to account for lost sales.It is also assumed that although all the customers travel to the location of facilities, due to capacity restriction, just a proportion of them would wait in the queue to be served and the rest are dismissed.The current study makes two contributions to the edge-based facility location literature.First, it proposes a new problem definition by incorporating the idea of lost sales into the congested location problems in a multi-objective context.Second, several modified solution approaches are developed to solve the problem.Besides, it is tried to present a comprehensive review of the related literature.For a fuller view, the main contributions of the this study could be listed as follows: -The number of opened facilities is not known and fixed in advance, and each facility has a threshold for the number of assigned servers to account for the available space restrictions.-Lost sales and waiting time for satisfied customers are optimized.
-A new mathematical formulation is developed to tackle the complications of the problem.
-A multi-objective version of the study with conflicting functions is introduced.
-New modified solution algorithms are introduced.
The objective functions that are considered include the aggregate customers' traveling distance, the aggregate covered customers waiting time, the total cost of the system, and the total number of lost sales.Here, a method based on the weighted distance metric is used for transforming the developed multi-objective problem into a single-objective one.
The complexity of traditional facility location problems has been thoroughly analysed in the literature.It has been shown that center problems, maximal covering, and median problem are NP-hard [42][43][44].Set covering location problems are also known to be NP-complete [42].Therefore, the application of metahuristics in solving location problems is prevalent [45].It is clear that the contributions mentioned in our work heightens the complexity of the studied problem and the concomitant solution methods.As the solution methods, this study proposes modified versions of some of the existing heuristics and metaheuristic algorithms widely used to solve NP-hard location problems.Memetic algorithm which is known as an effective metaheuristic to solve the edgebased location problems [11,37,38,40] is applied in this study.Add and Drop algorithms that are popular greedy heuristics able to find high-quality solutions efficiently [46][47][48] are also used as solution methods.Due to the complexity of the developed mathematical model and the need for optimizing the number of servers with the number and location of facilities, the performance of the original Add and Drop heuristics were not promising.Thus, this study develops modified versions of applied algorithms to improve the quality of obtained results.
The remainder of this study is organized as follows.Section 2 describes the problem and the mathematical model.Solution methods are described in Section 3. Numerical results are reported in Section 4 and the conclusion is drawn in Section 5.

Problem definition
The edge-based facility location problem considers the idea of distributed demands along the network edges.Due to the uncertain nature of demand generation, the location and the rate of demands on each edge follow a probabilistic distribution.Any point on network edges could be a potential point for demand occurrence.Thus, the uniform distribution could be used to model the position of demand generation points along the network edges.On the premise that the demands are uniformly distributed along the network edges, this paper studies a multiple-server multi-objective congested immobile facility location problem on the discrete plane.As a continuous distribution commonly used to measure the expected time for an event to occur [49], it is assumed that the time interval between the arrival of any two consecutive demands and the service time of each server follow exponential distributions.Therefore, the demand occurrence rate on each network edge follows a Poisson distribution with a rate dependent on the population density of the area and the length of the edge.These arrival and service time distributions are widely used in queuing systems 4 .Since the customers located on the same network edge (homological demands) have different distances to their closest open facility, the time interval between demand generations differs from the time interval between the arrivals of customers and as a result, the demands generated earlier may wait longer.The number of the opened facilities and their corresponding number of assigned servers are the decision variables limited by a restriction on the maximum number of servers at each facility.Since the customers travel to the location of the closest open facility, their movement along the network edges could be considered as M/G/∞ queuing systems with Poisson outputs [50].Therefore, as the demand enters the open facilities, the queuing system at each facility would be according to M/M/C.It is assumed that customers are served based on the First-in-First-out policy.The capacity of each facility is restricted by its total number of assigned servers and their associated service quality coefficient.The total number of demands that can be satisfied through each facility is limited and the excess demands that arrive at the facility but fail to receive services are considered as lost sales.It should be noted that the lost demand does not affect the calculations of the waiting time.Minimizing the aggregate customers' traveling distance, the aggregate customers' waiting time, the total cost of the system, and the total number of lost sales are the four objective functions addressed in this study.

Model formulation
In order to develop the mathematical model, the sets, parameters, and variables are defined as follows: Sets: -  : the minimum distance between node  and facility  -   : the cost of establishing a facility at candidate location  -  : the unit cost of servers - max  : the maximum number of servers that can be assigned to location  -: the common service rate of servers - max : the maximum time that each customer can wait in the system -: the coefficient of service quality level -: an infinitesimal positive number - : a large positive value

Decision Variables:
-  : 1 if node  is assigned to facility ; 0 otherwise -  : 1 if location  is selected for opening a facility; 0 otherwise -  : the number of servers assigned to facility  Output Variables: -  ′  ′ : the distance between node  and decomposing point of the edge (,  ′ ), where  and  ′ are assigned to opened facilities  and  ′ , respectively -  : the expected waiting time at opened facility  - 0  : the probability that opened facility  is idle -  : the demand entrance rate of opened facility  -  : the number of customers served at opened facility  According to Figure 1, if the closest opened facilities to the endpoints of a network edge are not identical, the edge would be decomposed into two different segments, each assigned to its closest opened facility.
As it is illustrated in Figure 1, on the premise that  and  ′ are respectively the closest opened facilities to the endpoints  and  ′ , the edge (,  ′ ) is decomposed into two segments.Customers located on the first segment with the length of   ′  ′ are assigned to facility , whereas the rest are assigned to the facility  ′ .Since the decomposing point has the same distance to the closest opened facilities of the endpoints (  +   ′  ′ =   ′  ′ +   ′  ′  ) and having the length of the edge (  ′  ′ +   ′  ′  =   ′ ), the length of each segment is calculated using equation 1: It is incontrovertible that   ′  ′ has a non-negative value less than or equal to   ′ , since on the contrary case, facilities  and  ′ would not be the closest open facilities to nodes  and  ′ respectively, and it is in contradiction with the considered assumptions.It is clear that if the closest facility to nodes  and  ′ is identical, the entire edge is assigned to the same facility.
Using the general formula for the probability density function of the uniform distribution, the expected traveling distance of customers located on the edge (,  ′ ) for arriving at the corresponding endpoint is calculated as: Based on the splitting property of Poisson distribution, in the case of decomposing edge (,  ′ ), the average traveling distance for the demands located on the segment   ′  ′ to arrive at the location of facility j is calculated using equation 3: As the result, the developed multi-objective mixed-integer non-linear programming (MINLP) model can be formulated as follows: Min Min Subject to: Objective functions (4), ( 5), (6), and (7) minimize the customers' aggregate expected traveling distance, the customers' aggregate expected waiting time, the total cost, and the total number of lost sales, respectively.Constraint (8) guarantees that each customer is assigned to only one opened facility.Assigning customers to closed facilities is refrained by constraint (9).Constraint (10) assures that each customer is assigned to its closest opened facility.Constraints (11) and ( 12) calculate the length of partitioned segments of each network edge and the demand entrance rate at each facility, respectively.Constraint (13) ensures that the servers are assigned just to the opened facilities.While assuring that the queuing system will reach a steady-state mode, constraint (14) reflects the capacity of each opened facility based on the number of assigned servers and their service quality level.Constraint (15) categorizes the customers to covered and lost sales.Constraints ( 16) and ( 17) measure the idle probability and expected waiting time at each opened facility, respectively.The maximum waiting time threshold is considered by constraint (18).Constraint (19) sets an upper bound for the number of servers assigned to each facility.Constraint (20) assures reaching the steady state mode.Constraints ( 21), (22), and (23) reflect the binary and non-negative decision variables.

Solution methods
Since the developed mathematical model is categorized as a constrained mixed-integer nonlinear programming model, it is considered as an NP-hard problem, so achieving the exact solution for such a problem is difficult, if not say impossible, specifically when the size of the problem is big [36].However, the first attempt to solve the problem was to code the model in General Algebraic Modeling System (GAMS) and solve it using exact MINLP solvers such as BONMIN, DICOPT, and KNITRO.In order to transform the four objective functions given in equation ( 4)-( 7) into a single one, the   -metric method is used as a weighting scheme in this study [51].Using equation ( 24), the method minimizes the differences between objective functions and their pertinent desired values ( *  ∈ R  ): Subject to: Where   , which is defined by the decision-maker, is the weight of the objective function   (),  is a positive integer, and  *  is the optimum value of objective function , if considered alone.If  = 1, which is the case of this study, the method is called the Manhattan metric.
It is noteworthy to mention that it is not possible to solve the model in its current form using GAMS.It is due to having a  function over the decision variables in constraint (14), a decision variable (  -1) as the upper-bound of the first summation operator in constraint (16), and the factorial of a decision variable (  !) in constraints ( 16) and (17).To get rid of these problems, a constraint manipulation was done by introducing a new decision variable called   where   : 1 if node there are  servers assigned to facility ; 0 otherwise.Now, constraints (13) is replaced with: Introducing a new set of binary decision variables called  1 and  2 , such that  1 +  2 = 1, constraint ( 14) could be replaced with the following set of constraints: In a similar way, constraints ( 16), ( 17), (20), and ( 21) are replaced with the following set of constraints, respectively.It is noteworthy to mention that constraint (19) is not required anymore, as the total number of servers is calculated using the summation of   over , where  runs from 1 to  max  .
The edge-based facility location literature suffers from the lack of benchmark instances required to check the validity of developed models, or to analyze the performance of the proposed solution algorithms.Except for the real-life case studies, all the previously mentioned works in this domain have generated random instances to discuss the performance of their solution methods.This study also uses some randomly generated instances to compare the results obtained from different methods.To generate semi-realistic instances, a desired number of network nodes are randomly positioned in a network.Using a random number of connecting links passing through each network node, the generated nodes are linked to form a connected graph.To simulate more realistic transportation networks, it is assumed that the number of connecting links varies from 1 to 8 while considering a higher probability for numbers between 2 and 4. A random number in range [1,3] is multiplied by the euclidean distance of all pairs of connected network nodes to account for the length of the corresponding network edge.The candidate locations are randomly selected among the existing network nodes.The cost of establishing facilities at each candidate location is a random integer in range [100K, 200K].The demand rate on each network edge is randomly chosen in range (0, 4].The cost of each server is fixed to 5000.The common service rate and maximum allowed waiting time are fixed at 25 and 0.5, respectively.Using this procedure, several small instances are randomly generated to check the performance of developed algorithms and compare the results with optimal solutions obtained from exact MINLP solvers of GAMS compiler. The results obtained from solving the model in the GAMS environment using different solvers and a time limit of 7200 seconds are shown in Table 2.The reported values are obtained from the Manhattan-metric method by considering equal weights for normalized objective functions.The normalization is done by dividing each objective function by its corresponding maximum value.Assigning each customer to its farthest facility would yield the maximum value of the first objective function.Considering the maximum allowed waiting time for each customer returns the maximum value of the second objective function.Opening all the candidate sites and assigning the maximum number of allowed servers yields the maximum value of the third objective function.The maximum value of the fourth objective is estimated by considering all demands as lost sales.According to Table 2, the obtained results revealed that the maximum solvable size of the problem is limited to 70 network nodes (endpoints of the edges).Comparing this size with the real city networks having thousands of network nodes is proof of the incapability of exact methods in solving such a complex problem for large instances.Since all the MINLP solvers of GAMS failed to solve the model for realistic instances of the problem, the genetic algorithm hybridized by a local search technique, called Memetic Algorithm (MA), along with its modified stochastic version called Stochastic Memetic Algorithm (SMA), as well as modified versions of two well-known heuristics in solving location problems, called modified Add and modified Drop algorithms are used to solve the proposed model.It is noteworthy to mention that the same results of solving small instances using GAMS solvers were also obtained by all of the developed solution algorithms (see Tab. 2).

Memetic Algorithm (MA)
Evolutionary algorithms combined with local search techniques are generally called memetic algorithms [52].This combination is used to merge the biological evolution with the individual learning procedure to mimic the cultural evolution [53].The individual learning procedures are exerted via local search techniques [54].In this paper, MA equips the Genetic Algorithm (GA) with a local search technique.The algorithm starts with a population of   number of randomly generated solutions called chromosomes.Each chromosome (S) is represented as a matrix of size 2 × , where  is the total number of the candidate locations.The first n genes of the encoding matrix form a binary vector representing open facilities by 1 and closed facilities by 0. The second part determines the number of assigned servers by a randomly generated integer number between 1 and the maximum number of servers for the corresponding facility.It is noteworthy to mention that the second part alleles corresponding to closed facilities is set to be equal to zero.Since the satisfaction of constraints ( 18), (19), and (20) could not be guaranteed using the mentioned structure, a repairing operation that may be followed by applying a penalty function is used to refrain from dealing with infeasible solutions in the initial population.The repairing operation tries to satisfy constraints (18) and ( 20) by iteratively adding one server to the facility for which the the current number of assigned servers cannot satisfy the steady state and waiting time threshold restrictions.If satisfaction of constraints (18) and (20) results in violating constraint (19), a penalty is added to the objective function.After selecting parents using the roulette wheel procedure, the two-point crossover operation depicted in Figure 2a is applied to merge the parents and generate two offspring.
In this study, one of the three different types of mutation operators is randomly selected for being applied to each generated offspring (Fig. 2b).The mutate operator selects a random number of genes in the first part of the chromosome and changes the status of the corresponding candidate location.The corresponding genes in the second part are filled in a manner described in solution representation.After selecting two random points along the chromosome, the reverse operator inverts the order of genes trapped between them.Finally, the swap operator considers two randomly selected candidate locations and exchanges their corresponding genes.The crossover and mutation operations are controlled by the probabilities called   and   , respectively.The repairing operator and the penalty function are considered for infeasible solutions generated.
Through successive iterations, the local search technique generates a predetermined number of neighborhood solutions () for each offspring generated via mutation operation.Neighborhood solutions are generated using the mutate technique described for the mutation operator.It means that not only the number of assigned servers of the selected genes may be changed, but also the status of the candidate location is reconsidered.It is noteworthy to mention that the repairing operator and penalty function are considered for infeasible solutions.In each iteration of developed MA, the generated solutions through the crossover, mutation, and local search operators are added to the population.After sorting the population based on the fitness value of the solutions, the best solutions are kept to resize the population to the initial  .The algorithm stops after reaching a predetermined number of iterations ( ).The pseudo-code of the applied MA is illustrated in Figure 3a.

Stochastic Memetic Algorithm (SMA)
The difference between the applied SMA and MA is in the structure of the local search technique.All the neighborhood solutions of MA are generated by applying the local search method to the solution produced by the mutation operator.On the contrary, each generated neighborhood solution in SMA has the chance of being the basis for applying the local search technique and generating other solutions.In other words, the local search in MA is applied on a fixed solution while inspired by the Simulated Annealing algorithm, the solution that is undergoing the local search technique is stochastically updated in SMA.To do so, first, the offspring generated via mutation operator is set to be the best solution of the local search technique (  ).
The generated neighborhood solution is compared to the best solution in terms of the objective function.If the generated neighborhood solution has a better objective value,   is updated and the next local search operator is applied on the updated   .Otherwise, based on a predetermined probability (  ) used for refraining from trapping in local optima, a decision is made to apply the next round of local search technique on the current solution or the   .The pseudo-code of the applied stochastic local search technique is illustrated in Figure 3b.Considering the applied mutation, crossover, and roulette wheel selection operations, the complexity of GA is ( (  × || +   × || +  )), which is equivalent to (  ×   × ||) [55].The overall time complexity of the applied memetic algorithm as a hybridization of GA with the local search operation is (  ×   × || × ) in the worst case.

Modified Add & Drop Heuristics
In the local search-based Add [56] and Drop [57] heuristics, facilities are successively placed in or deleted from locations, dependent on some measures of suitability.In this way, different combinations of selected locations are checked in the search procedure [25].
Drop heuristic initiates by opening all the candidate locations with the maximum number of servers.The objective function is calculated and then the algorithm determines if closing any of the facilities leads to cost savings or improvement of the objective function.If so, the selected facility will be closed and the algorithm continues to search for the next facility to close.In case there are several facilities that their drop from the solution leads to the improvement of the objective function, the one with the highest impact is selected.If no cost saving is possible, the algorithm terminates.
Add heuristic starts with all the facilities closed.Then the algorithm tries to open facilities one at a time.The facility that is selected is the one with the highest gains in the objective function.The candidate location which yields the best objective function is opened and the algorithm proceeds to the next step where it would look for the second candidate location to open beside the first one.The algorithm continues to increase the size of the selected candidate locations while the objective function is improved compared to the previous step.Thus the algorithm terminates in a finite number of iterations.
In this study, the maximum number of servers is assigned to the facilities.In the next step with the set of opened facilities at hand, a greedy algorithm is utilized to justify the total number of assigned servers for each facility.In this respect, facilities are selected one by one and the total number of servers is decreased continuously while both constraints (18) and (20) are still satisfied and the quality of objective function improves, otherwise, the greedy algorithm moves to the next facility.The complete solution comprised of opened facilities and their respective number of servers is returned by the algorithm.The pseudo-codes of the modified Add and Drop heuristics are illustrated in Figure 4.
The add heuristic starts with the determination of the cheapest facility.First, each facility is considered independently and as the only open facility.Having || potential facilities, the cost evaluations are repeated || times in the first step to select the cheapest facility (lets say  1 ).In the next step, we will consider the combination of  1 with the remaining || − 1 facilities.The process goes on as described in pseudocode 4. Thus, the total number of cost evaluations for selecting open facilities is ||(|| − 1)/2.For each facility we consider different number of servers ∈ {1, . . .,  max  },  max  being the maximum number of servers for  th facility.Considering  max = max ∈  max  , the maximum number of evaluations for the Add heuristic will be  max ||(||− 1)/2 .Accordingly, the complexity of this algorithm is equivalent to ( max ×|| 2 ).The drop heuristic complexity is equivalent to the ADD heuristic.

Parameter Tuning
The applied algorithm parameters affect the performance of developed metaheuristic solution algorithms.Therefore, tuning the parameters before the implementation will improve the quality of obtained solutions.This study applies the Taguchi method to tune the algorithm parameters.Taguchi uses orthogonal arrays to investigate a large number of controllable factors through a small number of experiments [32].Taguchi determines  the best level of parameters via minimizing the noise effect.The variation of response is evaluated by a signal to noise ratio / that is calculated by using equation (35), in which Y denotes the response value and n represents the number of orthogonal arrays.
The applied MA and SMA parameters with their corresponding ranges are shown in Table 3.Using each parameter combination shown in Table 4, five randomly generated instances with different specifications are solved five times and the average of obtained objective values are considered as the response value of the applied algorithms.The algorithms are coded in Python 3.8 and implemented on an Intel Core i7-6500U CPU @ 2.50 GHz laptop with 8 GB RAM, 6 MB L3 Cache, and 1 MB L2 Cache.The Minitab 19 is used to perform the Taguchi method.The obtained / rations for applied MA and SMA are shown in Figure 5.As it is illustrated in Figure 5a, the SMA parameters should be set at their third, third, second, third, second, and third levels respectively.

Experimental results
All of the existing FLP benchmark datasets are designed for node-based location problems.To translate such benchmarks to edge-based instances, for each node one should consider a random number of edges to select directly connected nodes, random edge lengths, and random demands on each generated edge.There is no doubt that this would completely change the nature of the considered benchmark dataset.Thus, generating appropriate instances to evaluate the performance of applied algorithms would be inevitable in edge-based location problems.Using the random instance generation procedure and the Manhattan-metric method for normalized objective functions of equal weights explained in Section 3, the results of implementing the proposed algorithms on 16 randomly generated test cases of different sizes and specifications are shown in Table 5.The average objective value and required CPU times of 10 runs of each proposed algorithm on each test case are shown in this table.In order to discuss the potentials of implementing the solution algorithms on real-life large-scale networks, all the test cases are generated in large size networks with a plethora of nodes and edges.To simulate a more realistic image, the number of edges passing through each network node varies between two and five.Generally, the downtowns are the main crowded areas of real-life case studies.To simulate this behavior, it has been tried to consider higher demand rates for the central edges, in comparison with those that are closer to the network barriers.
As the GAMS optimization compiler failed to solve the developed MINLP model even on small instances, the performances of the proposed algorithms are compared together.Using the same number of evaluations, Figures 6 and 7 provide the graphical comparison of applied algorithms in terms of obtained mean objective value and required CPU times, respectively.Additionally, the one-way analysis of variance (ANOVA) for the As it is shown in Figure 8, the worst results are obtained by the Add heuristic.Although there is no significant difference between the results obtained by MA and Drop algorithms, MA could be considered as a more efficient algorithm due to requiring fewer CPU times.The developed SMA slightly outperforms MA and Drop algorithms.It could be said that for the case of emergency applications such as determining shelter or relief distribution centers, MA could be considered as the most efficient algorithm to determine reasonable solutions.However, in strategic decision-making procedures where the required CPU time is not among the main criteria to evaluate the performance of applied algorithms, SMA could be considered as the best solution approach to find high-quality solutions.

Conclusion
Although the literature is replete with studies focusing on different types of node-based facility location problems, there are just a few studies concerning the idea of distributed demands along the network edges.Since the edge-based location problems could present a better image of the real-life problems, the necessity of enriching the literature by scrutinizing this specific type of location problem is more highlighted.Thus, this study presents a novel multi-objective mixed-integer non-linear programming model for a congested capacitated facility location problem with uniformly distributed demands along the network edges.The demand occurrence rate on each network edge follows a Poisson distribution, while the service time of each server follows an exponential distribution.It is assumed that the number of servers at each location is limited by a predetermined threshold.Therefore, just a fraction of demand could be satisfied and the remaining is considered as the lost demand.This paper considers an M/M/C queuing system with a maximum waiting time threshold for the satisfied demands.This study aims to determine the number of facilities, their location, and the number of servers assigned to each facility such that the traveling distance, the waiting time, the total cost, and the number of lost sales are minimized simultaneously.
Due to the NP-hardness of the proposed mathematical model, a stochastic memetic algorithm, a memetic algorithm, and modified add and drop algorithms are developed as the solution methods.After tuning the parameters, algorithms are implemented to solve four randomly generated large-scale test cases and the obtained results are compared using the ANOVA test.According to the obtained results, SMA outperforms the developed MA, Add, and Drop algorithms.
Future studies may check the possibility of obtaining LP-relaxed optimal solutions as a bound to compare the solution algorithms.They also may use different queuing systems by considering different distributions for the demand generation and service time, or by restricting the length of the queue to account for the capacity of servers based on the number of customers arrive at each center at the same time.Furthermore, using fuzzy inputs such as the cost of establishing the facilities or implementing the servers and the rate of demand or service time could be of great interest.Also, one may use other metaheuristic or heuristic solution algorithms to obtain better results.

Figure 2 .
Figure 2. a) The crossover operation; b) The mutation operation.

Figure 3 .
Figure 3.The pseudo-code of a) MA; b) stochastic local search technique.

Figure 4 .
Figure 4.The pseudo-code of a) Drop Heuristic; b) Add Heuristic.

Table 2 .
Computational results for small instances.

Table 4 .
Computational results for tuning MA and SMA.