SIMULATION AND OPTIMIZATION OF ANT COLONY OPTIMIZATION ALGORITHM FOR THE STOCHASTIC UNCAPACITATED LOCATION-ALLOCATION PROBLEM

. This study proposes a novel methodology towards using ant colony optimization ( ACO ) with stochastic demand. In particular, an optimization- simulation-optimization approach is used to solve the Stochastic uncapacitated location-allocation problem with an unknown number of facilities, and an ob- jective of minimizing the ﬁxed and transportation costs. ACO is modeled using discrete event simulation to capture the randomness of customers’ demand, and its objective is to optimize the costs. On the other hand, the simulated ACO ’s parameters are also optimized to guarantee superior solutions. This approach’s performance is evaluated by comparing its solutions to the ones obtained using deterministic data. The results show that simulation was able to identify better facility allocations where the deterministic solutions would have been inadequate due to the real randomness of customers’ demands.

1. Introduction and literature review. In the continuous uncapacitated location allocation problem (also referred to as multi-facility weber problem (MWP)), the aim is to locate m new facilities that will serve n known demand points with an objective of minimizing the fixed costs of opening facilities and variable costs of transportation. The literature reports that most of the models developed for the facility location problem are very hard to solve to optimality and classified as NP-hard [15].
Due to the computational complexity of the problem, the literature offers limited exact methods for its solution. Earlier studies were able to solve for 2 facilities and 30 customers using branch and bound [10]. The problem size increased to cover 2 1216 JEAN-PAUL ARNAOUT, GEORGES ARNAOUT AND JOHN EL KHOURY to 100 facilities and 287 customers with Krau [9] that used a combination of B&B, global optimization, and Column generation.
Since exact solutions were not appropriate for practical scenarios where the problem size is usually large, researchers have focused more on heuristics and metaheuristics for the deterministic version of the problem. Bischoff and Klamroth [5] presented several heuristics to solve the capacitated version of the problem with rectilinear distances. Aras et al. [2] tackled the problem by partitioning the feasible region into a finite set of domains and solving the corresponding mixed-integer subproblems. Bischoff et al. [4] provided two heuristics for the MWP with barriers, and reported that their algorithms can attain solutions of reasonably sized multi-facility locationallocation problems with barriers, both with regard to computation time and solution quality. Brimberg et al. [6] compared between Genetic Algorithms (GA), TS, and Variable Neighborhoud Search (VNS), noted that the heuristics' solutions deteriorate when the number of facilities increases, with VNS performing better. Arnaout [3] introduced an Ant Colony Optimization (ACO) algorithm for the same problem with deterministic demand. ACO was compared to GA, VNS, TS, and the previous' superiority proven. For more on deterministic location-allocation literature, the reader can refer to Jabalameli and Ghaderi [8], Aras et al. [2], Salhi and Gamal [18], Gamal and Salhi [7], Liu and Xu [11], and Pasandideh and Niaki [17].
As shown above, most of the literature has focused on the deterministic version of the problem. Unfortunately, and in realistic scenarios, customers demand are always changing; i.e. customers have a stochastic demand and not a deterministic one. Logendran and Terrell [12] were the first to introduce the stochastic uncapacitated version and they modeled price sensitive stochastic demands with an objective of maximizing expected profits. Zhou [21] tackled the same problem with stochastic demand using an expected value model, chance-constrained programming and dependent-chance programming. In a later study, Zhou and Liu [20] generated new stochastic models for the capacitated variant of the problem. Finally, Mehdizadeh et al. [13] introduced a new hybrid algorithm for the capacitated version also.
More recent studies by Ozkisacik et al. [16] and Altinel et al. [1] addressed the probabilistic version of the problem, where they considered the customer locations to be randomly distributed. Their work differed from this research as the latter considers random customer demand but known locations.
In this study, we modify the ACO that was introduced for the deterministic problem in [3] to account for the stochastic nature of the problem, where customers demands follow uniform distributions. Following this, the stochastic ACO (referred to hereafter as ACO S ) is modeled using discrete event simulation, and then the algorithm's parameters are optimized using a combination of metaheuristic procedures in order to reach superior solutions. In other words, our novel approach optimizes the simulated ACO S parameters in order to optimize the solution of the problem at hand.
To the best of the author's knowledge, there does not exist published research that addresses the problem with unknown m and stochastic demand, nor research that models ant colony using discrete event simulation and uses optimization with simulation to determine the algorithm's parameters.
2. ACO S application to MWP. The objective function minimizes the TotalCosts (transportation and fixed costs) and is formulated as: where, n is the number of demand sites, d(i, j) is the distance between demand point i and facility j, k j is the index of new facilities (k j = 1 if facility j is open, 0 otherwise), and m = n is a suitable upper bound on the number of facilities as highlighted in [3]. Necessary parameters and variables are defined as follows: • x Fj , y Fj : decision variables that indicate respectively the x and y coordinates of new facility j. • x Di , y Di : respectively the x and y coordinates of demand point i.
• F: fixed cost of opening a new facility.
• T: transportation cost per unit of distance per unit demand.
• D i : stochastic demand / number of trips to demand point i.
As discussed earlier, the ACO S will be modeled using simulation. Furthermore, the MWP with m unknown can be solved in three ACO S stages: first, we determine the number of facilities m; second, we decide on initial temporary estimates for the facilities' coordinates (x Fj , y Fj ); and third, we determine the assignment of demand points n to the facilities m based on the latter's estimated coordinates, and we compute the final exact coordinates of each facility. The algorithm simulation logic is depicted in Figure 1, where at the initialization phase an entity is created to assign appropriate values to the pheromone trails of each of the three stages: τ I j , τ II ij , τ III ij . For instance, if we have 10 customer points (n = 10), then the following vector: S1 = [1 2 3 4 5 6 7 8 9 10] will represent the possible edges that an ant could move to. If an ant moves to edge [3] in (S 1 ), it means that m = 3 facilities in this ant's tour. The Pheromone trail (τ I j ) is defined for this stage to indicate the favorability of choosing the number of facilities from (S 1 ). In addition to the pheromone, we assist the algorithm with (η I j ) to solve the problem. The probability to select m is calculated as shown in (2) and (η I j ) is calculated using (3), which suggests the greedy heuristic of minimizing the number of facilities in order to minimize the fixed costs.
Following Equations (2) and (3) and as highlighted in Figure 1, the cumulative distribution function (CDF I (j)) is populated. Next, a random variable (rv) between 0 and 1 is generated based on which m is determined. Using this probabilistic approach in choosing m (versus solving for all its possible values) is better in terms of computational time, especially as m increases.
All m Facilities assigned coordinates?
Loop from j=1...n to assign the following:   Note that the m facilities should be located far from each other. Having two facilities next to each other is inefficient because they are uncapacitated; i.e. it makes more sense to replace them with one single facility. The coordinates' estimates are obtained by putting them equal to the ones of m demand points, where the latter are chosen as follows. The first facility (j = 1) coordinates are set to be equal to the ones of a randomly selected demand point. Now that one facility has coordinates (x F1 , y F1 ), (η II ij ) for the remaining facilities are computed following (4), where d(i, k) is the Euclidean distance between demand point i and facility k, with k referring to the predecessor facility (with an already assigned temporary coordinates). The probability to set the coordinates of facility j to be equal to the one of demand point i is calculated as shown in (5).
The Pheromone trail (τ II ij ) is defined for this stage to indicate the favorability of choosing the coordinates of demand point i for facility j's temporary coordinates. Similar to Stage 1, CDF II (i) is populated and rv is generated, based on which each facility is assigned temporary coordinates. The second stage output can be represented by a vector (S 2 ) that contains m entries representing the demand points from which we obtained facility j's temporary coordinates. For example, if we determined m = 4 from Stage 1, then S 2 =[2 10 4 6] indicates that the first facility's coordinates were set to be equal to the coordinates of the second demand point (i = 2), the second facility to i = 10's coordinates, the third to i = 4, and the forth facility to i = 6's coordinates.
2.3. Stage 3: Allocation of demand points to facilities. Now that the m facilities have initial coordinates, we determine in Stage 3 the assignment of demand points n to the facilities. A Pheromone trail (τ III ij ) is defined for this stage to indicate the favorability of assigning demand point i to facility j. The probability to assign i to j is calculated as shown in (6) and (η III ij ) is calculated using (7), which suggests the greedy heuristic of minimizing the distances between the demand points and the facilities.
The third stage output can be represented by a vector (S 3 ) that contains n entries representing the demand points. Each entry is populated with the assigned facility number (j = 1, , m). Note that this assignment was based on the temporary facilities' coordinates from Stage 2, and it will be used to determine the actual facilities' coordinates as highlighted in Section 2.4.

2.4.
Facilities coordinates, total costs, and pheromone update. After the third stage, each facility has a group of towns assigned to it. We use here a modification of the iterative procedure proposed by Weiszfeld [19] to come up with the facilities' coordinates. For every facility j, the steps below are implemented, where n j refers to the demand points that are assigned to facility j.
Step 1. Set the inital values of x Fj and y Fj according to 8: Step 2. For i = 1, ..., n j , {if(x Fj = x Di and y Fj = y Di ) → (x = x Di ; y = y Di )} else, exit for-loop and go to Step 4.
Step 3. For each demand point i, evaluate x and y as defined in 9: Step 4. If(x = x Fj and y = y Fj ), STOP Else, set(x Fj = x and y Fj = y ) and go to Step 2.
After all ants finish their paths, we update the pheromone amounts in each link locally by reducing the amounts due to evaporation. We also update the pheromone amounts in each link globally by increasing the ones constructed by the ant that produced the least T otalCosts. This is estimated according to 10, 11, and 12. T otalCosts Best , if arc(i, j) is used by best ant 0, Otherwise The innovative aspect about our solution methodology is the use of three successive trails to solve the problem: one for determining the number of facilities, another one for generating the facilities' temporary coordinates, and one more for assigning the demand points to the facilities. Note that D i follows uniform distribution as described in Section 4.

3.
Optimization of ACO S parameters. As can be seen from Section 2, ACO S has several parameters that will regulate the algorithm's performance. In particular, before simulating ACO S , we need to find the optimal values for the following parameters (along with their ranges): Ants : (5, 60); ρ : (0.01, 0.3); ϕ : (0.01, 0.3); α, β : (1, 2); τ I j , τ II ij , τ III ij : (1, 10). The ranges were adopted from Arnaout [3]. In order to decide on the optimal ACO S parameters, an optimization tool is used with simulation as shown in Figure 2. In particular, OptQuest from OpT ek Systems, Inc. is a computer software system that allows users to automatically search for optimal solutions to complex systems, and is integrated within Arena, the Simulation software that was used in this study. The optimization software based on its integrated heuristics, decides on values for the ACO S parameters. The latter is used in the simulation run and the T otalCosts are recorded. T otalCosts is input back into the optimization tool to decide on appropriate values of ACO S parameters for the next simulation. This loop is repeated for as many replications and simulations needed   to guarantee a minimum T otalCosts with a 95% Confidence Interval that is at most 5% from the average T otalCosts.
4. Computational tests. The proposed ACO S was modeled on Arena 13.0 Simulation Software from Rockwell Systems, running on Windows XP with a Pentium 4 processor at 2.33 GHz and 2 GB of RAM. The algorithms were run on small (5 and 20 stochastic demand points) and large (60 and 100 stochastic demand points) problems and each size was tested with 5 instances of the problem. The reason stochastic data was used for the customers' demand is to ensure a more realistic representation of the supply chain environment. The demands are stochastic following uniform distributions: where D di is the demand of customer i in the deterministic case as generated in Arnaout [3]. In particular, the data used for the demand (D i ) and coordinates of the demand points came from randomly generated values from a uniform distribution U [1, 100]. The reason uniform distributions were used is due to their high variances, ensuring that the presented heuristics are being tested under unfavorable conditions. Furthermore, the algorithms were compared under different dominance of F and T as follows. The selection of the fixed cost of opening a facility (F ) and the transportation (variable) cost (T ) values determines the level of dominance. That is, when the fixed and variable costs are balanced (denoted by F, T Balanced), they are set to F = $30, 000 and T = $30. When fixed costs are dominant (denoted by F Dominant), then the values of F and T are $30, 000 and $10 respectively. When the transportation costs are dominant (denoted as T Dominant), then F and T are set to $30, 000 and $50 respectively. Following this, a total of 60 problem instances were solved by each algorithm. This dominance approach was used to cover the different scenarios that might exist in practice, where the ratio of fixed and variable costs is highly dependent on the industry type and its location. Furthermore, the values of F and T were adopted from Arnaout [3] to have a fair comparison between ACO S and ACO D .

4.1.
Deterministic ACO (ACO D ). As one of the objectives of this study is to highlight the advantage of simulating stochastic demand versus using directly the deterministic ACO, the latter was implemented as follows.
The MWP with stochastic demand can be solved using the deterministic ACO (ACO D ) in Arnaout [3] that was modeled using Visual C++ and can only handle deterministic data. This is done by using the average of the demands' stochastic distributions (i.e. D di ) as input, running ACO D to output the number of facilities, their locations and allocations. Next, these outputs are simulated under the stochastic demand in order to generate the actual T otalCosts.

4.2.
Results. Table 1 and Figure 3 show the average and maximum T otalCosts for all problem structures and all algorithms. In particular, Table 1 shows the average and maximum T otalCosts attained, along with the number of new facilities recommended by both ACO D and ACO S . It can be seen from the latter that the number of facilities is different when the simulation was used to model ACO and capture the stochastic demand, leading to lower T otalCosts. It can be seen that ACO S performed better than ACO D in all the 60 instances. This is due to the fact that a more suitable number and allocation of facilities was attained after simulating the randomness of the demand.  5  3  91418  93468  2  77430  77658  18  20  20  4  245229  248614  3  220412  220894  11  13  60  7  627527  644356  6  502022  503403  25  28  100  10  945643  971370  8  716396  724903  32  34   F, T Balanced   5  4  120563  123091  3  99451  99746  21  23  20  8  388249  400121  10  363961  367935  7  9  60  18  1108033  1138302  23  872467  27882405  27  29  100  29  1833518  1884079  38  1378585  1385353  33  36   T Dominant   5  4  127561  131693  4  123327  123328  3  7  20  10  460391  472180  10  427723  434819  8  9  60  26  1459744  1499851  33  1186784  1199881  23  25  100  42  2534624  2604796  55  1964825  1973330  29  32 Even though ACO S outperformed ACO D in all problem instances, and as the latter's solutions appeared to be close to ACO S (see Figure 3), the percentage deviation (δ) of ACO D from ACO S was used as a measure of performance for each problem instance. That is: Figure 4 shows the average values of δ for the three dominance combinations for all problems. It can be clearly seen that ACO D 's deviation from ACO S increased with the problem size. Another important note about the δ values is that the average and maximum deviations of ACO D for all problem instances were 20% and 36% respectively. 4.3. Computational times. As highlighted earlier, ACO D was modeled using Visual C++ for deterministic demand, while ACO S was modeled using discrete event simulation in order to account for the stochastic demand. Consequently,      comparing the computational times for the two ACO approaches would not be a reasonable assessment. Furthermore, the solution of ACO D is simulated in order to generate the T otalCosts, and ACO S 's parameters are optimized to generate the T otalCosts; i.e., both will require an additional computational time beyond the algorithms' scope.
Having said this, and to give the reader a sense of the computational performance, Table 2 depicts the average times (minutes) for the algorithms (excluding ACO D 's simulation and ACO S 's optimization). As can be seen, both ACO D and ACO S have similar performance; this is logical because both follow similar algorithmic steps. Furthermore, it can be noted from Arnaout [3] that ACO D requires less time to converge to more superior solutions than Genetic Algorithms (GA); subsequently, the same would hold for ACO S .

5.
Conclusions. While lots of research has tackled the deterministic continuous facility location-allocation problem, few works addressed the stochastic version of the problem, and no previous published studies were found on the problem with unknown number of facilities with stochastic demand, nor research that modeled ACO using discrete event simulation and used optimization with simulation to determine the ACO parameters. In this paper, we have modeled a previously introduced ACO ([3]) using discrete event simulation to account for the randomness of customers' demands for the Stochastic Euclidean facility location-allocation problem with unknown number of facilities. The differences between this study and its predecessor ( [3]) are as follows: 1. In Arnaout [3], ACO was modeled using Visual C++. On the other hand, in this study ACO was modeled using simulation. Up to the author's knowledge, no literature exists that models Ant Colony Optimization using simulation; in particular, discrete event simulation. This allowed for the introduction of stochastic demand in this study.

The computational results indicated a significant reduction in costs when
ACO was simulated in comparison to its previous deterministic implementation in Arnaout [3]. In particular, a reduction of up to 36% in total costs was attained. Furthermore, better facilities allocations were achieved when the randomness of data was taken into account. 3. One of the challenges of ACO is finding suitable parameters. In Arnaout [3], design of experiments was used to find suitable parameters. In this study, an optimization approach was used with Simulation to find the optimal ACO parameters. This by itself is another addition to the body of knowledge. In particular, optimization decides on values for the ACO S parameters. The latter is used in the simulation run and the T otalCosts are recorded. T otalCosts is input back into the optimization tool to decide on appropriate values of ACO S parameters for the next simulation. This loop is repeated for as many replications and simulations needed in order to guarantee a minimum T otalCosts.
In other words, we are optimizing an optimization method (ACO S ) using Simulation.