1. INTRODUCTION
In many businesses and organizations, facilities are categorized into multiple levels, each offering a different service level. These multilevel hierarchical structures are common in both public and private sectors, for example in systems that deliver health and education services, bank branches, ATMs, etc. (Marianov and Serra, 2001; Amouzad Mahdiraji and Yousefi Talouki, 2021; Nasiri and Ahmadi Daryakenari, 2021; Sharma and Sharma, 2021). Facility location is a very important issue for these systems. However, in many realworld cases of hierarchical location, the service time of facilities does not match the volume of demand, resulting in considerable queuing or congestion. To address this prevalent problem, the generated queue systems need to be analyzed in terms of criteria such as queue length or average waiting time of customers. Because of the uncertainty involved in the time that it might take to serve arriving customers, the problem of location of congested facilities could be difficult to model and solve. The present study aims to bridge the gap in the literature on the multilevel congested facility location and the impact of prioritizing services compared to using nonprioritized (turnbased) queues. The rest of this paper is structured as follows. In the next part of this section, we review the literature on the subject of hierarchical facility locationallocation and congested facility locationallocation. In section 2, we present the problem statement and the proposed formulation and describe the model indices, parameters, and assumptions. In section 3, we provide a detailed description of how the model is solved and report the results for a numerical example. Finally, in Section (4), we present the conclusions and offer some suggestions for future research.
In a study by Mitropoulos et al. (2006) they formulated a hierarchical facility location model with the objective of maximizing coverage and equity in a healthcare system with a twolevel structure comprised of local health centers and hospitals. They assumed that efficiency and equity of distribution of facilities among citizens are the two most important factors for the optimal design of healthcare systems. The efficiency was optimized in the form of an objective function seeking to minimize the total distance between patients and facilities. The equity of service distribution was optimized by an objective function that minimizes the maximum distance between facilities and patients. Hodgson (1986) modeled the location of primary healthcare facilities in a threelevel structure with the objective function of minimizing the total distance between demand points and the closest facility with the appropriate service level. The model of this study is a logical yet practically unrealistic hierarchical model in which demand at every level is met at the nearest appropriate facility. This study also presented an interactive model in which the objective is to minimize the distance of individuals from the higherlevel facility they anticipate to be referred to from the first place they visit, and help choose the first place based on the sum of distances in question over the entire hierarchical system.
Since the early 1940s, many researchers have proposed mathematical formulations for more realistic modeling of customer travel behaviors. These include researchers such as Mitropoulos, Hodgson, and Jacobsen, who have point out the unrealistic nature of the assumption that customers will visit the nearest facility. According to Mitropoulos et al. (2006) the attractiveness of facilities depends not only on proximity but also on the quality of services provided. Hodgson and Jacobsen proposed an interactive model called the expected referral distance, which models several types of customer travel behaviors based on their expected referrals. Teixeira and Antunes (2008) also stated that in the location of public facilities, it is very important to pay attention to the limitations of customerfacility allocation and analyze the spatial patterns of this allocation.
Fard and HajaghaeiKeshteli (2018) modeled the problem of multilevel hierarchical locationallocation with reverse and forward flows between different levels. In this study, researchers considered a series of customer zones, distribution centers, and recovery centers in the model structure and formulated it with the help of Stackelberg game theory. The model was solved with three wellknown algorithms called Tabu search, particle swarm optimization, and variable neighborhood search, and also with two new algorithms called the Keshtel algorithm and water wave optimization. Since the performance of metaheuristic algorithms highly depends on their parameters, the optimal parameter settings for these algorithms were determined with the Taguchi method.
In a study by Mousazadeh et al. (2018) they proposed a model for the problem of designing a health service network with a threelevel hierarchical structure. The proposed objective functions were focused on minimizing the cost of building the facilities and minimizing the distance between the demand points and the facility at each level. The proposed model was formulated in the form of a nonlinear programming problem and solved with four robust optimization approaches. Finally, the performance of these four approaches was analyzed in a case study.
Zarrinpoor et al. (2018) developed a hierarchical location allocation model for designing a health service network. This model was constructed with a twolevel multiflow hierarchical structure that takes into account the uncertainty of demand and services and geographical access. In this model, patients are prioritized based on urgency in receiving services and service quality is measured in terms of the expected waiting times considering the priority queuing system and risk of unexpected disruptive events.
One of the earliest studies explicitly focused on congested facility location is the research carried out by Wang et al. (2002). In this study, they proposed a facility location model with immobile servers where a maximum of one server can be deployed in each facility. The queue model developed in this study is M/M/1.
Berman and Drezner (2006) expanded the model of Wang et al. from singleserver to multiserver mode and investigated the performance of M/M/c systems in congested facility location with immobile servers.
The model presented by Pasandideh and Niaki (2012) is one of the few models where both customer and service aspects are simultaneously considered in the objective function. In this study, a biobjective model was developed to simultaneously minimize the total time spent by customers (travel times plus waiting times) and the idle time of service providers. In this model, it is assumed that only one server can be deployed in each facility, which makes the proposed queue model of the M/M/1 type.
Chambari et al. (2011) also presented a biobjective model in which customers and service providers are both simultaneously considered. The objective functions of this model simultaneously minimize customer travel and waiting times and the idle rate of service providers. The main assumption that distinguishes this model from the previous model is the limited capacity of waiting space in each facility, which makes the queuing system of the M/M/1/k type.
In a study by Hamaguchi and Nakade (2010), they modeled the problem of location of fixed facilities with the M/G/1 queuing system, in which times between consecutive customer arrivals follow an exponential distribution, but service times have a general distribution.
Harewood (2002) modeled the problem of location of multiple emergency service/ambulance facilities with multiple objective functions including the maximization of the covered population and minimization of the demand coverage costs. Also, the Queuing Probabilistic Location Set Covering Problem approach was used to estimate the probability that all service providers in the area are busy. For this purpose, the facilities were assumed to behave according to the M/M/m/K model, and having the minimum required confidence level and the probability that all employees of each facility are busy, the minimum number of employees required to cover the demand points with a specific probability was calculated.
Aboolian et al. (2009) formulated a new model for server facility location with the objective function of minimizing the maximum time spent by customers, including travel time and waiting time. In this study, the facilities were considered as M/M/m queuing systems and service allocation was also incorporated into the model formulation.
In another study, Aboolian et al. (2008) modeled the locationallocation problem in a congested network with the objective of minimizing systemwide costs (including fixed facility construction costs and variable costs associated with travel and waiting). In this study, The queue model used for the facilities was of the M/M/c type.
Rahmati et al. (2013) developed a location model with an M/M/c queuing system and server allocation, in which the maximum number of servers per facility was limited. They also incorporated the cost of allocating staff to the facility and the service level into the model formulation.
In a study by Hajipour et al. (2016), they presented a multiobjective model for multilevel congested facility location within the M/M/1 queue model. This model has three objective functions: minimizing the total travel time and waiting time of customers, minimizing the maximum probability of facility idleness, and minimizing the number of facilities (or in other words, minimizing the cost of constructing facilities). In this model, the system is structured in such a way that the facilities of different levels provide different services and each customer must receive service from facilities of all service levels. They solved the proposed model with two metaheuristic methods, namely the multiobjective vibration damping optimization and the multiobjective harmony search algorithm, as well as the nondominated sorting genetic algorithm and multiobjective simulated annealing and compared the obtained results.
TavakkoliMoghaddam et al. (2017) presented the new facility location model with the M/M/c/K queue model combined with pricing and facility capacity determination. In this work, in addition to location and allocation variables, the number of employees and queue capacity, and the price of service in each facility were included in the model and the number of employees and capacity in each facility was assumed to be limited. This model defines the satisfaction of customers as a function of the price of service/product and their distance from the facility, with different customers having different levels of sensitivity to this price and distance. Therefore, the queue length in this model is controlled by the pricing policy and the number of service providers in each facility.
Pasandideh et al. (2013) formulated a multiobjective facility location model with batch arrivals and an Mx/M/1 queuing system. In this model, each batch of customers is assigned to a facility and the size of these batches is a random variable. To balance the goals of customers and service providers, they considered three objective functions that minimize the total travel time and waiting time of customers, minimize the maximum idleness of each facility, and minimize the total fixed cost of constructing facilities. They considered different weights for the terms in the first objective function. These researchers stated that minimizing the average probability of idleness does not necessarily lead to minimizing the probability of idleness of all facilities, but minimizing the maximum idleness reduces the likelihood of idleness of all facilities. To improve the quality of service, they placed a coefficient alongside the service rate of servers in the capacity constraint.
In another study by Hajipour et al. (2014) they modeled a facility location problem with the allocation of service providers to facilities. In this study, they introduced the notion of reliability to the M/M/c queuing system and developed a model with the objectives of maximizing system reliability, minimizing system costs, and minimizing waiting time. The model was solved with multiobjective vibration damping optimization algorithms, multiobjective genetic algorithm, and simulated annealing algorithm. The results showed that the vibration damping optimization algorithm is more efficient than the other two metaheuristic algorithms.
Araz et al. (2014) proposed a pmedian location model for a certain number of facilities to provide urban public health emergency services in Arizona. The model was developed for the location of dispenser centers deployed in response to an anthrax attack, the allocation of demand points and human resources to these centers. Each dispenser center in this network was considered as an M/G/c queuing system with the objective of minimizing the total average waiting time and the total travel time of customers, which were assumed the main criteria for measuring system performance.
Khezerlou et al. (2021) developed a model for optimizing a biomasstobiofuel supply chain for reliability under different risk items.
2. METHODOLOGY
The considered locationallocation problem is a hierarchical problem where each level provides its own services and the facilities at the same level provide the same type of services. In other words, the facilities at each level are similar to each other. To be served completely, customers must visit one facility from each level. Customers are assumed to be positioned in specific predetermined locations, which are divided into several demand points. Customer demand is probabilistic and follows a Poisson distribution with a known mean. There are a certain number of potential locations for establishing facilities of each level. In this problem, the goal is to select suitable places for establishing a suitable number of facilities from among potential locations and then allocating customers to the facilities of each level using two queue modes: prioritybased queue and turnbased queue. In this model, the possibility of congestion in facilities is also taken into account. Since customer demand follows the Poisson distribution, the demand allocated to each facility, which is the sum of multiple Poisson distributions, also has a Poisson distribution. In addition, each facility is assumed to have a server with an exponential service time and therefore an M/M/1 queuing system. The facilities do not have capacity constraints and the longterm efficiency coefficient of the system is assumed less than one. The conceptual model of the considered problem is illustrated in Figure 1.
In the following, we first present a multiobjective mixed integer nonlinear programming model, where the first objective function minimizes the total waiting time of customers in the facility queues and the second objective function minimizes the maximum facility idleness. Obviously, these objectives are conflicting, because as the idleness probability decreases, the waiting time increases and vice versa. These objective functions create a balance between the goals of customers and service providers because the first function works in the customers’ favor while the second function works in the favor of service providers.
2.1 Problem Assumptions

 Each customer receives each level of service only once and cannot revisit a facility or opt out of receiving a service.

 Each facility operates as an M/M/1 queuing system.

 The time between customer arrivals and service times are independent of each other.

 The efficiency coefficient of each facility is less than one.

 Service providers (facilities) are immobile and customers must visit them to receive the offered services.

 All demands must be fully satisfied (demand loss is not allowed).

 In the prioritybased queue, customers with higher priority are placed in front of the queue, but in the turnbased queue, services are provided based on the firstinfirstout (FIFO) policy.
2.2 Model Sets and Indices

i : Set of customer points (demand); i=1,…,M

j,j’ : Set of potential locations for the facilities of each level; j,j’=1,…, N

l : Index of facility level; I=1, …, K
2.3 Model Parameters
2.4 Model Decision Variables

y_{jl}=1 if facility j is constructed at level l, and y_{jl}=0 otherwise.

x_{ijl}=1 if customer i is allocated to facility j of level l, and x_{ijl}=0 otherwise

z_{1} : Minimizes the total waiting time of customers in facility queues

z_{2} : Minimizes the maximum facility idleness

λ_{jl} : Average arrival rate of customers at facility j of level l

wq_{jl} : Waiting time of customers in facility j of level l

Wqjl^{(i)} : Waiting time of customer with priority i in facility j of level l
2.5 Mathematical Formulation of the Problem
Equation (1) minimizes the total waiting time of customers in the facility queue. Equation (2) minimizes the maximum facility idleness. Obviously, there is a conflict between these two objective functions, because as the service rate of a facility increases, although the total waiting time in the queue decreases, the probability of facility idleness also increases (and vice versa). Therefore, facilities need to be established in such a way as to reach a balance between these objectives. Equation (3) calculates the average customer arrival rate at each facility. Equation (4) guarantees the longterm stability of the system. Equation (5) computes the probability of a facility remaining idle over a long period. Equation (6) guarantees that at least one facility should be established at each level. Equation (7) ensures that customers should only be allocated to the facilities that are established, and that each established facility must have at least one customer and at most, the total number of customers allocated to it. Equation (8) ensures that each customer visits only one facility from each service level. Equation (9) ensures that customers receive all levels of service. Equation (10) defines the range and type of decision variables.
2.6 Analysis of Waiting Times in Different Queuing Systems
Given the mathematical complexity of the prioritybased queuing system, there are limited results for models of this system. In the model considered in this study, customers arrive at facilities of different levels based on the Poisson process with parameter (λ). Regardless of their priority, all customers receive the service at the same rate, which is equal to (μ). In this model, customers are prioritized in levels 1, 2, 3,..., M. Customers with priority 1 receive service ahead of customers with priority 2 and by the same token, customers with priority M1 receive service ahead of customers with priority M. The number of servers in each facility is 1. Therefore, based on the above assumptions, the waiting time of a customer with priority k is given by the following equation (Tirkolaee et al. 2020).
In the turnbased system, customers arrive at facilities of different levels based on the Poisson process with parameter (λ). The service rate is the same for all customers and is equal to (μ). In this model, customers do not have any priority over each other and receive the service based on the FIFO policy. The number of employees in each facility is 1. Based on these assumptions, the waiting time of customers can be obtained from the following equation (Mohtashami et al., 2020).
3. RESULTS AND DISCUSSION
To evaluate the model, in this section, several numerical examples of different dimensions are solved and the results are analyzed. Since Sherali and Nordai have proven that the multifacility locationallocation problem (where facilities can be located anywhere in Euclidean space) with definite parameters is an NPhard problem, the best approach to solving these problems is to use heuristic and metaheuristic methods. Therefore, a multiobjective nondominated sorting genetic algorithm is used to solve the proposed model.
Since the authors found no similar model in the literature, the data for testing the model were generated at random. In addition, since in most articles on locationallocation, the required data have been produced with a uniform distribution, this article also used a uniform distribution for this purpose. The stochastic parameters of the problem are given in Table 1.
The nondominated sorting genetic algorithm (NSGA) has been developed to solve multiobjective problems. The second version of this algorithm, which was developed by Deb et al. (2000) uses the concept of crowing distance instead of sharing function and employs an elitist approach to ensure diversity in Pareto solutions. To solve the proposed model in its different dimensions, this study used the second version of this algorithm, which is called NSGAII.
In this paper, a special method is used to increase the feasibility of solutions and satisfying more constraints. The chromosomes representing the solutions are structured in the form of a matrix where:

1. There are as many rows as there are customers.

2. There are as many columns as there are service levels.

3. Each element is a random number between 1 and the total number of facilities in the corresponding level.
The structure of this chromosome is shown in Equation (16).
In this chromosome, the number of rows is equal to the number of customers (M) and the number of columns is equal to the number of service levels (k). Each element of this matrix is a random number between one and the number of facilities in that level (J). To decode the chromosome, the variables related to the customers allocated to the facility are set to one. With this allocation, the facility will be constructed and its corresponding variable will be set to one.
Since this way of defining chromosomes satisfies most of the model constraints, the solutions generated in each iteration of the algorithm may not necessarily be feasible. This way of defining chromosomes ensures that constraints (6) to (9) are always satisfied. One of the common methods for dealing with such issues in constrained optimization problems is to use a penalty function to turns these problems into unconstrained problems. Since the model presented in this study is multiobjective, the penalty function should be added to both objective functions. The penalty functions used in this study are:
Where U is a large positive value, g(x) is the target constraint, and p(x) is the penalty considered for the infeasible chromosome. It should be noted that the above equation has been designed for constraints where (x) ≤ b.
Using a sorting algorithm, the solutions obtained in each population must be converted to a local Pareto front. In this research, the nondominated sorting algorithm was used for this purpose (Coello et al., 2007).
In the genetic algorithm, the members of the parent population were produced at random and the population of children was generated using the singlepoint crossover and mutation operators.
The solutions obtained from metaheuristic methods are sensitive to the values selected for the algorithm parameters (Rahmati et al., 2013). In this study, the optimal values of the parameters of the multiobjective NSGAII for problem instances were obtained using the Taguchi method at three levels. After reviewing the previous studies and some trial and error, the parameter levels were defined as shown in Table 2.
According to the standard Taguchi table, for four factors with three levels, one can use either L9 or L27 arrays (Montgomery, 2017). In this study, the L9 array was used as it involves fewer calculations. Finally, using the Minitab software, the optimal values of algorithm parameters were determined as shown in Table 3. The objective function values in the solutions of multiobjective NSGAII with these parameter settings for problems of different sizes with prioritybased and turnbased queue systems are reported in Tables 4 and 5, respectively.
As shown in Figure 2, using the prioritybased queues reduces not only the total waiting time of each customer in the queue but also the total waiting time of all customers.
4. CONCLUSIONS
The hierarchical congested facility locationallocation problem is a widely used type of location problem that has attracted the attention of many researchers over the past two decades. This study investigated the effect of considering prioritybased queues and turnbased queues in the serviceprovider facilities on this type of locationallocation problem. Since this problem is NPhard, the multiobjective NSGAII was used to solve the model. The optimal parameter settings for this algorithm were determined using the Taguchi method. Given the obtained results, it can be stated that if the goal of hierarchical facility location is to reduce the waiting time of a particular group of customers, it is best to establish a prioritybased queuing system or in other words, give that group of customers higher priority, as it will reduce the average waiting time of the entire system. Considering prioritybased queues in the service delivery system will mean that once arrived at a facility, customers that have a higher priority will be automatically placed in front of the queue, replacing the customer with lower priority and leaving his service halffinished. Therefore, this could be an attractive subject for future studies.