Modelling and solving a congested facility location problem considering systems' and customers' objectives

This paper presents a model for the location of facilities subject to congestion. Motivated by applications to locating servers in communication networks, bank branches, automatic teller machines (ATMs), and police services centres in order to ease the access of customers to service centres and to reduce total costs of both customers and service providers, the given model in this research is proposed for situations, in which immobile service facilities are congested by stochastic demand, which originates from nearby customer locations. We aim to locating an optimal number of single server facilities experiencing M/M/1 queuing policy and assigning a set of demand nodes to them in such a way as to optimise three objective functions: cost, time, and quality. Two heuristics in order to find near optimal solutions are proposed and compared using statistical methods in order to determine the best heuristic based on designing a set of numerical examples.


Introduction
We address a congested location-allocation model including a number of service providing firms and stores operating in a region; the stores can be considered as customer nodes to patronise the firms. The number and location of customer nodes (demand points) are known; each customer generates requests for service that constitute an independent Poisson distribution. There are a number of potential facility sites from which a subset of locations is to be selected in order to locate a maximum number of single-server facilities at the selected locations. Each server whose service time is randomly distributed is assumed to be immobile. The customer nodes should be served by the most accessible facility, which corresponds to the closest facility so that the utilisation rate of each server can not exceed a predetermined value.
The current research in this paper is an extension to Wang et al. (2004). The addressed paper gives several models for the location of facilities subject to congestion. In other words, there are three models in Wang et al. (2004) including: a location model from the service provider's perspective with unlimited and limited number of facilities b from the customer's perspective considering no waiting cost and considering waiting cost in the objective function c from customer's and system's perspective with considering no waiting cost and waiting cost in the objective function.
In the given models and most of the other models, the waiting cost is included with other costs of the system as a single objective function. We have developed a three-objective model in which waiting time is considered in a separate objective function. As mentioned earlier, Wang et al. (2004) propose a model with a single objective function optimising cost and waiting time; they present a simple heuristic in order to solve the problem and declare that it can be improved by combining with other heuristics such as tabu search (TS) or Lagrangian relaxation. We propose a model with three independent objectives including maximisation of total quality of the opened facility sites expressed as triangular fuzzy numbers, minimisation of customers' accessing cost to the opened facility sites and total cost of setting up and operating facilities, and minimisation of the sum of average waiting times of customers in different facility sites.
Since the problem belongs to nonlinear integer programming (NLIP) problems, using the common solution methodologies would be time-consuming, specially, as the problem size increases. This is why we propose some heuristics to solve the addressed problem. There are several examples on using meta-heuristics for the nonlinear problems. Lieckens and Vandaele (2007) proposed a genetic algorithm-based heuristic for solving an NLIP problem in designing a logistics network with stochastic lead times. Iannoni et al. (2008) adapted the well-known hypercube model in order to analyse emergency medical systems on highways.

Literature review
Optimisation models for location with (nearly) constant demand are derived from the location set covering problem given by Toregas et al. (1971), the p median given by Hakimi (1964) and Revelle and Swain (1970), and the maximal covering location problem presented by Church and Revelle (1974). Real situations very often have demand for their services that is both variable and random in nature. Although the facility may be able to cope with average demand, there are times of heavy demand when it can not cope; such a facility is called congested (Boffey et al., 2007).
The maximum expected covering location problem (MEXCLP) maximises the expected coverage by a free server (Daskin, 1983). Larson (1974) combined the concept of queuing model with facility location problem. Berman et al. (1985) developed a heuristic algorithm to optimally locate one server on a congested network. The probabilistic location set covering problem by Revelle and Hogan (1988) forced all demands to be covered with a pre-specified reliability. Revelle and Hogan (1989) presented a model, which maximises the population covered with a pre-specified reliability.
When dealing with designing service networks, such as health and emergency services, banking, or distributed ticket selling services, the location of service centres has a strong influence on the congestion at each facility and, consequently, on the quality of service. Marianov and Serra (1998) formulated several maximal coverage models with one or more servers per service centre and developed heuristics to solve the models for a 30-node service network. The main constraint of the models is that nobody should stand at queues for a time longer than a given time limit. An extension of the previous research is Marianov and Serra (2001), which seeks to cover all population and includes server allocation to the facilities. Based on the number of servers allocated to each facility being constant or variable, two different models are proposed and some heuristics with good performance on a 55-node service network are developed and tested.
In a recent research, Berman and Drezner (2007) introduce multi-server location problem in a stochastic environment. The objective is to minimise the sum of travel time and average waiting time for all servers. In Wang et al. (2004) at most, one server can be located at any potential location; thus, the system is modelled as an M/M/1 queue rather than an M/M/k queue. A similar model with additional constraints on demand lost because of insufficient coverage or congestion, is studied by Berman et al. (2006). Aboolian et al. (2008) introduced the multi-server centre location problem; in this problem, a given number of servers are to be located at nodes of a network. Demand for services is located at each node, and a subset of nodes is to be chosen to locate one or more servers in each. Each customer selects the closest server. The objective is to minimise the maximum time spent by any customer, including travel time and waiting time at the server sites. Applications for such a model include the location of medical facilities in which the number of medical care personnel must be determined at each location, location of post offices, and location of bank branches (Aboolian et al., 2008). Beraldi and Bruni (2009) formulated and solved a probabilistic model in order to determine the optimal location of facilities in congested emergency systems.
The research in this paper is closer to that of Wang et al. (2004) and Aboolian et al. (2008). Since we have used fuzzy logic in order to quantify the quality of the facility sites, we give a brief review of the corresponding literature regarding facility location. Fuzzy logic was initially introduced by Zadeh (1965). This concept was applied to location models by Canos et al. (1999). The major research in this area was suggested by Shavandi and Mahlooji (2006) presenting a fuzzy location-allocation model for a congested systems such as health and emergency services, public safety, fire fighting and so on. In the addressed paper, they utilise fuzzy logic to develop a queuing maximal covering location-allocation model considering fuzzified queuing parameters as well as fuzzified constraints. A genetic algorithm is developed to solve and test the model using up to 50-node networks.
Şakir and Tarık (2009) present a fuzzy clustering-based hybrid method for a multi-facility location problem. The capacity of each facility is assumed to be unlimited. Customers are grouped by spherical and elliptical fuzzy cluster analysis methods with respect to their geographical locations; then, facilities are located at the proposed cluster centres. Finally, each cluster is solved as a single facility location problem. Liu and Xu (2011) introduced the concept of hybrid variable and proposed a mixedinteger programming model for random fuzzy facility location-allocation problem. The model is reduced to a deterministic model by expected value and chance constraint programming technique. A priority-based genetic algorithm was proposed in order to solve the model. Sajjadi and Cheraghi (2011) presented a model integrating inventory control, as a tactical decision, with location-routing problem. The presented model considers the multi-product network under the fixed order interval inventory policy with stochastic demands for products. This paper also presented a solution approach based on simulated annealing (SA).
Pasandideh and Niaki (2012) developed a bi-objective facility location problem within M/M/1 queuing framework considering the p-median problem. They utilised the desirability function technique and genetic algorithm to solve the problem. Konur and Geunes (2012) developed formulation on location decisions for competitive firms serving a set of markets. Firms come across firm-specific transportation, congestion, and location costs. Market price is linear and decreasing with the amount shipped to the markets by all firms. Pasandideh et al. (2013) proposed a three-objective model with batch arrivals to service centres. They convert this model into a single-objective one using LP-metric method and then solve it using GA and SA algorithms. Lakshmi and Iyer (2013) gave a review regarding the contributions and applications of queuing theory in the field of healthcare management problems. This review provides classification of healthcare problems, which are modelled using queuing models. Jouzdani et al. (2013) proposed dynamic dairy facility location and supply chain planning in such a way as to minimise the costs of facility location, traffic congestion and transportation of raw/processed milk and dairy products under demand uncertainty. Jabal-Ameli et al. (2013) formulated cell formation and layout design as different mathematical models, which differed with respect to layout constraints and backward movements. A multi-objective scatter search is proposed to solve the integrated cell formation and layout design problem. They introduce a new solution approach. The solution approach is compared with the epsilon-constraint method and a well-known multi-objective evolutionary algorithm, non-dominated sorting genetic algorithm II (NSGA-II). The proposed algorithm performs better than NSGA-II. Rahmaniani et al. (2013) tackled the uncapacitated multiple allocation p-hub location problem under an uncertain environment developing the mathematical formulation of the considered problem with uncertainty in flows and distances parameters. Furthermore, based on the variable neighbourhood search, a heuristic is proposed to solve the problem. Wang et al. (2014) proposed queuing problems of bulk arrival and service with the balking and reneging behaviour of customers. This research formulates queues of this type using compound Poisson processes and gives some key probabilistic measures. Vidyarthi and Jayaswal (2014) proposed a nonlinear model for location-allocation problem considering the effects of congestions and queuing delays. The problem searches for simultaneously locating service facilities, equipping them with appropriate capacities, and allocating user demand to these facilities in such a way as to minimise the total system costs. Alenezi and Darwish (2014) integrated location model with risk pooling and mode selection problems in a supply chain consisting of one supplier, distribution centres, transient points, and retailers who come across random demand which follows a certain probability distribution. The model determines facility location decisions, inventory decisions, and intermodal transportation decisions. Lagrangian relaxation approach is used to develop lower and upper bounds on the optimal solution. This paper is structured as follows: in Section 3, notations used to model the problem are defined. Section 4 gives nomenclature and assumptions. In Section 5, the model is stated as a mathematical programming including objective functions and constraints. In Section 6, two SA and TS-based heuristics are proposed in order to solve the problem. In Section 7, a few numerical examples are designed and solved using the proposed heuristic. Finally, the conclusions resulted from the research are presented and further research in this area is recommended in Section 8.

Notation
The required notations to model the addressed problem are as follows: The decision variables are as follows: 1 if a facility is opened at site , 0 otherwise, Given the above notations, we can conclude . • There are a number of customers with a constant demand rate.
• The number of potential facility sites and their location is known.
• The maximum number of facilities which can be opened is known.
• The quality of each facility site is representing by a triangular fuzzy number.
• Access cost of customer to each facility site is known.
• Each customer is assigned to exactly on facility.
• Each customer is assigned to facility with lowest distance.

Problem formulation
The problem is to locate an optimal number of single server facilities, following M/M/1 queuing system policy, at some facility sites and to allocate a set of demand nodes to them in such a way as to satisfy three objective functions in terms of cost, time and quality. The model can be stated as follows: Subject to: The objective function Z 1 in (1) represents the total quality of the opened facility sites, which should be maximised. The objective function Z 2 in (2) consists of two sections; the first section represents customers' accessing cost to the corresponding facility site while the second represents the total cost of setting up and operating facilities per time unit. The objective function Z 3 in (3) represents the sum of average waiting times of the customers in different facility sites. Constraint (4) ensures that each customer node is allocated to exactly one facility site. Constraint (5) guarantees that each customer is assigned to an opened facility site. Constraint (6) declares that the number of opened facility sites should not exceed the predetermined maximum number . P Constraint (7) ensures that each customer at node i patronises the facility site at lowest accessing cost; K is a very big positive number. Constraint (8) ensures that the utilisation rate of each server at the corresponding opened facility site does not exceed the upper bound. Constraint (9) indicates that decision variables of the model are binary. Lai and Huwng (1992)'s proposed technique is used to convert the fuzzy objective function Z 1 into crisp one assuming that the triangular fuzzy number j q is represented as (s j , l j , u j ). According to the aforementioned technique, instead of maximising fuzzy value of Z 1 , a set of three objective functions is defined in such a way as to maximise the central and right values (i.e., l j y j and (u j -l j )y j ), and to minimise the left value (i.e., (l j -s j )y j ). As a result, (1) can be replaced with (10) to (12).
Now, there is a mathematical model with five objective functions of 1 1 , Z 2 1 , Z 3 1 , Z Z 2 , and Z 3 . Initially, we convert the model into a single objective one and then apply a few heuristics in order to solve it. We use desirability function method (DFM) as Pasandideh and Niaki (2006) in which such a method is used for solving a multi-response simulation optimisation problem; furthermore, Pasandideh and Niaki (2012) solve a congested location problem using a genetic algorithm within desirability function framework. According to DFM, after converting all objective functions to maximisation format, an upper and a lower bound is considered for each objective Z i (i.e., Lower i ≤ Z i ≤ Upper i ), and then a total objective function D is considered to be maximised as in (13) in which D 1 , …, D k are defined for each objective as in (14).
1, 2, , The behaviour of total objective function D depends on different values of S in (14); it increases linearly for s = 1, it is convex for s < 1, and concave for s > 1. We have assumed s = 1. To simply deal with the minimisation objective functions, we can define d l as in (15).
We let x ij and y j , , i M j N ∀ ∈ ∈ equal to one or zero in order to determine the upper and lower bounds for each objective function except for the last objective function Z 3 . For this objective function, we can easily find a lower bound equal to zero, however, the upper bound can be found as in (16) (17) to (22) subject to (4) to (9).

Proposed heuristics
Since the proposed model in this research belongs to NLIP, we have proposed two SA and TS-based heuristics to solve it. Closer to this model is that of Wang et al. (2004), which is proved to be NP-hard by Garey and Johnson (1979).

SA-based heuristic
The main idea of SA has been derived from metals annealing process, which was initially developed by Metropolis et al. (1953) and then was used by Kirkpatrick et al. (1983) for solving combinatorial optimisation problems. Anthony et al. (1994) studied the application of SA for operations research problems. SA is of no memory structure. There are a few important parameters in this heuristic. We have assumed the initial temperature (T 0 ) to be selected in such a way that the system does not get stuck on a subset of solutions at the beginning of the process. A proposed heuristic in order to find initial temperature is given by Heckmann and Lengauer (1995) in which the problem is solved for 20 times and the best objective function values are obtained; then, the standard deviation of the generated values for σ is calculated and T 0 is obtained from (23).
ρ is a random number and takes a value from 0.5 to 0.95; we have selected 0.95 for this purpose. Let us consider the following notation: The general outline of the proposed SA-based heuristic is as follows: The steps of the algorithm are as follow: Step 1 Initialising the annealing parameters including T 0 as given by equation (23), TDN, α, RTN; set P = 0, S = { }, t = 0, and il = 1.
Step 2 Generate a feasible solution not violating the constraints of the model in such a way that a site is opened randomly and all customers are allocated to this site; then, set P = 1 and obtain objective function value of the current solution (DS) from (17); let DS * = DS.
Step 3 Schedule annealing: Step 3.1 Generate a neighbourhood as follows: Initially, generate a random number r.
If r > P r -P facilities are opened randomly.
If r < P P -r facilities are closed randomly.
If r = P a facility is closed and another is opened randomly.
Finally, update S and calculate DS.
Let DS * = DS; Else, the new solution is rejected.
Step 3.3 If il ≤ RTN, then inner loop is continued; set il = il + 1 and go to step 3.1; Else, inner loop is terminated; set il = 1 and go to step 3.4.
Step 4 Print the best solution (DS * ) and stop.
The flowchart of the algorithm is given as in Figure 1.

TS-based heuristic
Initially, TS was introduced by Glover (1989) for solving vendor problem. High efficiency of TS has made it to be considered as a solution for other problems such as location problems. Rolland et al. (1997) used TS for solving p-server location problem. Filho and Galvao (1998), Gosh (2003), and Sun (2006) applied it to solve classic problem of server location. Arostegui et al. (2006) indicated that TS was more effective than the other meta-heuristic algorithms in solving location problems.
The most important point while applying TS is the definition of tabu list and determination of its length. Since tabu list is restricted, it is usually updated by the policy of first in first out (FIFO) or last in first out (LIFO).
The tabu lists are historical in nature and form the TS memory. Maintaining a list of solutions which must be avoided is the major role of this list.
Let us consider the following notations: The general outline of the proposed TS-based heuristic is as follows: Step 1 Let set of open facilities S be empty; TL = {}; Ins = 0.
Step 2 Generate a feasible solution that satisfies all constraints and obtain objective function value of the current solution (DT) from (17); let DT * = DT; update TL.
Step 3 Generate n P − neighbours at each iteration as the following: • Generate g 1 ; set G = {g 1 }; If 1 , G TL ⊂ ignore neighbourhood generation; otherwise, the new facility (corresponding to g 1 ) is closed if it is open and is opened if it is close; obtain DT 1 . Step 5 If |TL| > L update TL according to FIFO policy (by terminating the additional opened facilities), else continue.
Step 6 If Ins ≥ CSN, then stop (which means we have performed CSN successive swaps that have not improve the objective function value); otherwise, go to step 3.
Step 7 Print the best solution (DT * ) and stop.
The flowchart of the algorithm is given as in Figure 2. Example 1: Consider a network with four customer nodes (m = 4). Initially, there are five potential facility sites (n = 5), but the maximum number of facilities ( ), P which can be opened is assumed to be three. Average service rate at each facility site j (μ j ) is assumed to be 8, the upper bound on the server utilisation rate at each open facility site (r) is assumed to be 0.7, demand rate of service required by customers at nodes 1, 2, 3, and 4 is considered 2, 1, 3, and 1, respectively. Setup costs of facility sites per time unit (f j ) are assumed to be 5, 3, 4, 3, and 2 for the sites 1, 2, 3, 4, and 5, while the operating costs of giving service at the sites (b j ) are assumed 1, 2, 1, 3, and 2, respectively. The accessing cost matrix between customer nodes (rows) and facility sites (columns) is considered as follows: 1 2 1 3 2 2 3 3 1 4 4 1 2 3 3 We have used a few triangular fuzzy numbers as provided in Table 1 in order to specify the quality priorities of the facility sites 1 to 5. Supposing that we have used fuzzy pair wise comparisons, the values of the quality of the facility sites turn out to be Source: Atesh et al. (2006) Let TDN = 60, α = 0.95 and RTN = 40. Suppose we are at the initial temperature. If we only open facility 4, S = {4} with demand rate γ 4 = λ 1 + λ 2 + λ 3 + λ 4 = 7 since all customers are allocated to facility 4; the objective function value is calculated, DS * = 0.3634. In the next step, assume that two new facility sites should be opened. Consider that facility sites 1, 2 are selected; so, S = {1, 2, 4} and demand rates are turned out to be γ 1 = λ 1 = 2, γ 2 = λ 3 = 3, and γ 4 = λ 2 + λ 4 = 2 with DS = 0.5785; therefore DS * = DS = 0.5785. This process continues until il = 40. In the next step, we let t = t +1 and T t is updated; then, this process continues until t > TDN.
The final values of decision variables y j 's and x ij 's turn out as: (1 1 1 0 0) j y = 0 0 1 0 0 Which means three facilities are opened at sites 1, 2, and 3; customer nodes 1 and 4 are assigned to facility at site 1, customer node 2 is assigned to facility at site 1, and customer node 3 is assigned to facility at site 2. The value of the objective function given by (17) turns out to be 0.6140. The details of the value are as follows: In the rest of Section 7, we design and solve 16 numerical examples of various sizes in order to verify the model accuracy and to measure the accuracy of the proposed heuristics. Table 2 gives the values of customers nodes (m), potential facility sites (n), average service rate at each facility site j (μ j ), and upper bound on the server utilisation rate (r). We have considered the values of μ j 's to be identical for all j's and also to be greater than sum of customers' demand rates. Such an assumption guarantees that the system can serve all demands if only one facility site is opened. The upper bound on the server utilisation rate (r) is generated randomly between 0.5 and 1. The values of other parameters are generated randomly; we assume demand rates to be randomly generated from Alenezi and Darwish (2014) and Boffey et al. (2010), setup costs of a facility sites per time unit to be randomly generated from Atesh et al. (2006) and Church and Revelle (1974), and operating costs of giving service at a site to be randomly generated from Aboolian et al. (2008) and Berman et al. (2006). The accessing cost matrix values between customer nodes and facility sites are assumed to be generated from [1,18]. The value of the quality of a facility site, which is a triangular fuzzy numbers, is also generated randomly; initially, three random numbers are generated between 0 and 1; the first number is considered as the left number, the second number plus one is considered as the middle number, and the third number plus 2 is considered as the right number.
The results from the SA and TS-based heuristics are reported in Tables 2 and 3 respectively. We have also solved each example using Lingo solver for which Table 3 gives the results.
In order to explain how we have tuned the parameters of the heuristics, we give the parameters of the heuristics. For the SA, the Parameters include TDN, α, and RTN. For the TS, the parameters are CSN and LT.
We have used the well-known response surface method (RSM) in order to tune and determine the optimal values of the parameters. For the SA-based heuristic we have considered the initial values of the parameters as 30 ≤ TDN ≤ 90, 0.89 ≤ α ≤ 0.99, 10 ≤ RTN ≤ 70. The optimal values of the parameters for each numerical example turn out to be as given in Table 2. For an instance, we just give the estimated fitness function for numerical Example 1, which is represented by f (TDN, α, RTN) as in (24). We find optimal values of the parameters maximising f (TDN, α, RTN) For the TS-based heuristic, we have considered the initial values of the parameters as 40 ≤ CSN ≤ 160 and 2 ≤ L ≤ 16. The optimal values of the parameters for each numerical example turn out to be as given in Table 3.     Solving the numerical examples from various sizes, we can conclude that for smaller sizes where facilities considered for opening servers are of few values, the two given heuristics and Ling are of the same result, such as 3. For larger sizes, due to the existence of a lot of combinations for solutions and considering the fact that the obtained solution is local optimum in most cases, the generated results are different, however one method can not be considered superior than the others. Furthermore, for the heuristics, SA has worse objective function values than TS, except for Examples 1 and 2. Using analysis of variance (ANOVA), where hypothesis '0' is equality of the results and hypothesis '1' representing significant difference among the results the results are as in Figure 3. Since the obtained value for P is 0.429, hypothesis zero is accepted which means there is no significant difference. The processing time of TS-based heuristic is less than that of SA; Lingo solver had higher processing time for the larger sizes numerical examples.

Conclusions and further research
In this paper, we presented a model for the location of facilities subject to congestion. The model can be applied for situations in which immobile service facilities are congested by stochastic demand originating from nearby customer locations. The model found optimal number of single server facilities, following M/M/1 queuing system policy, at some potential facility sites and assigned a set of demand nodes to them in such a way as to satisfy three objective functions in terms of cost, time and quality. The cost consists of two parts; the first part represents customers' accessing cost to the corresponding facility site while the second represents the total cost of setting up and operating facilities per time unit. The time represents sum of average waiting times of the customers in different facility sites. The quality represents the total quality of the opened facility sites. Two SA and TS-based heuristics were proposed in order to find near optimal solutions; Statistical analysis using T-test with n = 16 and t 0.05,15 = 1.753, gives the value of 0.63 for t. The results indicated that both heuristics gave solutions near optimal since there was no significant difference among the results obtained from the proposed heuristics and Lingo solver whose solutions were near optimal. The processing time of TS-based heuristic was less than that of SA; Lingo solver had higher processing time for the larger size numerical examples.
As further research, we can recommend considering other queuing policies for serving the customers at the facility sites besides developing other efficient meta-heuristic algorithms such as genetic, ant colony optimisation, particle swarm optimisation and so on in order to solve the addressed problem.