Combined Simulated Annealing Algorithm for the Discrete Facility Location Problem

Abstract

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.

1. 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, LP-rounding [3–8], primal-dual [9–11], dual-fitting [12–14], and local search [15–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–21] and genetic algorithm [20, 22].

In this paper, a general algorithmic framework of combined simulated annealing (CSA) is developed for the DFLP, 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.

2. 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, 𝐼={1,2,…,𝑛} is used to represent the set of the customers and 𝐽={1,2,…,𝑚} is used to represent the set of the potential location sites. 𝑓_𝑗 is used to represent the fixed cost of opening a facility at site 𝑗∈𝐽, 𝑣_𝑗 is the service capacity of facility at site 𝑗, 𝑑_𝑖 is the demand of customer 𝑖∈𝐼, and 𝑐_𝑖𝑗 is the cost of serving one unit of demand at customer 𝑖 from site 𝑗, in other word, the unit variable shipping cost between customer 𝑖 from site 𝑗. We reasonably assume that 𝑐_𝑖𝑗≥0, 𝑓_𝑗≥0, and 𝑣_𝑗≥0 for all 𝑖∈𝐼, for all 𝑗∈𝐽.

And two binary variables are set as: 𝑋_𝑗=𝑌1ifafacilitysetuponsite𝑖;0otherwise;_𝑖𝑗=1iffacility𝑗servescustomer𝑖;0otherwise. (1)

Then the general model of DFLP could be stated as the following linear mixed-integer program: minΦ=^𝑚_𝑗=1𝑓_𝑗𝑋_𝑗+^𝑛𝑚_𝑖=1𝑗=1𝑑_𝑖𝑐_𝑖𝑗𝑌_𝑖𝑗Subjectto^𝑚_𝑗=1𝑌_𝑖𝑗=1∀𝑖∈𝐼_𝑛𝑖=1𝑑_𝑖𝑌_𝑖𝑗≤𝑣_𝑗𝑋_𝑗𝑌∀𝑗∈𝐽_𝑖𝑗−𝑋_𝑗𝑋≤0∀𝑖∈𝐼,∀𝑗∈𝐽_𝑗∈{}𝑌0,1∀𝑗∈𝐽_𝑖𝑗∈{0,1}∀𝑖∈𝐼,∀𝑗∈𝐽. (2)(3)(4)(5)(6)(7)

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.

3. Simulated Annealing Algorithm for DFLP

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.

In last 30 years, SA has been studied widely and used extensively in many optimization problems [24–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, 𝐒, 𝐒^′, 𝐒^′′, and 𝐒 are used to represent the different feasible solutions of the model; 𝐷(𝑇_𝑖) is the cooling function of temperature, in which 𝑇_𝑖 is the current temperature value. 𝑇_𝑓 is the stop temperature value. 𝑁 is the maximum iteration number at each temperature value. Φ(𝐒) is used to represent the objective function value of the solution 𝐒. According to Jayaraman and Ross [19], the SA for DFLP could be given as follows.

Step 1 (initialization). Set iteration counter 𝑖=1. Generate an initial feasible solution 𝐒 and regard 𝐒 as the optimal solution. Set the initial temperature 𝑇_𝑖 and the stop temperature 𝑇_𝑓 are specified. Define the cooling function 𝐷(𝑇_𝑖).

Step 2 (generate a feasible neighboring solution). Perform the neighboring function on current solution 𝐒 and get the new neighboring solution 𝐒^′.

Step 3 (evaluate current solution with neighboring solution). If the objective function value of the new solution 𝐒^′ is no less than that of the current solution 𝐒, namely, Φ(𝐒^′)≥Φ(𝐒), then proceed to Step 4; otherwise, if Φ(𝐒^′)<Φ(𝐒), then 𝐒=𝐒^′, proceed to Step 5.

Step 4 (examine metropolis condition). Determine the difference Δ𝐶 between the incumbent solution 𝐒 and the neighboring solution 𝐒^′, as ΔΦ=Φ(𝐒^′)−Φ(𝐒). Generate a random number 𝜌 from the interval (0,1), if 𝜌<exp(−ΔΦ/𝑇_𝑖), then 𝐒=𝐒^′. Proceed to Step 5.

Step 5 (check increment counters). Set 𝑖←𝑖+1. If 𝑖≤𝑁, then return to Step 2. Otherwise proceed to Step 6.

Step 6 (adjust temperature). Adjust temperature by the cooling function, Mathematically this is 𝑇_𝑖←𝐷(𝑇_𝑖).

Step 7 (convergence check). If 𝑇_𝑖≥𝑇_𝑓, then reset 𝑖=1 and return to Step 2. Otherwise, stop and output the optimal solution 𝐒.

4. Combined Simulated Annealing (CSA)

The solution of DFLP includes two parts: 𝑋_𝑗s and 𝑌_𝑖𝑗s. 𝑋_𝑗 denotes whether or not open the facility at site 𝑗, while 𝑌_𝑖𝑗 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 𝑌_𝑖𝑗 is subject to the variable 𝑋_𝑗. 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 𝑋_𝑗. The external layer subalgorithm (ELSA) would perform the allocation optimization under fixed 𝑋_𝑗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 𝐒, 𝐒^′, 𝐒^′′, 𝐒 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].

4.1. 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 𝑀(∑^𝑚_𝑗=1𝑋_𝑗<𝑀), then select a candidate site 𝑖 which satisfies 𝑋_𝑗=0 in current solution 𝑆 randomly and set 𝑋_𝑗=1, namely, there would locate a new facility.(2)If the number of the locate facility no less than 1(∑^𝑚_𝑗=1𝑋_𝑗≥1), then select a site 𝑗 randomly which satisfies 𝑋_𝑗=1 in solution 𝑆 and set 𝑋_𝑗=0, namely, there would close a open facility.(3)Select a site 𝑗 which satisfies 𝑋_𝑗=0 and another site 𝑗^′ which satisfies 𝑋_𝑗′=1 in current solution 𝑆, then set 𝑋_𝑗=1 and 𝑋_𝑗′=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 𝑌_𝑖𝑗 and 𝑌_𝑖𝑗′ that satisfy 𝑌_𝑖𝑗=1 and 𝑌_𝑖𝑗′=0; then set 𝑌_𝑖𝑗′=1, 𝑌_𝑖𝑗=0, namely, it allocates the demand of customer 𝑖 from facility 𝑗 to facility 𝑗^′.(2)Exchange the facilities which serve two customers re spectively with each other. To be specific, select four allocation variables as 𝑦_𝑖1𝑗1=1, 𝑦_𝑖2𝑗2=1, 𝑦_𝑖1𝑗2=0, 𝑦_𝑖2𝑗1=0, then set 𝑦_𝑖1𝑗2=1, 𝑦_𝑖2𝑗1=1, 𝑦_𝑖1𝑗1=0, 𝑦_𝑖2𝑗2=0.

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.

4.2. The Procedure of CSA

The following is a step-by-step description of the procedure of CSA.

Step 1 (initialization). Set iteration counter 𝑖=1, 𝑘=1. Set the initial temperature 𝑇_𝑖, and stop temperature 𝑇_𝑓, the initial feasible solution is S, and let S be the optimal solution at the same time. Define the cooling function 𝐷(𝑇_𝑖). 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 𝐒^′.

Step 4 (perform ILSA). Step 4.1. Set 𝑘=1. Regard solution 𝐒^′ as the initial solution and the current optimal solution.
Step 4.2. Perform the internal layer neighboring function on solution 𝐒^′, and generate a neighboring feasible solution 𝐒^′′.
Step 4.3. If Φ(𝐒^′′)<Φ(𝐒^′), then set 𝐒^′=𝐒^′′; otherwise, generate a random number 𝜌 from (0,1), if 𝜌<exp(−(Φ(𝐒^′′)<Φ(𝐒^′))/𝑇_𝑖), then set 𝐒^′=𝐒^′′.
Step 4.4. Consider 𝑘←𝑘+1. If 𝑘≤𝑁₂, then return to Step 2. Otherwise proceed to Step 5.
Step 4.5. Stop and output the optimal solution 𝐒^′, return to the ELSA.

Step 5 (save the global optimal solution). If Φ(𝐒^′)<Φ(𝐒), then 𝐒=𝐒^′.

Step 6 (evaluate current solution). Evaluate current solution with neighboring solution. If Φ(𝐒^′)<Φ(𝐒), then let 𝐒=𝐒^′ proceed to Step 7; otherwise, proceed to Step 7.

Step 7 (examine metropolis condition). A random number 𝜌 is generated from (0,1), if 𝜌<exp(−(Φ(𝐒^′)−Φ(𝐒))/𝑇_𝑖) then 𝐒=𝐒^′. Proceed to Step 8.

Step 8 (increment counters). Set 𝑖←𝑖+1. If 𝑖≤𝑁₁, then return to Step 3; otherwise, proceed to Step 9.

Step 9 (adjust temperature). Adjust temperature by the cooling function: 𝑇_𝑖←𝐷(𝑇_𝑖).

Step 10 (convergence check). If 𝑇_𝑖≥𝑇_𝑓, then reset 𝑖=1 and return to Step 3. Otherwise, stop and output the optimal solution 𝐒.

𝑁₁,𝑁₂ 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.

5. 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 𝑇₁ 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 𝑁₁=𝑚, 𝑁₂=5𝑛, the stop temperature value 𝑡_𝑓=0.001. The SA is used to solve the instances too, and its iteration maximum under the same temperature is 5𝑚𝑛; 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.

Table 1 

Comparison for OR-Library instances.

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.

Figure 1 

OFV versus iteration number in Cap101.

Figure 2 

OFV versus iteration number in Cap131.

Example 2 (The Ghosh Instances). The Ghosh instances are 90 instances in total, and with 𝑛=𝑚 on all cases. These instances of the same size are divided into two categories, symmetric, where the transportation cost 𝑐_𝑖𝑗=𝑐_𝑗𝑖 holds, and asymmetric, where 𝑐_𝑖𝑗=𝑐_𝑗𝑖 does not necessarily hold. Each group contains three values of 𝑚×𝑛: 250×250, 500×500, and 750×750. The instances in each category are further divided into three instance classes, and they differ in the range of values from which opening costs 𝑓_𝑗 are drawn, which could be chosen from the [100,200] (Range A), [1000,2000] (Range B), and [10000,20000] (Range C) randomly.
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.
As shown in 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 𝑚×𝑛=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 𝑚×𝑛=250×250 in type C.

Table 2 

Comparison for the Ghosh instances.

6. 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 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.

Acknowledgment

This research was supported by the National Natural Science Foundation of China (Grant no. 71101155).