An Efficient Algorithm for Capacitated Multifacility Location Problems

: In this paper, a squared-Euclidean distance multifacility location problem with inseparable demands under balanced transportation constraints is analyzed. Using calculus to project the problem onto the space of allocation variables, the problem becomes minimizing concave quadratic integer programming problem. The algorithm based on extreme point ranking method combining with logical techniques is developed. The numerical experiments are randomly generated to test efficiency of the proposed algorithm compared with a linearization algorithm. The results show that the proposed algorithm provides a better solution on average with less processing time for all various sizes of problems.


INTRODUCTION
Facility location problem, searching for the optimum locations of facilities to attain minimum cost or maximum profit, is a significant problem that almost every organization has to face with. The appropriate locations for the facilities, such as warehouses, factories, service centers, outlets, clinics, and so on, to provide the minimum total distances between facilities and customers are completely required to enhance the organization's performance. Conversely, if the facilities are not located in appropriate locations, the company will not only be suffered from large investment but it will be suffered from low manufacturing and service performance for a long period as well. Thus, it is undeniable that the excellent decision on facility location is indispensable for all types of organization to enhance the core competence of the company. As a consequence, this decision problem has been studied continuously for many decades in both wider and deeper area of research.
Capacitated version of multifacility locationallocation problem (CMLP), which requires locating a set of facilities and simultaneously allocating to these facilities demands for service from a set of customers in order to optimize some performance criteria, is a specific class of facility location problem proven to be NP-hard. This means that it is hard or impossible to be solved by the exact methods when the problem is large.
Unsurprisingly, most of the developed algorithms for CMLP are heuristic algorithms, which although cannot provide the best solution, they give good solutions with much less computational effort than exact algorithms. CMLP has been studied in wild diverse versions, which can be classified as follows. The classical CMLP studied by Nauss [1] , Sa [2] , Akino and Khumawala [3] considered the problem on network and the products are separable or can be fraction. The p-median problem studied by Mulvey and Beck [4] , Koskosidis and Powell [5] , Lorena and Seene [6] considered the CMLP on network with inseparable products. The last class which was studied by Sherali and Tuncbilek [7] is to consider CMLP on plane with separable products. Therefore, the heuristic algorithms for CMLP have wild varieties depending on defining the problem.
In this paper, an efficient algorithm for another specific class of CMLP, which considers the problem on plane with inseparable product, is developed under balanced transportation constraints. The objective is to minimize total distance measured by squared-Euclidean metric between the facilities and customers. It will be hereafter called Capacitated Multifacility Location Clustering Problem (CMLCP). The applications of CMLCP can be usually found in computer network or electronics system setup such as finding the appropriate locations of computer servers and allocation of the clients to these servers in order to minimize losses due to distance between servers and clients.

MATERIALS AND METHOD
This part describes the proposed algorithm to solve CMLCP which can be divided into 3 sections as follows. Define the problem and formulate mathematical model: The studied problem is described as follows. There are 1 m > new facilities to be located on the continuous plane with a certain capacity. They have to serve n customers in their responsibilities whose locations and inseparable products or demands are known and deterministic. The objective is to find the good locations of these new facilities and allocation of customers to them so as to minimize the total distance measured in squared-Euclidean metric with respect to facility capacity. The problem can be mathematically formulated as follows. The objective function above gives the total transportation cost, while the first constraint set ensures that all customer demands are satisfied and the second constraint set ensures that all facility capacity limitation are respect. If the allocation variable ij z is fixed, the unconstrained minimum of the strictly convex objective function is readily obtained by partial derivatives of equation (1) to i x and i y at the solution Substituting (2) into the objective function of (1), the objective function can be projected onto the space of ij z variables as follows.
The equation (3), which is the reduced objective function (1), is equivalent to the following equation The problem (1) with objective function (4) can be written in the matrix form as follows.   Algorithms development: Observing that problem (6) is minimizing concave function or maximizing convex function subject to convex set of constraints, the optimal solutions occur at the extreme point of the convex set [8] . Therefore, the proposed algorithm is based upon extreme point ranking method, which is a specific algorithm for this kind of problem. Moreover, since the problem is 0-1 quadratic integer programming, then logic based method will be combined to reduce the number of variables by fixing value of some variables. After that, the next adjacent vertices will be explored by exchanging cluster of customers corresponding to the balanced transportation constraints and then be ranked using extreme point ranking approach. These methods and techniques added into the algorithm can be described as follows.
Extreme Point Ranking Method: The basic idea of extreme point ranking is ranking the vertices of the polytope defining the feasible region in order of importance regarding the global solution. Starting from one of the vertices of the polytope, the nearby vertices are ranked using an extreme point approach. This provides a new vertex to move to and the process continues until no adjacent vertices can be found with a decreasing objective function value. The initial vertex is found using the same method as Gupta et al.'s method [9] , while the next adjacent vertices are found using the proposed techniques. At each step, linear integer programming problem P2 shown below is solved to provide lower bounds on the objective function values of the quadratic integer programming problem P1 while the upper bound can be easily updated by substituting this solution in ( ) f z .
The above equation shows that, for all feasible solutions in S, P2 will not give the objective function value over that of P1. Thus, P2 is the lower bound of P1. Proposition 2, proven below, explains how ranking feasible solutions for P2 can provide the optimal solution to P1. Proposition 2: Let r z be r th extreme points ranked in ascending order of the objective function value of ( ) g z and k T be the set of all extreme points collected from the 1 st to the k th set of adjacent extreme points.
This means that ( *) f z is the least among the values of ( ) f z at all of the integer feasible solutions in S. Thus, * z is an optimal solution of problem P1.
To find c U and ( ) g z , Hessian matrix generally found by calculus, is needed. In this paper, the algebraic method is used to construct the closed form solution to each component of the Hessian matrix. The algebraic method arises from finding the closed from solution of the partial derivative terms corresponding to each , 1, ..., ; 1, ..., The closed form solutions Then, plugging these closed form solutions into their related position in the original Hessian matrix shown below. 1 Logic Based Method: The objective of using logic here is to tighten bound of some variables. Regarding to 0-1 quadratic problem, tightening bound of variables is equivalent to fix value of variables. This means that the number of variables to branch will be reduced. Therefore, the computational time is expected to be decreased. The logic used here is fixing 0 Observe that this logic will work only when some facilities have capacity limitation less than maximum amount of products shipped: Exchanging Method: Although the existing method for finding the next adjacent extreme points in quadratic integer programming problem is cutting plane method [9] , it was shown to be non-convergent due to cycling or need an infinite sequence of cutting planes by Zwart [10] . The reason for the cycling behavior as well as non-convergence of these approaches lies in the fact that although the approaches generate cones during the algorithm, they failed to explicitly incorporate these cones into the remaining step. This is essential to avoid the reemergence of vertices that have already been considered. To avoid the cycling, the exchanging method is proposed here. The problem here can be classified into clustering problem with balanced transportation constraints. Changing value of allocation variables between 0 and l, which means changing cluster of customers, affects both demand and supply constrains. To conserve the balance of the constraints a customer can move from a current facility to other facility only when there is another customer requiring the same amount of products or a group of some customers whose summation of amount of products equal to that of leaving customers to exchange with. The exchanging method proposed here exchanges customers only one pair of facilities at a time not consider crossing of the pair to avoid exponentially growth of computational time. The customers to be exchanged will be considered in order of appearing in vector of variable. Therefore, no cycling emerges. The exchanging method can be summarized as follows. Let t be number of customers considered to be moved out at a time. At the k th adjacent extreme point, t customers running from 1 to k customers of a current facility will be exchanged with k customers of the other facility whose summation of amount of products equal to that of the t customers. Observing that the maximum value of k is the maximum number of customers assigned to each facility.
The proposed heuristics, called EPR algorithm, can be summarized step by step as follows. Otherwise, go to step 2.
Step 2: Search for the new "current best solution", which will be the best incumbent solution to be searched for its next adjacent vertices, According to the exchanging method, not all possible (but some high possibility to be an optimal solution) extreme points are considered. Therefore, this optimal condition may not be satisfied and then 2 following additional stopping rules are constructed.
• No more possible exchanging pairs of customers exist. • There is no improvement of u f within q=2 consecutive sets of adjacent extreme points. Note that q can be more than 2. But, the higher value of q, the more computational effort is required. If at least one of these 2 additional stopping rules is satisfied, the existing current best solution is the final solution. And, * Verification and Validation of the Solutions: To verify and validate the proposed algorithm the numerical experiments are done. With MATLAB, both the proposed algorithm and exact algorithm will be applied to these sets of data. The solutions and computing time will be compared. The selected exact algorithm is linearization algorithm proposed by Glover and Woolsey [11] . To enhance the efficiency of exact algorithm, logic based method will be combined. With linearized objective function of (3), the problem becomes ( ) ( ) Observe from problem (9) that the last constraint is the least effective constraint. The experiments for solving the problems with and without the last constraints were done with various numbers of facilities and customers and the results show that the solutions for both of them are equivalent. Therefore, the last constraints will be ignored to reduce size of constraints when the problem is solved with "bintprog". This constraint is constructed to force ijk z to be one when both ij z and ik z are one. The value of ijk z should be set at one if there is no restriction on it to gain better objective function value. Therefore, if both ij z and ik z are one that allows ijk z to be one, the value of ijk z will be automatically one. As a result, the last constraint will be cut off and the constraints control value of ijk z will remain only There are ( 1) / 2 mn n− additional variables compared with problem (6) solved by the proposed algorithm. Problem (9) with ( 1) / 2 mn n + variables will be solved under branch and bound approach using command "bintprog" of MATLAB. The solutions obtained from this algorithm will be sequentially compared with ones obtained from the proposed algorithm.

RESULTS AND DISCUSSION
This part shows the results of the numerical experiments and discussion on these results. It is organized as follows. The first section analyzes the effect of Hessian construction time on processing time in order to ensure that most of processing time is devoted to solving the quadratic integer programming problem not to developing Hessian matrix. For second section, the results of numerical experiments with various problem sizes will be shown and discussed. . And, the number of variables, which equal the multiplication of m and n, runs from 10 to 240 variables. The measure of processing time will be divided into two parts: Hessian construction time and execution time. These times corresponding to problem sizes are shown in Table 1.

Number of variables = mxn
Observe from Table 1 that time used to construct Hessian matrix is still less than execution time if number of variables are not over 150 variables. For problems whose number of variables is not less than 150 variables, most of process time are devoted to constructing Hessian matrix. The Hessian matrix is theoretically developed by using calculus method, which is doing double partial derivative corresponding to each variable. Unsurprisingly, it takes much longer time developing Hessian matrix when even small number of variables increase. Time used for developing Hessian matrix can be reduced by using algebraic method instead of calculus method. To verify and measure the efficiency of constructing Hessian matrix using algebraic method, the same set of experiments will be solved by the same algorithm but new Hessian construction method. Process time using new Hessian construction method compared with the old one (from Table 1) can be shown in Fig 1. According to Fig.1, even the problem size grows, time to construct Hessian matrix by using algebraic method does not appear. Thus, process time is used only for solving problem and it seems to be equal to execution time of using calculus method to construct Hessian matrix. On contrary, time to construct Hessian matrix by using calculus method increases rapidly when problem size grows. Using algebraic method to construct Hessian matrix is very efficient. As a result, it will be used to develop Hessian matrix from now on. The results of evaluating EPR algorithm: The numerical experiments are constructed with various problem sizes, which number of variables corresponding to EPR algorithm = m n × no more than 1000 variables. The problem sizes vary from ( , ) m n = (2, 500) to (20, 50). For each problem size, 10-100 sets of data are generated and then solved by two algorithms: EPR algorithm and linearization algorithm. Both algorithms are coded with MATLAB and use command "bintprog" at solving step. The processing time and the following % of error for every case of EPR algorithms are collected, but only average of them will be reported.
OFV of EPR algorithm-OFVof exact algorithm % of error 100 OFVof exact algorithm = × OFV abbreviates from objective function value. Since the number of variable of linearization algorithm grows nonlinearly when n increase, there are some cases taking too much time. Therefore, the cases whose number of variable corresponding to the linearization algorithm = ( 1) / 2 m n n × × + over 400 will be premature terminated by time limitation. The level of time limitation is determined as follows.   Table 2. In Table 2, the blank cells represent the unsolvable cases due to number of variables over limitation of command "bintprog", while underlined values shows that there are some cases in the problem size are premature terminated. There are 10 sets of data for these both cases while 100 sets of data were generated for the other problem sizes. Processing time and % of error of EPR algorithm are separately measured into three parts: processing time to obtain initial solution ( 0 z ), incumbent solution ( 0 * z ), and final solution ( ; k r z z T ∈ ), so that the improvement rate of solution can be observed.
According to Table 2, EPR algorithm not only provides better solution, but also utilizes much less computational time than linearization algorithm for all problem sizes. Observe that, for all sizes m of facility processing time of linearization algorithm grows exponentially when n increases even there are premature terminated cases, while that of EPR algorithm almost disappear when m is small and grows linearly when m is larger.