Queueing Maximal Covering Location-Allocation Problem : An Extension with M / G / 1 Queueing Systems

We consider the queueing maximal covering location-allocation problem QM-CLAP with an M/G/1 queueing system. We first formulate the problem as a binary quadratic programming problem and then propose a new solution procedure based on decomposition of the problem into smaller binary quadratic sub-problems. The heuristic procedure GRASP is used to solve the subproblems, as well as the entire model. Some computational results are also presented.


Introduction
An extensively studied and widely used model in location theory is the maximal covering location problem MCLP .It can be stated as follows.Given a set of m demand points and a set of n potential new facility locations in the plane, find the location of p < n new facilities such that the demands covered are maximized.A demand point i is said to be covered by a facility j if the distance or time between i and j is no more than a given fixed value, R. In this simple form, the problem is attributed to Church and ReVelle 1 who have modeled it as a binary integer program.Extensions and applications of this model have been discussed by many including ReVelle 2 , Schilling et al. 3 , Marianov and ReVelle 4 , Brotcorne et al. 5 , Marianov and Serra 6 , and Goldberg 7 .A bibliography for some discrete location problems, including maximal covering location problem, can also be found in ReVelle et al. 8 .Erkut et al. 9 present a computational comparison for five maximal covering models.
A variant of this problem, usually called location in congested systems, takes into account the stochastic nature of demand and response times.Berman et al. 10 are among the first to incorporate such systems into location models.They extend Hakimi's 1-median Advances in Decision Sciences problem by assuming that demands at nodes follow a Poisson process, and the service time has a general distribution, hence, a network with an M/G/1 queueing system.They propose an algorithm to find the optimal location of a server on such networks.Berman and Mandowsky 11 present a location-allocation model for districting a region with congestion.Berman et al. 12 consider the problem of locating a single facility with M/G/k queueing system.Batta and Berman 13 consider locating a single facility on a network with an M/G/k queue.Batta 14 considers the problem of locating a single server on a network operating as an M/G/1 with nonpreemptive priority queue.Marianov and ReVelle 15 consider the queueing probabilistic location set covering problem with M/G/s/s queueing system.They 4 , also use an M/G/s/s queueing system for siting emergency facilities.Marianov and Serra 16 consider probabilistic maximal covering models with constraints on waiting time and for queue length.Jamil et al. 17 consider the problem of locating a single facility with a "center-type" objective in an M/G/1 queueing environment.Silva and Serra 18 consider the MCLP with an M/M/1 queueing system and with different priority levels.They consider both a "directed choice environment" model, where a demand point may be allocated to one center for one priority and to another center for a different priority, and a "user choice environment," where it is assumed that the customer always chooses the closest center.They propose a heuristic procedure to solve this problem.Moghadas and Kakhki 19 consider the QM-CLAP with M/M/k system and side constraints.In this model, k is unknown and some constraints on the number of servers at each center, as well as constraints on the total costs of establishing a center and locating servers, are imposed.Berman and Krass 20 is a thorough review of the models with stochastic demands and congestion at facilities.
In this paper, we consider the maximal covering location problem in a congested system.A single mobile server resides at each center, and demands for service occur in time as a Poisson process.If the server is available, it is immediately dispatched to the demand point.After providing the service, the server returns to its base.If the server is busy, the customer waits in a queue with an M/G/1 system.The objective is to choose the location of at most p service centers and to allocate demand points to those centers so that the population covered is maximized.Each demand point should be serviced by at most one established service center, and the average waiting time at each service center must not exceed a given threshold.We first formulate the problem as a binary quadratic programming problem with linear and quadratic constrains.Then we propose two solution procedures.The first procedure is based on a decomposition of the problem into smaller binary quadratic knapsack-type subproblems.The second method is the metaheuristic procedure GRASP.Finally some computational results are presented.

Problem Formulation
In order to formulate the problem, we use the following notations and variables: -I: the set of all existing demand points incident locations |I| m , -J: the set of all possible locations of new facilities centers |J| n , -a i : population at point i, -N i : the set of points in a pre-specified neighborhood of i; that is, and propose a decomposition approach to solve it.Silva and Serra 18 considered the priority queueing covering location problem PQCLP for an M/M/1 system and proposed a GRASP heuristic procedure to solve their problem.
In this paper, we consider a problem in which requests for services at each demand point occur with a Poisson process.If the server is available at the time of a call, it travels to the demand point to provide on-scene and perhaps off-scene services.A service is completed only when the server returns home.If the server is unavailable, the customer is entered into a queue with an M/G/1 system, with an infinite queue length and a FIFO discipline.If we assume that the arriving calls from a demand point i have a Poisson distribution with intensity f i , then the requests for service at center j is a Poisson process with an intensity λ j , λ j i∈I:j∈N i f i x ij see e.g., Marianov and Serra 16 .i ∈ I : j ∈ N i denotes the set of all demand points i which are in a neighborhood of candidate center j.
For an M/G/1 queueing system, the waiting time is given by the Pollaczek-Khintchine formula Kleinrock 24 : where S j and S 2 j are the first and the second moments of the service time at center j, respectively.Hence, if the stability condition 1 − λ j S j > 0 holds, then constraints 2.5 can be written as Now, if we define T ıj and T 2 ıj as the first and second moments of service times for the customer i at center j, then we have see, e.g., Jamil et al. 17 where h i is the fraction of calls originating from demand point i which is defined as

2.10
Therefore, the problem with an M/G/1 queue can be written as Max 2.2 -2.4 , 2.10 and 2.6 .

2.11
To determine T ıj and T 2 ıj , we use the same assumptions as in Berman et al. 10 ; that is, the average service time for demand point i from service center j consists of travel time to the scene, on the scene service time, travel time back to the center j, and possibly additional off-scene time.T ıj and T 2 ıj are defined as where the parameter β is greater than or equal to one, d i, j is the distance between demand point i and center at j, v is the speed, and Z ı is the average on-scene plus off-scene "nontravel related" service time associated with node i.This problem is NP-hard, and, in general, as the problem is a nonconvex programming problem, even small instances cannot be solved with the existing commercial softwares.A semidefinite programming relaxation for the problem is recently presented by Kakhki and Moghadas 25 .In Section 3 we discuss an approach for decomposing the problem into smaller subproblems and present some solution procedures.

A Solution Approach for Problem P
The basic idea for our proposed algorithm is to decompose the problem into smaller subproblems which is inspired by the solution methods suggested in the literature for the capacitated p-median problem see, e.g., Baldacci et al. 26 and Lorena and Senne 27 .
Suppose that we have determined the location of k − 1 k ≤ p centers and want to determine the location of the kth center from among the n−k 1 remaining candidate centers.If we remove the demand points assigned to the previous k − 1 centers, then constraints 2.2 would be satisfied for any customer.In addition we consider the constraints 2.3 and 2.4 implicitly.Therefore, the problem can be decomposed into smaller binary quadratic subproblems BQSP as follows: where I is the set of remaining demand points, not assigned to the k − 1 previous centers.Now starting with k 1, we can solve n − k 1 binary quadratic subproblems P j for all j not in the set of selected centers, J * .We then find the optimal solution, j * , delete all the demand points assigned to j * from I, and add j * to J * .The procedure is continued until k ≥ p, or all the nodes are covered I φ .This procedure is outlined in Algorithm 1, below.
Solve P j , I solves the binary quadratic subproblem P j for a candidate center j.If a demand point i i ∈ I is assigned to j, then i is added to I * j .Update I removes all demand points assigned to j * from I.
This procedure can be easily extended to solve the more general case of fixed charge facilities; that is, when the objective is Max i∈I:j∈N i a i x ij − j∈J c j y j , where c j is the cost of establishing a center at site j.Notice that problem P is nonlinear and can be decomposed into smaller binary quadratic subproblems.Treating the quadratic terms x ij x tj 's as a single variable and using linearization methods, such as that of Sherali and Adams 28 will considerably increase the dimension of the problem.Instead we propose to solve the subproblem P j using GRASP.In the next section, we discuss this procedure in more detail.

Solving Subproblem P j with GRASP
GRASP Greedy Randomized Adaptive Search Procedure has been developed in late 1980s by Feo and Resende 29, 30 .Since then, it has been successfully applied to many problems, including the maximal covering location problem Resende 31 .A bibliography of different applications is compiled by Festa and Resende 32, 33 .GRASP consists of two main phases: construction and local search.In the construction phase, a solution is built using a greedy function and randomization.In the local search, an optimal solution in the neighborhood of the solution found in the construction phase is obtained.Here we describe the details of the procedure for solving subproblems P j with GRASP.Algorithm 2 shows the GRASP algorithm for solving BQSP P j .
Algorithm 3 shows the details of the preprocess and CGRS-LS procedures.Here W I * j is defined as W I * j i∈I * j a i .The procedure starts with k 1, key 1.The best solution is stored in BestSolution.The steps of the algorithm are repeated for Maxitr iterations.If a new solution is constructed in the construction phase, then key is set equal to one, the local search is performed, and the best solution is updated.
The restricted candidate list RCL for our problem includes all demand points that can improve the objective while maintaining the waiting time constraint.If this list is empty, then key is set equal to zero and no local search is performed; otherwise, a candidate point is randomly selected and key is set to one, the local search is performed, and the best solution is updated.
Algorithm 4 shows the procedure for constructing the restricted candidate list.In this procedure, the so-called candidate parameter α ranges from zero to one.α 0 indicates that the  In this algorithm S1 is actually the left hand side of 2.10 .If demand point i i ∈ I \ I * j is added to I * j , then the left hand side of this equation is increased by In the local search, we use the 2-exchange neighborhood structure.The exchange between two demand points i ∈ I * j and s ∈ I\I * j is only possible if a s > a i , and, in addition, adding s to I * j and deleting i from I * j would not violate the constraints.The procedure for local search is shown in Algorithm 5.
If s ∈ I\I * j is added to and i ∈ I * j is deleted from I * j , then S1 is increased by R ssj t∈I * j R tsj R stj and decreased by R iij t∈I * j , t / i R itj R tij .In the above procedure, i ∈ I * j is replaced by a demand point r r ∈ I\I * j which has maximum population and would not violate the constraints.
Finally, in construction and local search phases, adding i ∈ I to I * j is possible only if the stability condition is not violated.

GRASP for Problem P
Our procedure for solving P is similar to that proposed for the maximal covering problem by Resende 31 .The main difference is in selection of the greedy function for construction of the restricted candidate list.The construction of the restricted candidate list is illustrated in Algorithm 6.
We use the GRASP procedure to solve the subproblem P j for any candidate facility j j ∈ J \ J * , as indicated in line 5.
There are few differences between our implementation of GRASP and that of Silva and Serra 18 for the M/M/1 case with priorities.In Silva and Serra 18 , the allocation is based on a set D ij which is constructed for each candidate center j.D ij contains the indices of all demand points i for which d i, j is less than R. Entries in D ij are ordered according to their distances from j. Demand points from D ij are assigned to j as long as the waiting time constraint is not violated.Another difference is in the selection of the greedy function.
In Silva and Serra 18 , the greedy function is the total costumer arrival rate, while ours is the amount of increase in the objective function.Finally their local search is a comprehensive search over all feasible solutions which improves the objective, while we use a 2-exchange neighborhood search.

Computational Results
The queueing maximal covering location-allocation problem with an M/G/1 is considerably more difficult to solve than the one with an M/M/1 system, since the waiting time quality of service constraints are quadratic with possibly indefinite coefficient matrices.Hence, even small-sized problems cannot be solved with the state-of-the-art software packages of today.Here we have tried to solve some randomly generated test problems by the proposed algorithms.We first considered the solution of problems with heuristic algorithm where subproblems P j were solved with GRASP.Then we used the GRASP to solve the entire model P .All calculations were performed on a Pentium IV processor with 2.80 GHz and 2.50 GB of RAM.
The number of points in the instances ranged from 20 to 200.In these instances, the number of candidate locations, n, was taken to be equal to the number of points, m.We assume that each demand point is also a potential server location, and the distances are considered to be Euclidean.
f i and τ j , namely, the daily call rate and the average time limit, were taken to be 0.005 times the population and 12.75 minutes for each candidate center.The covering radius, R, was taken to be 1.5 miles.Parameters α and β were set to 0.85 and 2; respectively, and v was set to 5 miles/hour .Z ı is assumed to be the same for all demand points and was taken to be 0.5 hour .The maximum number of iterations was set to 20.Each problem was run 5 times.Percent of coverage listed are the averages.Tables 1 and 2 show the results for instances solved with the suggested heuristic and with the GRASP procedure, respectively.
Comparison of the results for the two approaches reveals that solving the subproblems with GRASP takes far less time than does solving the entire model.The CPU time increases dramatically for the second algorithm as the number of points increases, while the coverage does not improve substantially.For example, as indicated in Tables 1 and 2, for the 200 point test problem, the first approach provides a solution with 100% coverage by selecting p 22 in 468.40 seconds, while the second approach gives a solution with 100% coverage with p 21 in 9905.76 seconds.
Finally, in order to have some idea about our solution methods, we tried to solve some small instances by enumerating all possible outcomes.The results for 5 and 10 points with 1 and 2 centers are shown in Table 3.For these limited instances, there were no discrepancies between the solutions.However, this of course is no proof that the algorithm always obtains an optimal solution, but we hope and expect it to at least obtain a rather good solution.

Conclusion
In this paper we considered the queueing maximal covering problem with an M/G/1 queueing system.A quadratically constrained integer programming model was presented for the problem.To solve the problem, we first considered a heuristic algorithm based on decomposing the problem into smaller knapsack-type subproblems.Then we solved the binary quadratic subproblems with GRASP.We also used GRASP to solve the entire model.Our limited computational results showed that using the proposed heuristic algorithm provide better solutions.
queueing systems and proposed heuristic methods to solve them.Corrêa and Lorena 21 proposed a constructive genetic algorithm, and Corrêa et al. 22 proposed a clustering search procedure to solve this problem.Corrêa et al. 23 convert the problem into a covering graph problem time at each center j does not exceed a predetermined amount τ j .The objective maximizes the population covered.Marianov and Serra 16 considered this problem for the M/M/1 and M/M/m Heuristic algorithm for solving P .

Algorithm 3 :
GRASP procedure for solving BQSP P j .Preprocess and CGRS-LS procedures.points are randomly selected, while α 1 yields the greedy selection.S1 and R itj are defined as

Table 1 :
Computational results for test problems solved with the heuristic procedure.

Table 2 :
Computational results for test problems solved with GRASP.

Table 3 :
Comparison of the results for enumerative method with the heuristic procedure and GRASP.