1. INTRODUCTION
In most cases, there are several levels of facilities, each providing a different level of service, compared to the other levels. Such systems with a multilevel or hierarchical structure are common in the public and private sectors, including school systems, bank branches, ATMs and public health services (Marianov and Serra, 2001; Roozitalab, 2020). Therefore, it is crucial to evaluate this type of location problem. However, facility service time compared to consecutive customer visits is considerable and the formation of queues or crowds is inevitable in most realworld applications of hierarchical location problems. One of the important aspects of considering this significant factor in most realworld applications is analyzing the queue systems generated and criteria such as mean waiting time in the queue and system (Goli and Malmir, 2020). Given the uncertainty in times between customer arrivals and service times, the location problem of crowded facilities is recognized as a complicated issue for modeling and solution. With regard to studies performed in this area, we consider the research gap in the area of multilevel facilities in the location of crowded facilities along with the comparison of the effect of priority in customer service provision with nonprioritized service provision (shift system) in the present study. The remainder of the research is structured as follows: the next section evaluates the relevant literature and reviews multiple studies in the field of locationhierarchical allocation problems and crowded facility locationallocation problems. Section 2 expresses the problem and the problem model and its indices, parameters, and premises. Section 3 provides details of the model solution and reports results, and Section 4 concludes and makes suggestions for future studies.
2. LITERATURE REVIEW
2.1 LocationHierarchical Allocation Problems
Mitropoulos et al. (2006) modeled (nested hierarchical location by considering maximum covering the distance in the health field with a twolevel structure, including health center levels and local hospitals. They considered two factors of efficiency and equal distribution of facilities among citizens to be important for an optimal design of a healthcare system and expressed the efficiency factor in the form of the objective function of minimizing the total distance between patients and facilities and the equal distribution factor in the form of the objective function of minimizing the maximum facility distance between patients and facilities. Hodgson (1986) focused on the modeling of primary healthcare facilities with a threelevel nested structure and considered the objective function of minimizing the total distance between the demand points and the nearest appropriate facility for all demand levels. The result was a logical hierarchical system, in which demand at every level is met at the nearest appropriate facility, but is unrealistic in practice. In addition, they proposed an interactive model, the objective of which was to minimize the distance predicted by individuals while considering the possibility of being referred from the first place to another higherlevel location. Moreover, decisionmaking was carried out based on the sum of the expected reference distances in the whole hierarchical system.
Since the early 1940s, several models have attempted to rationally justify the travel behavior of customers. For instance, researchers such as Mitropoulos et al. (2006), Hodgson (1986) have referred to the unrealistic assumption of the referral to the nearest facility according to the theory of spatial interaction. In fact, Mitropoulos et al. (2006) believe that facility attractiveness to people depends on the quality of services provided by various facilities in addition to the proximity criteria. On the other hand, Hodgson proposed an interactive model entitled expected referral distance, which exhibited several types of customer travel behaviors and was designed based on the expected referral distances. In a research, Teixeira and Antunes (2008) emphasized the importance of attention to constraints of customer allocation to facilities in general facility location problems, and analyzed a spatial model for user allocation to facilities. Fard and HajaghaeiKeshteli (2018) addressed a hierarchal nested allocation problem while simultaneously considering reverse flow between different levels and forward flow. They considered customer zones, distribution centers, and recovery centers in their model’s structure and carried out the modeling process using the Stackelberg game. To solve the model, they introduced two new models of the Keshtel algorithm and water wave optimization algorithm in addition to three algorithms of tabu search, particle swarm algorithm and neighbor search and solved the model using the two new algorithms. Since metaheuristic algorithms are highly dependent on their parameters, the optimal values of the algorithms were determined using the Taguchi design. In another study, Mousazadeh et al. (2018) focused on a hierarchical threelevel health service network design problem with the objective of minimizing the total establishment cost and total weighted distance between patient zones and health facilities at each level. The model was presented in a mixedinteger nonlinear programming framework, and four robust optimization approaches were applied to solve the model. In the end, the performance of the approaches was evaluated and analyzed. Zarrinpoor et al. (2018) focused on a hierarchical allocation problem in health service network design. Their model had a twolevel and multiflow structure and was solved while considering uncertainty in demand and services and geographical access. In this model, patient prioritization was considered, meaning that the priority was to serve the emergency patients. In addition, service quality was considered in the average expected time of priority queues and the risk of disruptive events.
2.2 Crowded Facility LocationAllocation Problems
In research, Wang et al. (2002) studied a crowded facility location problem with immobile servers for modes in which a maximum of one server is based in each facility. Notably, the M/M/1 queueing system was used in the foregoing study. By developing Wang’s model into a multiserver facility, Berman and Drezner (2006) evaluated M/M/c systems in congested facility locations with immobile servers. In research, Pasandideh and Niaki (2012) presented a rare model that simultaneously considered both aspects of customer and server in the objective function. They presented a biobjective model for simultaneous minimization of total time spent by customers (travel times and waiting times in the system) and idletime percentage. The model’s premise was the establishment of only one server in each facility, which meant that the queued model was of M/M/1 type. Chambari et al. (2011) proposed a biobjective model with simultaneous considering of two customer and server aspects. The objective functions of the model included the simultaneous minimization of travel times and waiting times of customers and idle percentage of servers. The main premise of the model, which distinguished it from previous models, was considering constrained waiting space capacity in each facility, which led to the use of M/M/1/k queue systems.
Hamaguchi and Nakade (2010) considered a facility location problem with fixed facilities based on an M/G/1 queuing system, and intervals between consecutive inflows of customers followed an exponential distribution but the serving times were considered as a whole distribution. Harewood (2002) conducted a research to model ambulance location problems by modeling the partial coverage problem with a specified number of facilitation and multiple objective functions. The objective functions included the maximization of the covered population and minimization of demand covering costs. In addition, a possible queueing approach was used for calculating the probability that all servers in a given region are busy. To this end, each facility’s behavior was considered to be in the M/M/m/K format and the minimum number of servers required to cover the demand points with a certain possibility was calculated following determining the minimum confidence level and calculating the probability of business of all servers in each facility.
Aboolian et al. (2009) introduced a new model for the network location of the center with the objective function of minimizing the maximum time spent by customers, including the travel time and waiting time. Each facility was considered as an M/M/m queue system, and server allocation was taken into account in the problem. Moreover, evaluated the locationallocation problem in a congestion network based on the objectives of minimizing the total costs of the system (including the fixed facility development costs and variable costs related to customer travel and waiting expenses). The M/M/c model was considered for each facility in the mentioned research. In another study, Rahmati et al. (2013) proposed a location model within a multiserver queuing framework, in which facilities behaved as M/M/m queues. They also considered server staffing costs and service level costs in the model. Hajipour et al. (2016) suggested a model with multiple objective functions for a multilevel congestion facility problem within the framework of the M/M/1 model. Their proposed model included three objective functions of minimizing the total customer travel and waiting time in the system, minimizing the maximum possibility of facility idle time, and minimizing the number of facilities the cost of facility development. The system’s structure was such that the facilities at different levels provided various services and each customer should receive service from facilities of all levels. They solved the model with a metaheuristic multiobjective vibrationdamping optimization (MOVDO) and the multiobjective harmony search algorithm (MOHSA), nondominated sorting genetic algorithm (NSGA) and multiobjective refrigeration simulation and compared the obtained results.
Using the M/M/c/K queue model along with adding pricing and determining the capacity of facilities, TavakkoliMoghaddam et al. (2017) presented a new model, in which the number of servers, queue capacity and service delivery price of each facility were considered in addition to the location and allocation variables. In addition, the number of servers and capacity of each established capacity was constrained. Their proposed model was such that it considered customer utility as a function of service or product cost and their distance to the provided facility. In addition, they considered the sensitivity of different customers to the distance and price to be diverse. Therefore, the pricing policy and the appropriate number of servers in each facility controlled the queue length in the problem. In another research, Pasandideh et al. (2013) considered the batch arrival queuing framework and proposed a model within the framework of M^{x}/M/1 queue systems. In this system, each group of customers is allocated to one facility, and the size of the groups is a random variable. Three objective functions considered in the foregoing study included minimizing the weighted sum of the waiting and the traveling times, minimizing the maximum idle time pertinent to each facility, and minimizing the total cost associated with the opened facilities. According to these researchers, minimizing the average idle time does not necessarily lead to the minimization of possible idle time of all facilities. Therefore, minimizing the maximum idle time decreases the possibility of idle time of all facilities. To improve the service delivery quality level, they used a coefficient along with the service rate of each service provider in capacity constraint. In a research, Hajipour et al. (2014) evaluated a server locationallocation problem and entered appropriate levels of redundancy or reliability to an M/M/m queuing system. They solved the model to maximize the system’s reliability, minimize the system’s costs, and minimize the waiting time with the help of a MOVDO, multilevel genetic algorithm (GA) and refrigeration simulation. According to the results, the MOVDO had a higher performance compared to the other two algorithms.
In a study, Araz et al. (2014) presented a middle location model with a certain number of facilities in the field of urban public health emergency services in Arizona, the United States. They evaluated the distribution location problem for providing medical services related to anthrax and allocating demand points to these centers by considering human resource allocation to each distribution center. Each distribution center in the network was considered an M/G/c queuing system. The objective functions of minimizing the total mean waiting time in the system and the total customer travel time were considered the most important system performance measurement criteria in the model. Basciftci et al. (2021) presented a twostage mathematical model for facility location under uncertainty. They used a robust optimization approach and random programming, and the results were indicative of a significant improvement in the robust optimization approach. Saif and Delage (2021) focused on the optimization of facility location problems using the datadriven robust optimization method. In addition, column generation was applied to solve the model. Silva et al. (2021) presented a heuristic approach to solve a locationallocation problem with different capacity modules. In this research, programming was considered periodically, and demand was introduced as a dynamic parameter. It is notable that GA was exploited to solve the problem.
3. MATERIALS AND METHODS
In this study, the locationallocation problem is hierarchical. Each level presents its certain services and facilities existing in each level provide a similar type of service. In other words, facilities in each level are similar to each other. Customers or applicants for service must go through all levels to complete their service and must receive service at each level of one of the facilities. The customers of these services are located in specific and predetermined places, which are divided into demand points. Customer demand has a potential behavior that follows the Poisson distribution with a certain parameter. A certain number of potential points exist for locating the facilities of each level, and the objective of the problem is to choose an appropriate location among potential locations for the establishment of a proper number of facilities at each level. In addition, customer allocation to each facility of each level is carried out in two modes of priority system and shift system in serving (Barzamini and Ghassemian, 2019). The phenomenon of congestion and queue formation is also considered in the present research. Given that each customer point follows the Poisson distribution, the demands allocated to each facility, which is the total number of Poisson distributions, have Poisson distribution. Furthermore, each facility has a server with an exponential service time. Therefore, each facility is an M/M/1 queue system. However, the facilities have no capacity constraint, and the system productivity coefficient is considered below one in the long term. The conceptual model of the problem is presented in Figure 1.
A multiobjective nonlinear integer programming model is presented below, the first and second objective functions of which is minimizing the total waiting time of customers in the facility queue and minimizing the maximum idle time of facilities, respectively. Evidently, the objectives are contradictory since the system’s waiting time improves with decreasing the possibility of idle time. In addition, idle time decreases with increasing the waiting time in the system. These objective functions create a balance between the objectives of the customers and the owners of the system since the first objective function works for the benefit of the customer and the second objective function works for the benefit of the system and the owners.
3.1 Hypotheses

 Each customer receives service from each level only once and it is not possible to return to the past station and withdraw from receiving the service.

 Each facility is considered in the format of the M/M/1 queue model.

 The interval between the entrance of customers and service provision times are independent of each other.

 The productivity coefficient of each facility is below one.

 The server is considered to be immobile and customers refer to the facilities to receive services.

 All demands are met completely and there is no missing demand.

 In the priority mode, customers with a higher priority are placed at the beginning of the facility serving queue to receive services sooner than others. On the other hand, service is provided based on the FIFO policy in the shift system.
3.2 Model Indices and Sets

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

j. j’ : Set of potential facility points at each level j, j’=1, …, N

1 : Facility level index 1=1,…,K
3.3 Model Parameters

D_{i}: Average demand of the ith customer

μ_{1}: Rate of serving in each facility at the lth level
3.4 Model Decision Variables

y_{ij}=1, if the jth facility of the lth level is established; otherwise, y_{jl}=0.

x_{ijl}=1, if the ith customer is allocated to the jth facility of the lth level; otherwise, x_{ijl}=0

z_{1}: Minimization of the total customer waiting time in the facility queue

z_{2}: Minimization of maximum facility idle time

λ_{jl}: Average rate of entrance into the jth facility of the lth level

wq_{jl}: Customer waiting time in the queue of the jth facility of the lth level

${\text{wq}}_{\text{jl}}^{(\text{i})}$: Customer waiting time with ith priority in the queue of the jth facility of the lth level
3.5. Mathematical Model of the Problem
Equation 1 shows the minimization of the total customer waiting time in the facility’s queue. Equation 2 indicates the minimization of maximum facility idle time. Obviously, there is a contradiction between the two objective functions above since, while the total waiting time of customers in the queue decreases with increasing the service rate of the facility, the possibility of facility idle time increases and vice versa. Therefore, facilities must be established in a way that balance would be created between the mentioned objectives. Equation 3 calculates the mean rate of entrance into each facility, whereas Equation 4 guarantees the condition of longterm system stability. Equation 5 calculates the probability of facility idle time in the long term while Equation 6 guarantees that a minimum of one facility is established at each level. Equation 7 ensures that customers must be allocated to the established facilities, and a minimum of one customer and maximum of the entire customers must be allocated to each facility. Equation 8 guarantees that each customer must receive service from only one of the facilities at each level. Equation 9 guarantees that customers must pass all levels. Finally, Equation 10 shows the range and type of model decision variables.
3.6 Analysis of Customer Waiting Time
3.6.1 Customer Waiting Time in Priority System
Since the priority system is mathematically tricky and complicated, the obtained results related to different models of the system are limited. The model considered in the present study is such that the customers enter the facilities of different levels based on the Poisson process with the $(\text{\lambda})$ parameter, and the service delivery rate is equal to $(\text{\mu})$ for all customers with different priorities. In this model, the priority system is established and customers have 1, 2, 3, …, M priorities. Customers 1 are assumed to have the highest priority, followed by customers 2 and others. Ultimately, the priority of (M1) customers is higher than customers (M). In addition, the number of servers in each facility is equal to one person. According to the premises, the waiting time of the customer with the kth priority in the queue is calculated using the equation below (Basciftci et al., 2021):
3.6.2 Customer Waiting Time in the Shift System
In the shift system, the desired model is such that the customers enter the facilities of various levels based on the Poisson process with (λ) parameter, and the service delivery rate is considered similar and equal to (μ) for all customers. In this model, the shift system is established and customers have no priority over each other. In addition, the service delivery process is based on the FIFO policy. Moreover, the number of servers in each facility is equal to one person. Therefore, the customer waiting time can be calculated by the following equation using the mentioned premises (Roozitalab, 2022;Basciftci et al., 2021):
3.7 Proposed NSGAII Algorithm
The modified version of NSGA and GAbased, the NSGAII is one of the conventional methods used to solve multiobjective problems. Figure 2 shows the flowchart of the NSGAII algorithm.
The NSGAII has the following steps (Goli et al., 2021):

Step 1: a random initial solution with the size of $i=1,\mathrm{...},pop$ is generated, k number of NSGAII algorithm iteration is set to 1 (k=1).

Step 2: particles are sorted based on dominance and are divided in fronts. The lower the number of fronts, the higher the number of particles dominated by those existing in the front. Accordingly, the following steps are carried out for each particle, such as P particle:

Step 21: SP is considered as the total number of population members dominated by P particle and its value is considered to be zero (SP=0).

Step 22: N_{p} is considered as the frequency of dominance of P particle, compared to other particles, and its value is considered to be zero (N_{p}=0).

Step 23: the following steps are carried out for each population member n=1,…,pop_size, similar to q.

Step 231: q will be added to S_{P} if P particle is able to dominate q particle.

Step 232: a unit will be added to N_{p} if q particle is able to dominate P particle.

Step 3: if N_{p}=Ø following the assessment of all particles, then it could be concluded that P is dominated by none of the other particles. Therefore, P is added to f_{1}. In other words, ${f}_{1}={f}_{1}\cup \left\{P\right\}$.

Step 4: all of the following steps are continued until the number of particles existing in i front is not equal to zero (${f}_{i}\ne 0$).

Step 41: 1 the set of particles considered in i+1 front is recognized as Q and is considered equal to zero (Q=0). Afterwards, the following stages are carried out for each P particle existing in fi.

Step 42: the following steps are taken for each particle (such as q) that exist in S_{P} in f_{i}. (It is notable that S_{P} is a set of particles dominated by P particle in the previous stage).

Step 421: a unit is subtracted from N_{q}, which shows the frequency of dominance of the q particle.

Step 422: 2 if N_{q}=Ø shows that q particle is in f_{i+1}, then Q must be moved with q (${f}_{i+1}=Q\cup \left\{q\right\}$).

Step 43: a unit is added to i (i=i+1)

Step 5: after particle fronting based on the degree of dominance of other particles, a number of particles are selected to create the next generations. The feasibility of the solutions is determined based on the following two criteria:

Step 51: 1 rank priority: in this priority, those solutions that have lower ranks or fronts are selected because the particles of these fronts can dominate most particles.

Step 52: in some cases, the two selected particles might have the same rank. In other words, they might be in the same front. In this case, a criterion known as CD is applied, which is explained below.

Step 521: for each n_{i} particle, f_{i} is considered as the number of particles existing in the particle.

Step 522: the distance between particles in fronts is called d1, and the distance of all particles from each other is considered to be zero (${f}_{i}({d}_{i})=0$).

Step 523: each of the objective functions of the problem (such as m) is considered for each particle (such as j) in f_{i}, and the following steps are taken:

Step 5231: in f_{i}, all particles are sorted based on the m objective function. in other words, the particles existing in f_{i} are sorted based on their objective functions separately.

Step 5232: following sorting particles in f_{i} based on the m objective function, I(d) crowing distance of the first and last particles are considered equal to infinity ($I({d}_{1})=I({d}_{n})=\infty $)
This is mainly due to the lack of existence of another particle near the particles to cover them. The $I({d}_{k})$ crowding distance is determined for particles 2 to n1 based on equations 16 and 17.
In Equation 16, $I({d}_{k}).m$ means the crowding distance from the m objective function. In order to calculate the total crowding distance, $I({d}_{k})$ must be calculated and summed up for all objective functions, as shown in Equation 17.

Step 524: following the calculation of crowding distance (CD), those particles that have a higher CD will be selected.

Step 6: after particle selection in the previous stage, a pool is created called the selected population. Afterwards, genetic operators are used to create the population of offspring. The genetic operators used in this article are crossover and mutation operators.

Step 7: after determining the population of offspring obtained from P_{t} genetic operators, this population is integrated with the main Q_{t} population. Each pool has n capacity, and the particles that are integrated with each other must be eliminated. To this end, the following steps are taken to achieve n capacity.

Step 71: particle fronting is performed based on the method explained in Step 2.

Step 72: CD of each particle in fronts is determined.

Step 73: the process starts from f_{1}, where particles are selected based on their CD and are added to the new population pool (K+1). This step continues until the capacity of the new pool reaches n.

Step 8: we move on to Step 2 after forming (K+1) population, and the stages are repeated until reaching the determined size.
Figure 3 presents a schematic illustration of NSGAII algorithm.
4. COMPUTATIONAL RESULTS
In this section, a number of numerical examples are solved and analyzed to evaluate the model. According to SheraliNordai, the multifacility locationallocation problem (facilities that can be used in each point of Euclidean space) with certain parameters is considered an NPhard problem. Therefore, the rational solution for such problems at different levels is to use heuristic and metaheuristic approaches. Therefore, the NSGA is used to solve the proposed model. Given the lack of finding similar studies in the literature, data are generated randomly to solve the model. In addition, since most locationallocation articles have created the required data with uniform distribution, the same approach is applied in the present study to create the desired data. Table 1 shows the random parameters of the problem.
4.1 The Second Version of the NonDominated Sorting Multiobjective GA
This section uses a nondominated sorting GA (NSGA) to solve the multiobjective problems. In the second version of NSGA (NSGAII), which Deb et al. (2000) introduced, the crowding distance approach replaces the sharing function, and an elitism approach is used to provide diversity in Pareto solutions. The NSGAII is proposed to solve the model at different dimensions.
4.1.1 Solution Display Structure
In order to increase the feasibility of the solutions and meet more constraints, we use a specific method in the present research. The structure of the chromosomes showing the problem’s solution is in the form of the following matrix:

The number of rows in this matrix equals the number of customers.

The number of columns in this matrix equals the number of service delivery levels.

Each entry of the matrix is a random number between one and the total number of facilities at that level.
The structure of the mentioned chromosome is shown in Equation 18.
In this chromosome, rows are considered in the number of customers (M) and columns are considered in the number of different serving levels (k). In addition, each entry of the matrix is a random number between one and the number of facilities at each level (J). To decode the chromosomes, the variables related to each customer that is allocated to each facility are considered to be one. By this allocation, the facility is established and the related variable is considered to equate one.
4.1.2 Controllability of Solutions
Given that this way of displaying chromosomes estimates most of the constraints of the model, the solutions generated in each iteration of the algorithm may not necessarily be feasible. This type of chromosome display guarantees the feasibility of constraints 69. The penalty function method is one of the common methods in dealing with constrained optimization problems, which turns these problems into unconstrained problems. Given the fact that the model presented in this research is of multiobjective type, the penalty function is added to both objective functions. The penalty method used in this study is as follows:
Where U is the large positive value, g(x) shows the desired constraint, and p(x) expresses the penalty allocated to the feasible chromosome. It is worth mentioning that the expression above is designed for constraints that are in the form of $\text{g}(\text{x})\le \text{b}$.
4.1.3 Solution Sorting
The solutions obtained in each population must be converted into a local Pareto front using one of the sorting algorithms. In this regard, the nondominated sorting approach is exploited in the present study (Coello et al., 2007).
4.1.4 Genetic Operators
In this section, the random method is used to generate the members of the parent population. In addition, the substitution mutation and the substitution mutation and the singlepoint crossover operator are applied to create the population of the offspring.
4.2 Tuning the Parameters of Multiobjective NSGA
The solutions obtained from the metaheuristic approaches are sensitive to the values of the parameters selected for the algorithm (Coello et al., 2007). In this study, the optimal values of the parameters of the multiobjective NSGA are tuned by the Taguchi method for the example samples at three levels. The values proposed for the parameters are considered based on previous studies and the trialanderror method according to Table 2.
According to the Taguchi standard table, two L9 and L27 levels can be used by considering four threelevel factors (Qazani et al., 2021). In the present study, the L9 level is used due to fewer calculations. Ultimately, the optimal values of the algorithm parameters are tuned by Minitab software, as shown in Table 3. In addition, the values of the problem’s objective functions are solved for the shift and priority systems at different levels using the multiobjective GA, the results of which are presented in Tables 4 and 5, respectively.
According to Figures 4 and 5, it could be concluded that that establishing a queue priority system not only reduces the total waiting time in the queue for each customer but also reduces the total waiting time of all customers in the queue.
5. CONCLUSION
The crowded facility hierarchical locationallocation problem is one of the most widely used location problems that have attracted the attention of many researchers of the field in the past two decades. The present study compared the effect of the queue priority system on facility service provision to the impact of a shift system in this regard. Since the problem was of NPhard type, the second version of a multiobjective NSGA was used to evaluate the efficiency of the model. In addition, the Taguchi design was exploited to tune the parameters. According to the results, the priority system must be used in the design of hierarchical systems and facilities if the goal is to reduce the waiting time of a specific class of customers. In other words, that class must be prioritized in order to decrease the mean waiting time in the queue of the entire system. Considering the priority mode in the queue and service provision process in the system means that the priority system is established, and customers with a higher priority are placed at the beginning of the queue. After entering the facility, if the customer receiving service has a lower priority, their service is halffinished. The server will provide service to the customer with a higher priority, which can be an attractive area for future studies.