Combined Simulated Annealing Algorithm for the Discrete Facility Location Problem

The combined simulated annealing (CSA) algorithm was developed for the discrete facility location problem (DFLP) in the paper. The method is a two-layer algorithm, in which the external subalgorithm optimizes the decision of the facility location decision while the internal subalgorithm optimizes the decision of the allocation of customer's demand under the determined location decision. The performance of the CSA is tested by 30 instances with different sizes. The computational results show that CSA works much better than the previous algorithm on DFLP and offers a new reasonable alternative solution method to it.


Introduction
The classical facility location problem (FLP) is one of the most important models in combinatorial optimization, which is to determine the number and locations of the facilities and allocate customers to these facilities in such a way that the total cost is minimized. The FLP may be the most critical and most difficult decision in the designing of an efficient supply chain for the facilities are costly and difficult to reverse after being located. The problem is also encountered in other areas such as material distribution, transportation network, and telecommunication network.
The FLP can be classified in two categories as discrete problem and continuous problem according to whether the sets of demand points and facility locations are finite. The discrete facility location problem (DFLP) assumes that the solution space is discrete and generally the facilities are located on the nodes of the network, which brings a lot of complexities to the problem. And many practical problems without facilities to locate, such as cluster analysis, machine scheduling, economic lot sizing, portfolio management, and computer network design, can also be modeled as DFLP [1]. Due to its strategic nature, DFLP has been widely studied by researchers over many years, especially developing solution methods for the DFLP has been a hot topic of research for the last 30 years. Many successful contributions of the DFLP have been reported both in theory and in practice. The Daskin [2] and Melo et al. [1] provided thorough reviews of the DFLP.
The DFLP is a NP-hard problem, and this nature makes the exact algorithms only for small problems and heuristics the natural choice for larger instances. Therefore much attention has been focused on designing heuristics algorithms with good performance. Following the first heuristics algorithm presented by Shmoys et al. [3], there was a long list of work on designing heuristics algorithms for this problem over the years. As a result, four basic algorithms with different features have been emerged, namely, LProunding [3][4][5][6][7][8], primal-dual [9][10][11], dual-fitting [12][13][14], and local search [15][16][17][18]. Although LP-rounding could be used to design algorithm with better results than other three, it is noncombinatorial in nature and needs more CPU time. Primal-dual and dual-fitting method could be adapted to solve variants of the FLP but are less robust. Recently some applications have proved that local search is more powerful on the hard DFLP [16]. And along with the development of the computation method, more and more researchers use the heuristic algorithms based on the local search to solve the problem, such as simulated annealing [19][20][21] and genetic algorithm [20,22].
In this paper, a general algorithmic framework of combined simulated annealing (CSA) is developed for the DFLP, The Scientific World Journal and its performance is compared with other existed solution methods.
The rest of this paper is organized as follows. In Section 2, we describe the general model of DFLP. The basic procedures of the SA and CSA are presented in Sections 3 and 4 respectively. Section 5 shows the computational results on test instances, and Section 6 concludes the paper.

The Formulation of DFLP Model
We start by giving the mathematical formulation of a general model for DFLP with a minimized objective function. In the model, I = {1, 2, . . . , n} is used to represent the set of the customers and J = {1, 2, . . . , m} is used to represent the set of the potential location sites. f j is used to represent the fixed cost of opening a facility at site j ∈ J, v j is the service capacity of facility at site j, d i is the demand of customer i ∈ I, and c i j is the cost of serving one unit of demand at customer i from site j, in other word, the unit variable shipping cost between customer i from site j. We reasonably assume that c i j ≥ 0, f j ≥ 0, and v j ≥ 0 for all i ∈ I, for all j ∈ J.
And two binary variables are set as: Then the general model of DFLP could be stated as the following linear mixed-integer program: The objective function (2) is to minimize the total system cost, including the location cost and the shipment cost. Constraint (3) is the demand constraint, which makes the demand of each customer be met; (4) is the variable upper bound constraint; (5) is the capacity constraint of facility; (6) and (7) are standard binary integrality constraints. Kirkpatrick et al. [23] introduced the concept of simulated annealing (SA) algorithmin 1983, which is a stochastic optimization technique. To be specific, SA is a probabilistic heuristic for the global optimization problems of finding a good approximation to the global optimum of a given objective function in the search space. It is often used when the solution space is discrete. In the searching process, the SA accepts not only better but also worse neighboring solutions with a certain probability. This means that the SA has the ability to escape from the local minima. Therefore, it can find high-quality solutions that do not strongly depend upon the choice of the initial solution compared to other local search algorithms. And its another advantage over the other heuristic algorithms is the ease of implementation. So we adopted SA as the basic solution method to solve the DFLP.

Simulated Annealing Algorithm for DFLP
In last 30 years, SA has been studied widely and used extensively in many optimization problems [24][25][26][27][28][29], which have proved that SA is an effective tool for approximating globally optimal solutions to many NP-hard optimization problems.
In order to describe the procedure of the SA, S, S , S , and S are used to represent the different feasible solutions of the model; D(T i ) is the cooling function of temperature, in which T i is the current temperature value. T f is the stop temperature value. N is the maximum iteration number at each temperature value. Φ(S) is used to represent the objective function value of the solution S. According to Jayaraman and Ross [19], the SA for DFLP could be given as follows.
Step 1 (initialization). Set iteration counter i = 1. Generate an initial feasible solution S and regard S as the optimal solution. Set the initial temperature T i and the stop temperature T f are specified. Define the cooling function D(T i ).
Step 2 (generate a feasible neighboring solution). Perform the neighboring function on current solution S and get the new neighboring solution S .
Step 3 (evaluate current solution with neighboring solution). If the objective function value of the new solution S is no less than that of the current solution S, namely, Φ(S ) ≥ Φ(S), then proceed to Step 4; otherwise, if Φ(S ) < Φ(S), then S = S , proceed to Step 5.
Step 5 (check increment counters). Set i ← i + 1. If i ≤ N, then return to Step 2. Otherwise proceed to Step 6.
Step 6 (adjust temperature). Adjust temperature by the cooling function, Mathematically this is T i ← D(T i ).
Step 7 (convergence check). If T i ≥ T f , then reset i = 1 and return to Step 2. Otherwise, stop and output the optimal solution S.

Combined Simulated Annealing (CSA)
The solution of DFLP includes two parts: X j s and Y i j s. X j denotes whether or not open the facility at site j, while Y i j denotes the service demand allocation. The two variables are interdependent and interactional. As each demand must allocate to an opening facility, we could can conclude that the variable Y i j is subject to the variable X j . This relation also can be seen from the constrain (4) in the model in Section 1.
CSA works in two layers as internal layer and external layer to solve the problem. The internal layer subalgorithm (ILSA) would optimize the facility location decision variable X j . The external layer subalgorithm (ELSA) would perform the allocation optimization under fixed X j s which determined in the internal layer. In each layer the SA is used and they make up of the CSA. The method that divides the problem into two layers could make the search in the procedure explores in smaller solution space each time, so it increases the probability of obtaining the global optimal solution.
Some parameters should be initialized before the performance of CSA, including the initial temperature and stop temperature, cooling ruler of the temperature, iteration maximum in each temperature. The initialization of the parameters could be determined in similar ways to these in the SA. And the S, S , S , S are also used to represent the different feasible solutions here.
In addition, CSA must start with an initial solution or with a solution produced using a heuristic. In this work, we use the randomly generated initial solution, which proposed in Qin et al. [30].

Neighboring Functions.
Similar to the SA, the CSA algorithm is an iterative search procedure based on the neighboring function. The quality of the optimal solution is very sensitive to the way that the candidate solutions are selected. Thus, the neighboring function is crucial to the good performance of the CSA algorithm.
The ELSA would optimize the facility location decision. So the neighboring function of the ELSA to modify the configuration of the current solution and generate a neighboring solution could perform three different operations: (1) If the number of the located facilities is less than the allowed maximum M ( m j=1 X j < M), then select a candidate site i which satisfies X j = 0 in current solution S randomly and set X j = 1, namely, there would locate a new facility.
(2) If the number of the locate facility no less than 1 ( m j=1 X j ≥ 1), then select a site j randomly which satisfies X j = 1 in solution S and set X j = 0, namely, there would close a open facility.
(3) Select a site j which satisfies X j = 0 and another site j which satisfies X j = 1 in current solution S, then set X j = 1 and X j = 0.
In the implementation of the ELSA, we could select only one operation from the above three operations to perform the neighboring function each time. And after implementing the operation, it should allocate the customers' demand to the opened facilities again, as generate an initial solution again.
ILSA is to determine the demand allocation decision. According to the features of the allocation decision, there are two operations to generate the neighboring solutions in the ILSA.
(1) Select two allocation variables Y i j and Y i j that satisfy Y i j = 1 and Y i j = 0; then set Y i j = 1, Y i j = 0, namely, it allocates the demand of customer i from facility j to facility j .
Similarly, the neighboring function of the ILSA could select only one operation to perform each time. In addition, it must ensure that demands of all customers must be satisfied and the facilities have no capacity violations exist. Otherwise, it should return to the old solution and reselect an operation to perform.

The Procedure of CSA.
The following is a step-by-step description of the procedure of CSA.
Step 1 (initialization). Set iteration counter i = 1, k = 1. Set the initial temperature T i , and stop temperature T f , the initial feasible solution is S, and let S be the optimal solution at the same time. Define the cooling function D(T i ). Generate a feasible solution by randomly allocation with capacity restricted.
Step 2 (check feasibilities). The method now checks the demand allocations to ensure that no capacity violations exist. The demands of the customers are also checked. If the solution S is not feasible, we should return to Step 1.
Step 3 (generate a neighboring solution). Perform the external layer neighboring function on solution S, and generate a neighboring feasible solution S .
Step 4.1. Set k = 1. Regard solution S as the initial solution and the current optimal solution.
Step 4.2. Perform the internal layer neighboring function on solution S , and generate a neighboring feasible solution S .
Step 4.4. Consider k ← k + 1. If k ≤ N 2 , then return to Step 2. Otherwise proceed to Step 5.
Step 4.5. Stop and output the optimal solution S , return to the ELSA.
Step 9 (adjust temperature). Adjust temperature by the cooling function: T i ← D(T i ).
Step 10 (convergence check). If T i ≥ T f , then reset i = 1 and return to Step 3. Otherwise, stop and output the optimal solution S. N 1 , N 2 are the given maximum iteration number in ELSA and ILSA respectively. The Step 5 is to save the global optimal solution that has been found so far in the CSA. This operation does not take the acceptance probability of worse solution into consideration. So it could help the algorithm avoid losing the global optimal solution.

Computational Experiments
To assess the practical effectiveness of the proposed CSA algorithm, we use 30 instances with different sizes as benchmark problems. Twelve "Capacitated warehouse location" instances were given by Beasley [31], which are publicly available from the OR-Library and could be downloaded directly from the website: http://people.brunel .ac.uk/∼mastjjb/jeb/orlib/capinfo.html. The other 18 instances were proposed by Ghosh [32].
To perform the CSA, the initial temperature T 1 could be set equal to the objective function value of the initial solution. The cooling ruler is equal-ratio cooling, and the cooling ratio α = 0.95. The iteration maximum N 1 = m, N 2 = 5n, the stop temperature value t f = 0.001. The SA is used to solve the instances too, and its iteration maximum under the same temperature is 5mn; other parameters are same as CSA. So the total iteration numbers in the SA and CSA are equal. We also use the randomly generated method to find an initial feasible solution.
The model and CSA were implemented in Visual C# 2010. A personal computer with Intel E5800 CPU, 2G RAM and Windows XP Profession operating system was used for all tests.
Example 1 (OR-Library Instances). The instances in OR-Library are more than 150, and we selected 12 instances with different size to be solved by the CSA and compared with theirs optimal solution (all instances have been exactly solved by the Lindo software and provide their optimal solutions by the author).
Computational results of these OR-Library instances are reported in Table 1. Each row of the table gives the results of one individual instance. The instance names are originally used in the OR-Library. The gap in the table represented the relative error between the result and the optimal solution.
It can be observed from Table 1 that the CSA found optimal solutions for all these instances without any exception, while the SA didn't find anyone. The CPU time used by SA for each instance is from 3 to 8 times of that used by CSA.
The computation process of instances Cap101 and Cap131 with CSA and SA are depicted in Figures 1 and 2. The OFV is the objective function value. As shown in the figures, we can find that the convergence speed of CSA is more quickly than that of SA. The iteration time for converging to the optimal solution used by SA for each instance is from 10 to 25 times of that used by CSA.
Example 2 (The Ghosh Instances). The Ghosh instances are 90 instances in total, and with n = m on all cases. These instances of the same size are divided into two categories,  The optimal solutions of the Ghosh instances were not presented, but there used to be several methods to solve the instances and compared the results with each other in the literature [33]. And we compare the results of the CSA with them as reported in Table 2, in which we reported the results of the instances obtained by hybrid, CLM, GTS, TS, and CSA. Table 2, for all instances, CSA found solutions better than those of other methods with only five exceptions, and found solutions with the same result for symmetric instance with m × n = 500 × 500 in type C. Even in the five exceptions, the gaps between the solutions were found by CSA and the best solutions found by other methods are very small, and the maximal relative gap is only 3.1% in asymmetric instance with m × n = 250 × 250 in type C.

Conclusions
The DFLP is to determinate the facility location and demand allocation so as to minimize the total cost, based on the demands of customers satisfied without violating the 6 The Scientific World Journal capacity restriction of any facility. The CSA was proposed for the DFLP, which is a two-layer algorithm. Its external layer subalgorithm optimizes the facility location decision, while the internal layer subalgorithm optimizes the demand allocation under the fixed facility location which is determined by the external layer. The each local search process in CSA focuses on the smaller solution space, which could not only increase the probability of obtaining the global optimal solution, but also save computational time.
The performance of CSA is evaluated with 30 instances with different sizes. Those instances are solved by CSA in a very reasonable amount of time, and the solutions are compared with that of previous studies in the literature. It is showed that the new algorithm could give better results (or at least same) than the others for nearly all instances. Hence, the CSA works much better than the previous work and offers a reasonable alternative solution method to the DFLP.