A Simulated Annealing Methodology to Multiproduct Capacitated Facility Location with Stochastic Demand

A stochastic multiproduct capacitated facility location problem involving a single supplier and multiple customers is investigated. Due to the stochastic demands, a reasonable amount of safety stock must be kept in the facilities to achieve suitable service levels, which results in increased inventory cost. Based on the assumption of normal distributed for all the stochastic demands, a nonlinear mixed-integer programming model is proposed, whose objective is to minimize the total cost, including transportation cost, inventory cost, operation cost, and setup cost. A combined simulated annealing (CSA) algorithm is presented to solve the model, in which the outer layer subalgorithm optimizes the facility location decision and the inner layer subalgorithm optimizes the demand allocation based on the determined facility location decision. The results obtained with this approach shown that the CSA is a robust and practical approach for solving a multiple product problem, which generates the suboptimal facility location decision and inventory policies. Meanwhile, we also found that the transportation cost and the demand deviation have the strongest influence on the optimal decision compared to the others.


Introduction
The facility location problem (FLP) is one of the most important models in the combinatorial optimization problem, which is to determine the number and locations of facilities and allocate customers' demands to these discrete located facilities in such a way that the total cost is minimized. Due to its wide application in numerous contexts such as regional planning, telecommunication and transportation infrastructure deployment, and energy management, a lot of problem variants and theirs models, which capture features of the real facility location problems, have been developed in the literature since the original formulation in Weber's work [1].
The FLP is NP-hard, which can be classified in two categories as capacitated problem and uncapacitated problem according to whether it includes the capacity of facilities. The capacitated facility location problem (CFLP) assumes that each facility can produce or ship limited quantities of the product and the demand of each customer is satisfied without violating the capacity restriction of any facility, which makes the CFLP more widely used, but also brings a lot of complexities to the problem.
Daskin [2] and Melo et al. [3] provided thorough reviews of the CFLP models. And we could found that many past models have several common characteristics, such as deterministic demands, a single product. Obviously, these models are insufficient to cope with many realistic capacitated facility location problems. So many extensions to the basic problems have been proposed and studied extensively.
Stochastic capacitated facility location problems have been proposed to approach situations in which some parameters are uncertain or stochastic. And in reality the demands are typically stochastic, which would bring risk to the inventory management. In order to deal with the stochastic demands and improve the service level, managers often hold some products as stock, which can result in increased inventory cost. In many contexts where product safekeeping is expensive, the inventory cost may account for a significant portion of the total system cost. So measuring the trade-off between the inventory cost and the service level is naturally very important in the facility location decision. However, inventory costs were not typically considered in many FLPs.
Shen et al. [16] proposed a joint location model and restructured it into a set-covering model, and solutions only for two special cases were discussed; one is that variance of demand is proportional to the mean and the other is that demand has zero variance. Shu et al. (2005) [17] developed a more efficiency algorithm for the special cases. Shen [18] simplified the inventory cost of the products as the nonlinear function of the product quantity to the objective function and proposed the Lagrange algorithm to solve the problem; Shen and Qi (2007) [19] analyzed the transportation cost combined with the vehicle routing in the logistics network, but the order number was considered as a continuous variable in the formulation derivation. Kutanoglu and Lohiya [20] regarded the inventory cost as the linear function of product quantity. Ozsen et al. [21,22] and Sourirajan et al. [23] proposed two different extensions to the model in Shen et al. [16]. Chen et al. [24] optimized facility locations, customer allocations, and inventory management decisions when facilities were subject to disruption risks. Hui et al. [25] studied the locationinventory results with joint replenishment (JR) and independent replenishment (IR), respectively, and found that the JR policy can obtain better solutions than IR policy.
In this paper, we study the stochastic multiproduct capacitated facility location problem (SMCFLP) that seeks the optimal solutions to minimize the total cost, including the facility setup cost, order cost, day-to-day shipment cost, and inventory cost, in which the demands for multiproducts are stochastic and satisfy the normal distribution. An open facility could serve one or more customers, while each customer could only be allocated to exactly one facility. The present work is motivated by a study of a local delivery service at a machinery manufacturing enterprise.
The remainder of the paper is organized as follows. Section 2 presents some rational assumptions and notations of the SMCFLP. In Section 3, based on system cost analysis, we propose the model formulation that incorporates inventory control decision into multiproduct facility location problem with stochastic demand. In Section 4, an improved simulated annealing algorithm is developed to solve the model, and the setting of parameters is also discussed. In Section 5, the numerical example and the analysis of the results are given.

Notations and Assumptions
In formulating the optimization model, the following notations are used. denotes the candidate facility sites, = 1, 2, . . . , .
is the storage capacity occupation per unit product .
is the fixed setup cost of facility .
is the storage capacity of facility .
, are the mean and the standard deviation of the demand for product of customer , respectively.
, are the mean and the standard deviation of the demand for product of facility , respectively. is the lead time (in days); that is to say, the supplier takes the time to fulfill an incoming order from facility for product .
is the fixed order quantity per order.
ℎ is the inventory cost per unit of product of facility (per day).
is the fixed cost of placing an order for product in facility .
is the transportation cost per unit to ship product from the supplier to facility . is the elapsed time between two consecutive orders for product in facility .
denotes the distribution cost per unit product between the facility and the customer .
is the united service level in the system, 0 < < 1; namely, the fill rate of all demands for all products must not be less than in all facilities.
is the bank interest rate.
is discount rate (calculated by ).
is the planning horizon (in year).
is the maximum number of the facilities that are allowed to locate.
Note that if = 0, no products will pass facility and vice versa. According to the definition of the notations, we could know that = . (2) And the fixed setup costs of the facilities , for all , are all disposable, and the other costs are always invested per day, so in order to keep the consistency of all costs in time, the fixed setup costs should be shared in day in the planning horizon. And according to the above definitions, the discount rate could be calculated as follows: Some rational assumptions are proposed as follows.
(1) Similar to Shen et al. [16], we assume that all stochastic customer demands , for all , , are assumed to follow standard normal distribution and are independent of each other.
(2) Each customer demand point is directly served by one and only one facility; namely, the demand of a customer could not be partitioned.
(3) According to Shen et al. [16], we assume that all facilities have the same service level; in the other words, the fill rates for the demands in all facilities are identical in the system.
(4) According to Ozsen et al. [22] and Chen et al. [24], we assume facility , for all using an economic order quantity (EOQ) model ( , ) for inventory policy; that is to say, a fixed quantity is ordered to the supplier, once the inventory quantity falls to or is below a reorder-point .
According to assumption (1), for the rest of the paper, we could get the following evident relations: From assumption (1), we also could know that the demands of all facilities are also stochastic and follow the standard normal distribution.

Model Formulation
According to assumption (4), the facility performs ( , ) inventory policy to meet the stochastic demand pattern. But even if the order is triggered, the ordered products could be received after days. So once an order is submitted, the inventory products should cover the demand produced during lead time with a given probability . This probability is known as the level of service for the system or the fill rate for the demand. So the service level constraints in the facility can be expressed as follows: where ( ) is the stochastic demand for product in facility produced with the lead time .
According to Shen et al. [16,18], we could know the average inventory cost rate of product in facility as where is the value of standard normal variate, which accumulates a probability of .
The operation cost for product in facility during the period is We could get the operation cost rate of product in facility , based on expression (9) divided by ( = / ): So the total operation cost rate is The transportation cost rate of the entire system is The total fixed setup cost of the system is As mentioned before, the fixed setup cost is one-time expenditure in the planning horizon, but the other costs are counted by day, so the setup cost should be converted into day-cost in order to be consistent with all costs in unit: Therefore, the total cost rate of the system is 4 The Scientific World Journal The objective of the problem is to minimize the cost (15). But there are too many uncertain variables in it, which could make the solution approach difficult to be designed.
Assume that the cost function (15) is continuously differentiable on order quantity, then performing differentiation on it in terms of , for all , , and letting the derivation be equal to zero (minimizing the total cost in a centralized approach), we could obtain Then we could get the optimal order quantity of product per order in facility as follows: Replacing (4), (5), and (16) in expression (14), the cost function (15) can be expressed as follows: So, we could present the model to solve the SMCFLP considering the inventory cost as follows: ≤ , ∀ , , , The model would determine where to open facilities no more than in candidate sites and how to serve customers with stochastic demand under service level constraint. The objective function (18) is to minimize the total system cost, including location cost, inventory cost, order cost, and transportation cost. Constraint (19) states that each demand must be only allocated to one facility. Constraint (20) ensures that the total demands allocated to a certain facility should not exceed its capacity. Constraint (21) restricts a demand only to be allocated to an open facility. Constraint (22) is the maximum restriction of the opened facilities. Constraints (23) and (24) are the standard integrality constraints, and ∈ {0, 1} representing single-sourcing constraints, which means that all of the demands of a customer must be assigned to the same facility.

Solution Approach
The above SMCFLP model (19)-(24) is a large-scale nonlinear mixed-integer programming model, which is NP-hard problem [16]. We were able to easily compute optimal solutions for small problem instances. However, as the number of products, number of facilities, and number of customers approach some practical size, the attractiveness of this methodology to provide optimal solutions to such practical sized problems in a reasonable amount of time deteriorates considerably. Its NP-hard nature makes the exact algorithms only for small problems and heuristics the natural choice for larger instances. There was a long list of works on designing heuristics algorithms for this problem over the years.
The Lagrangian method (Miranda and Garrido, 2008 [26]; Nezhad et al, 2013 [27]) and branch and bound algorithm [28,29] are always used to solve the similar problem. But the Lagrangian method and branch and bound method are intuitive, inflexible, and complicated. In order to increase the adaptability of the solution method, we therefore began testing heuristic approaches to this combinatorial problem. Several approaches were tried including genetic algorithm (GA) and simulated annealing (SA) algorithm.
Kirkpatrick et al. [30] introduced the SA; then, it has proved to be an effective tool for approximating globally The Scientific World Journal 5 optimal solutions to many NP-hard optimization problems. We adopted the SA procedure because of its ability to quickly combine the facility location decision and inventory control decision into a single large problem. In the last 30 years, SA has been applied to many optimization problems in a wide variety of areas such as facility layout [31][32][33], job scheduling [34][35][36], and network design ( [37][38][39][40][41][42][43]).
The solution of SMCFLP includes two parts, and , where denotes whether or not to open the facility in site and denotes the service allocation of customers' demand. The two variables are interdependent and interactional. As each demand must be allocated to an opening facility, we could think that is determined by . This relation also can be seen from constraint (22) in the model.
According to the above characteristics of the problem, we could use the combined simulated annealing (CSA; [43]) algorithm to solve the model. The CSA works in two layers as the outer layer algorithm (OLA) and the inner layer algorithm (ILA), in which the OLA optimizes the facility location decision. Then, the ILA optimizes the demand allocation decision based on the determined facility location decision determined in the OLA. In each layer the SA is used and the combination of them is CSA.
It is well known that the neighboring function is crucial to the good performance of the SA. In the algorithm we should design different neighboring functions in OLA and ILA, respectively.
The OLA would optimize the facility location decision. So the neighboring function of the OLA to modify the configuration of the current solution and generate a neighboring solution could have three different operations.
(1) If the number of the opened facility is less than the allowed maximum , then select a candidate site which satisfies = 0 in current solution randomly and set = 1; namely, a facility is located in the site .
(2) If the number of the opened facility is no less than 1, then select a site which satisfies = 1 in current solution randomly and set = 0; namely, the opened facility as site is closed.
(3) Exchange two facilities with different status; namely, select a site which satisfies = 0 and another site which satisfies = 1 in current solution , and then set = 1 and = 0.
In the implementation of the OLA, we could only select one operation from the above three operations to perform the neighboring function each time. And after implementing this OLA operation, it should allocate the customers' demand to the opened facilities again.
ILA 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 ILA.
(1) Select two allocation variables and that satisfy = 1 and = 0; then set = 1, = 0; namely, it allocates the demand of customer for product from facility to another facility .
(2) Exchange the facilities which service two customers, respectively, with each other. To be specific, select four allocation variables as 1 1 = 1, 2 2 = 1, 1 2 = 0, and 2 1 = 0; then set 1 2 = 1, 2 1 = 1, 1 1 = 0, and 2 2 = 0. Similarly, the neighboring function of the ILA could select only one operation to perform each time. In addition, it should ensure that demands of all customers must be satisfied and the facilities have no capacity violations that exist in the neighboring solution. Otherwise, it should return to the old solution and reselect an operation to perform.
In addition, it is well known that the SA algorithm 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 is similar to that in Qin et al. [43].
Let Φ( ) denote the objective function value (total cost) of solution ; then, the steps of the CSA could be given as follows.
Step 0 (initialization). Set the initial temperature , the cooling rate , the stop temperature , and the maximum iteration number at each temperature value. Set counter number Step 1 (generate initial solution). The initial solution could be generated by randomly allocating customers' demands to the facilities. Let the global optimal solution = .
Step 2 (check feasibility). The method now checks the demand allocations to ensure no capacity violations exist. The demands of the customers are also checked. If the solution is not feasible, we should return to Step 1.
Step 3 (generate a feasible neighboring solution). Perform the outer layer neighboring function on , and get the new feasible solution .
Step 4 (perform ILA) Step 4.1. Regard solution as the initial solution and the current optimal solution in the ILA; then, generate the ILA neighboring solution based on solution .
Step 4.4. If 2 ≤ , then go to Step 4.1; otherwise, set 2 = 1, stop the ILA computation, and return to the OLA.
Step 8 (adjust temperature). Adjust temperature by the cooling rate . Mathematically this is ← ⋅ . And set 1 = 1; then return to Step 3.
Step 9 (convergence check). If the temperature is greater than or equal to the given stopping value , then go to Step 3; otherwise, stop the computation and output the optimal solution . Note that Step 5 is to save the global optimal solution that has been searched so far in the algorithm. This operation does not take the acceptance probability into consideration. So the algorithm could avoid the missing of the global optimal solution possibly.

Numerical Examples
To evaluate the performance of CSA algorithm on the problem, several problems with different sizes were developed. We varied the values of several parameters of problem and heuristic. The problem parameters are described by the number of candidate facilities , the number of customers , and the number of products . The heuristic parameter is described by the coiling rate .
Here we use 20 instances with different sizes as benchmark problems. The stopping temperature value = 0.0001. The setting of different parameters and computational results of these instances are reported in Table 1. Each instance was computed 10 times and the optimal solution is chosen as the results (the OFV is the objective function value Φ).
It can be seen from Table 1 that the CSA could find the near optimal solutions for all instances with different sizes. It is clear that the bigger value of cooling rate could provide more benefit to the larger size problem.
The quality of the optimal solution found by the local search algorithm may be dependent on its initial solution. In order to test the model and the algorithm, we solved ten optimal solutions with different initial solutions for the largest instance, in which the system includes 18 candidate facility sites, 60 customers, and 5 products. Table 2 shows the data.
Computational results based on the ten different initial solutions are reported in Table 2. We can find that even the initial solutions are different, the objective functions values Φ of the optimal solutions obtained by the CSA algorithm are very close, and the gap between the maximal value (7592152.387) and the minimum value (7390387.970) is only The Scientific World Journal 7   2.73%. Compared with the initial solutions, the optimal objective values have 15.28%∼34.94% cost saving. So it can prove the proposed optimal model and algorithm is rational and available. Figure 1 shows the change rate of OFVs when the service level varied from 0.5 to 0.99. It is obvious that the OFVs increase when the service level rises. Also the greater the service level is, the faster the OFV increases. This case accords with the fact perfectly. Figure 2 illustrates that the relative change rate of OFV, compared with all the unit transportation costs, the unit inventory costs, and the demand deviations, changes with the same magnitude, respectively. We could find that the transportation cost has the most impact on the system, and when its magnitude of change varied from -50% to 50%, the change rate of the OFV exceeds 90%, which is beyond the impaction of inventory cost and service level on the OFV by far. Meanwhile, the change rate of the OFV is about 12% when the demand deviation varied.
It should be noted that the result is found under the condition of the weight factor of the inventory cost 1 twice as large as the weight factor of the transportation cost 2 . So we could consider the transportation cost is the most important factor in the logistics system cost. From the above analysis, we could find that the service level, the demand deviation, the inventory cost, and the transportation cost all could have influence on the total system cost, which are directly proportional to the cost. But the demand deviation and the transportation cost have more influence on the system cost in contrast with other factors. So from the managers' point of view, if they want to save their logistics system cost, they should pay more attention to selecting the most economic transportation mode and predict the demand more accurately.

Conclusions
In this paper, we analyzed the stochastic multiproduct capacitated facility location problem with stochastic demand and service-level constraint. A nonlinear mixed-integer programming model is proposed, in which the objective is to minimize the total system cost, including fixed facility setup cost, 8 The Scientific World Journal inventory cost, order cost, and transportation cost, under the precondition of satisfying the certain service level.
The CSA solution method is developed to solve the largescale problem, which is divided into two layers: the outer subalgorithm and the inner subalgorithm. The outer subalgorithm optimizes the facility location decision, and the inner subalgorithm optimizes the demand allocation based on the determined facility location decision. This method divides the problem into two layers to be solved; that is to say, the solution space in each iteration random search of the algorithm is also divided into smaller parts; thus the probability of obtaining the optimal solution of the problem in each iterative computation is raised.
We offer two important contributions to the SA application. First, we extend the breadth of applications by studying a new facility location problem, involving multiproduct and stochastic demand. Second, we developed a bilevel SA algorithm to solve the problem, in which computational performance was systematically evaluated under a variety of problem scenarios and control parameter settings.
Computational results show the CSA is a robust and practical approach for solving a multiple product problem, which could provide scientific guidance for the stochastic multiproduct facility location problem. Furthermore, the influences of the parameters (service level, demand deviation, unit inventory cost, and unit transportation cost) on the total system cost are investigated. We found that transportation cost and the demand deviation could have more influence than the other factors on the total cost (or the optimal decision).