A Hybrid Approach Using an Artificial Bee Algorithm with Mixed Integer Programming Applied to a Large-Scale Capacitated Facility Location Problem

1 Escuela de Ingenierı́a Informática, Pontificia Universidad Católica de Valparaı́so, Valparaı́so 2362807, Chile 2 Department of Engineering Science, University of Auckland, Auckland 1020, New Zealand 3 Instituto de Estadı́stica, Pontificia Universidad Católica de Valparaı́so, Valparaı́so 2362807, Chile 4 CIMFAV Facultad de Ingenierı́a, Universidad de Valparaı́so, Valparaı́so 2362735, Chile 5 Universidad Autónoma de Chile, Santiago 7500138, Chile 6 Departamento de Computación e Informática, Universidad de Playa Ancha, Valparaı́so 33449, Chile 7 Escuela de Ingenierı́a Industrial, Universidad Diego Portales, Santiago 8370109, Chile


Introduction
Heuristics and bioinspired techniques have become efficient and effective alternatives for researchers in solving several complex optimization problems.These types of techniques are able to provide satisfactory solutions for most of the applied problems within acceptable computational times.However, in spite of their effectiveness, these techniques are not able to reach the optimal solution or ensure its optimality for large-scale combinatorial optimization problems.In contrast, mathematical programming techniques, particularly the Mixed Integer Programming MIP , have been studied and developed by scholars over several decades with the main goal of obtaining optimal solutions to difficult problems using as little CPU time as possible.In this case, researchers must face the tradeoff between computational time and the quality of the result.For these reasons, the combination of metaheuristics and various mathematical approaches has become a well-studied area.Interested readers can find two recent and comprehensive works on the hybridization of stochastic techniques and mathematical programming MP approaches in 1, 2 .
Swarm Intelligence SI , as well as other mathematical programming techniques, has been applied successfully to several difficult combinatorial optimization problems see 3-6 for a range of applications .Additionally, as mentioned above, hybrid strategies have been developed to improve the effectiveness of these techniques.A more complete review is provided in Section 2.
The Artificial Bee algorithm BA is a relatively new approach in SI.Originally proposed in 2006 7 , it was mainly inspired by the foraging behavior of honeybees and has been applied to a range of problems in both combinatorial and functional optimizations with highly successful results.The BA has also been applied to large-scale problems 8 and has outperformed other well-known swarm-based algorithms, including Particle Swarm Optimization PSO and Differential Evolution DE .However, as with almost all stochastic approaches, the BA cannot ensure optimality even when the optimal solution is found.For small-and medium-size instances, this limitation is not highly relevant because heuristics techniques have empirically demonstrated their ability to achieve convergence in quite acceptable time; therefore, solutions provided by stochastic approaches are likely to be optimal or a very good approximation in large-scale instances.However, when large-scale optimization is considered, it is not possible to know how "good" the solution provided by the heuristic.
In this paper, we propose a hybrid algorithm of the BA and the MIP for solution of several large-scale instances of the well-known Capacitated Facility Location Problem CFLP .The CFLP is one of the most important problems for companies that distribute products to their customers.The problem consists of selecting specific sites at which to install plants, warehouses, and distribution centers while assigning customers to service facilities and interconnecting facilities using flow assignment decisions.This paper considers a twolevel supply chain in which a single plant serves a set of warehouses, which in turn serve a set of end customers or retailers.Figure 1 shows the basic configuration of our supply chain.Therefore, we aim to solve this problem by finding a set of locations that allow us to serve the entire set of customers in an optimal way.As Figure 1 shows, each customer or cluster is served only by one warehouse.
Despite its good performance in several optimization problems, the BA is not able to provide an optimal solution for large-scale problems.Furthermore, it is possible to become trapped in a local optimum.To bypass these drawbacks, we propose a hybrid algorithm that selects a subset of promising centers using the BA and subsequently solves the subproblem using a simple MIP algorithm.The BA guides the process while the MIP provides the optimal values of the simplified problems.One distinctive feature of our algorithm is that it solves the primary problem directly using the MIP algorithm, which is possible due to the reduction of the search space produced by the BA algorithm.
Although several works that have proposed various hybrid approaches to solve the CFLP and its uncapacitated UFLP version exist in the literature e.g., 9, 10 , we are unaware of any prior publications that use the Bees Algorithm BA to solve a large-scale CFLP.Moreover, we have found no articles that hybridize the BA with an MP approach in any optimization problem.Therefore, one contribution of this paper is the presentation of a performance analysis for the hybrid BA-MIP algorithm.A second contribution is the application of the BA to a large-scale problem that provides optimal solutions rather than only locally optimal solutions.
The remainder of this paper is organized as follows.Section 2 presents an overview of the CFLP and BA concepts.The hybrid BA-MIP algorithm is covered in Section 3, and a detailed explanation of the algorithm is also presented in this section.Section 4 begins with a brief description of the benchmarks used in this paper, and the experimental results are subsequently presented and discussed.Finally, Section 5 outlines selected conclusions.

Literature Review
This section presents a literature review.Section 2.1 discusses the mathematical model for CFLP and provides relevant background for certain approaches presented in the literature.Section 2.2 provides an overview of the BA and highlights its main features.

Capacitated Facility Location Problem
The CFLP in this work contains a set of warehouses that supply a set of customers that are uniformly distributed in a limited area.The model considers the installation cost i.e., the cost associated with opening a specific warehouse and transportation or assignation cost i.e., the cost related to transportation of a specific amount of products from a warehouse to a customer .The mathematical model for the CFLP is presented as follows: Equation 2.1 represents the total system cost.The first term denotes the fixed setup and operating cost for opening warehouses, and the second term indicates the daily transport cost between the warehouse and the customers.Equation 2.2 ensures that the customer demands are completely served by the system.Equation 2.3 ensures that the customers are assigned to the installed warehouses X i 1 .Equation 2.4 states that the summation of the demand μ of each customer served by a particular warehouse i must be less than or equal to a threshold I CAP i , which can be different for each warehouse.Finally, 2.5 states the integrality 0-1 for the variable X i and sets the range of the variable Y ij .This model is NP-hard because it is clearly an extension of the UFLP, which is known to be NP-hard.Additionally, this model is one of the most basic examples related to location research, and several comprehensive surveys on location theory can be found in 11, 12 .
The CFLP is well known in the operational research literature, and several authors have tackled this problem using different techniques e.g., Genetic Algorithms are implemented in 13 , and the Tabu Search TS algorithm is used in 14 .A comparison of the performance of these heuristics is provided in 15 .Additionally, mathematical-based approaches have been extensively developed to solve the CFLP, and algorithms based on Lagrangian relaxation represent the most common math-based approach 4, 8, 16 .In 17 , the authors develop a column generation strategy to obtain the exact solution for largescale instances.Other mathematical approaches are revised in 18, 19 .Moreover, mixed approaches using MP techniques and heuristics have been previously proposed.For instance, in 20 , the authors develop a Lagrangian-based heuristic LH that provides lower bounds to the problem, and the TS algorithm is subsequently used to find the upper bound of the problem.In this case, the TS is initialized using the primary information provided by LH.Additionally, in 21 , the authors combine Lagrangian relaxation with Ant Colony Optimization ACO .Although the CFLP is one of the most studied models in combinatorial optimization, to the best of our knowledge, no BA study has tackled this problem.

Bees Algorithm
The Bees Algorithm BA is a nature-inspired approach that was originally proposed by Pham et al. 7 to solve complex optimization problems.Subsequently, several authors have applied different versions of the BA to tackle a wide range of combinatorial and continuous optimization problems.The algorithm has been demonstrated to be highly competitive, especially when compared with other swarm/population-based approaches such as the PSO 22 or Genetic Algorithms 23 .Karaboga and Akay 24 provide a comprehensive literature review and a comparison of the most important swarm/population-based approaches.A complete survey of various BA applications is also provided by Karaboga and Akay in 25 .
In the following section, we present a brief description of the general structures of our BA.This description is mainly based on the descriptions provided by 7, 24-26 .
One of the most important characteristics of swarm intelligence is the ability to exchange relevant information among individuals.This feature allows the swarm to generate and develop collective knowledge.In the case of the BA, this information addresses the recognition of promising sources of food found by any individual insect of the swarm.Another relevant feature of the bee swarm is the ability to intensify the search in certain patches that have been identified as promising.The two main attributes of most heuristics used in optimization are exploration and exploitation.The first provides a fast and wide search throughout the search space which is usually too large .The second allows an intensive search of certain reduced search spaces neighborhoods that have been identified as "good-quality patches" during the exploration phase.In the BA case, the scout bees provide the exploration characteristics.The scout bees seek the search space in a usually random manner.Each scout bee visits a patch and evaluates it, and the scout bees have the ability to communicate the quality of the patch fitness to the unemployed bees.Depending on the attractiveness of each patch, the unemployed foragers will follow the scout bees to exploit the most promising patches.Once the patch is no longer attractive, a subset of bees will continue to search while the others wait for another promising patch.As in other heuristics, the balance between exploitation and exploration is a notably important issue for the BA.If we prioritize the exploration phase of the BA, it is likely to suffer from rather slow convergence.However, if we prioritize the exploitation phase of the BA, it is likely to become trapped in local minima.In most of the swarm optimization algorithms as well as other heuristics e.g., Tabu Search , the BA uses memory structures to influence the next population swarm .In this case, two main strategies provide information from former population to the new population.The first strategy uses experienced foragers, that is, bees with notably good fitness are included in the next population.The second strategy includes the best bee from a subset of the high-quality patches.In this manner, the use of various elements will determine the balance between exploration and exploitation.

Hybrid BA and LP Algorithm
In this paper, we present a BA that is hybridized with a MIP algorithm provided by GUROBI.Our algorithm contains the following distinctive features.

Variable Neighborhood Sizes
We use three different neighborhood sizes.It is important to note that we only refer to the "size" and not the "structure" of the neighborhood.More specifically, the neighborhood structure does not change during the algorithm execution.Figure 2 shows the three sizes and their function in the experiments.Figure 2 a shows a small-size move in which only one location is modified in this case, warehouse 4 is opened , and Figure 2 b presents an example of a medium-size move.In this case, two locations are modified warehouse 1 is closed, and warehouse 4 is opened .Figure 2 c shows a large-size move.
Using notation from 27 , we describe our variable neighborhoods as follow: let S denote any subset of open warehouses S ⊆ M ; the solution space may be subsequently defined as all such possible subsets.The total number of solutions in is 2 m − 1.To define the neighborhoods, 27 defines a distance function as follows: let S 1 , S 2 be any two solutions in ; the distance between them is defined by Therefore, the distance between one solution and another will increase when allocation of a specific customer in S 1 is different from the allocation of the same customer made in S 2 .An intuitive consequence of the above is that the distance between solutions S 1 , S 2 is zero if and only if S 1 S 2 .Figure 1 shows an example of the distance between the different solutions.For instance, in Figure 2 a , ρ S 1 , S 2 1; in Figure 2 b , ρ S 1 , S 2 2; and in Figure 2 c , ρ S 1 , S 2 4. The values of ρ used in our algorithm are shown at the end of Section 3. It is important to note that we do not implement a variable neighborhood strategy as in 27 .Instead, we define different neighborhood structures that are only applied for quite specific tasks during the execution of our algorithm.
i Neighborhood-based start strategy: the set of scout bees is initialized with a neighborhood-based strategy that uses the "large" movement to move from one solution to the other.This strategy allows us to use the LP solver more efficiently without compromising the algorithm performance.
ii Intensification procedure: we implement an intensification procedure to exploit the promising patches found by the bees.When a promising solution is found, the intensification procedure is triggered, and a local search using a rather "small" neighborhood movement is carried out.The intent behind this procedure is to find a local and hopefully global minimum as fast as possible, and this minimum is likely to be located near the promising solutions.
iii Roulette-wheel selection procedure: we implement a procedure based on the wellknown roulette wheel to select the elite bees and to assign the follower bees to particular elite or selected nonelite bees.This procedure allows for the inclusion of selected diversification mechanisms during the algorithm execution.
iv Use of experienced foragers: we include the information of the best solution in the form of "experienced foragers" in each iteration of the algorithm, which allows us to include certain historical information to improve the convergence of the algorithm.
The algorithm begins with a set of ns scout bees, which are initialized using the neighborhood-based start strategy.Other techniques have found that certain "warm" start solutions could be implemented as well.However, in our practical experience, our approach is an easier and faster method of initializing the set of scout bees mainly due to the efficiency mechanisms of the solver.Once the ns scout bees are generated, they are sorted according to their fitness.The fitness is calculated using the objective function 2.1 , which corresponds to an attractiveness measure in relation to the other scout bees.Note that the cost of each bee corresponds to the optimal solution for the subset of warehouses.This optimal value is provided by the MIP solver in less than two seconds on average.Therefore, this approach allows us to carry out a "global" optimization in different subspaces of the problem.Equations 3.1 and 3.2 state this attractiveness measure: where S is the set of scout bees and f x is our objective function 2.1 .The total cost of the set of scout bees is calculated n 3.1 and in 3.2 , and the fitness is calculated based on the fitness of a specific bee and the total cost of the swarm.Because we are aiming to optimize the fitness, our results must be normalized such that the lowest cost becomes the most attractive.Once the scout bees have been sorted, the e bees elite bees are selected using the roulette-wheel selection procedure.We note that, if there are important differences among the costs of the ns scout bees, then the algorithm will likely select the e most attractive ones.In the same way, the ne nonelite bees are selected from the remaining scout bees.Subsets e and ne are sorted, and the attractiveness of each bee is calculated as explained above.Next, a set F followers is generated and assigned to the both elite and nonelite bees.
Because the roulette-wheel selection procedure is used again, the elite bees are likely to recruit more followers than the less attractive nonelite bees.A local search is carried out at this stage using medium-size movements.Once the local search is completed, the algorithm selects the best bee from each patch.If a new best solution is reached, the intensification procedure is triggered.Following that procedure, the entire swarm is sorted based on fitness.Next, the elite and nonelite classification is executed again.In this step, a few random bees are added to the swarm to diversify the search.Additionally, the best solution is included via an experienced forager.The pseudocode for our hybrid BA-MIP algorithm is presented as follows Algorithm 1.

Parameter Tuning
To obtain a parameter set for the BA algorithm, we performed several tests using different values for each parameter of the algorithm, and these values are presented in Table 1.The tests were applied on one of the medium-size instances.The selected values are shown in bold.

Computational Experiments and Discussion
This section explains the experiments, certain benchmarks from the literature included to validate our algorithm, and the generation of a set of large-scale instances.The results obtained from our BA-MIP algorithm are compared with those from a simple local search strategy which uses the MIP to solve the allocation problem as well as those from both the MIP algorithm and the BA applied independently.A comprehensive analysis and discussion is outlined at the end of this section.

Experiments and Computational Results
In this subsection, we present the benchmarks applied for performance comparison and the computational results obtained for the hybrid algorithm.Finally, we show a summary of the principal results obtained.The computational experiments were performed on an Intel Core Duo processor CPU T2700, 2.33 GHz with 2 GB of RAM and Windows XP operating system.The BA algorithm was implemented in the JAVA programming language using NetBeans IDE, and the MIP algorithm was modeled with the GUROBI solver.
To validate the algorithm and measure its convergence, we chose a set of medium-size instances that have known optimal solutions, namely, capa, capb, and capc.These instances were obtained from Beasley's ORLibrary http://people.brunel.ac.uk/∼ mastjjb/jeb/info.html .Additionally, a set of large instances 300 warehouses and 1000 clients, 500 warehouses and 1000 clients, and 1000 warehouses and 1000 clients was created using the strategy provided in 8 .The set of customers and the set of warehouses are uniformly distributed over a plane of 10 × 10 distance units.The Euclidean distance between a customer i and a warehouse j corresponds to the transportation cost Y ij .The demand d j is calculated using a uniform distribution U 5, 35 .The I CAP i is calculated using U 100, 1600 , and we amplify the capacity of the warehouse to obtain harder instances.Finally, the fixed cost of warehouse i is calculated by F i U 0, 90 U 100, 110 I cap i /10.This expression takes into account the economies of scale 8 .As proposed in 8 , we generate three different classes of problems: 300 × 1000, 500 × 1000, and 1000 × 1000 warehouses × customers .To avoid any instance-dependent effects, we generated 10 different instances for each class.We also executed the algorithm 30 times for each instance to assess and avoid outlier performance.The results presented in this section correspond to the average values from these experiments.The MIP algorithm was executed only once per instance due to its deterministic behavior, and the results obtained from the MIP algorithm are presented in Table 2.The stop criterion used in the three first instances capa, capb, and capc was GAP ≈ 0%, that is, we forced the algorithm to find the known optimal solution.Because the optimal solution value is not known for our generated large-scale instances, the algorithm was aborted after 2000 seconds or GAP ≤ 1%, whichever occurred first.
Figures 3 and 4 show the convergence of two algorithms MIP and BA MIP for the instances 500 3 and 1000 4, respectively.
As shown in Figure 3 for the 500 3 instance, the two algorithms are both able to find notably good solutions within the time allotted.As expected, the MIP requires more CPU time than the BA-MIP hybrid algorithm to find the solution, and the solid line shows the LB of this instance.A remarkable feature of our hybrid approach is its "warm" initial solution.
As shown in Figure 3, the first iterations produce very good solutions that are quite close to the best solution reached by the MIP algorithm after the 2,000 seconds.We assume that this observation is due to our neighborhood-based start strategy.
Figure 4 shows a similar situation for the 1000 4 instance.The main difference in this case is the rapid convergence of the MIP algorithm.Despite this observation, the MIP is not able to find a better solution than that obtained by our BA-MIP algorithm.
Table 2 shows the average results obtained by the four algorithms for the three medium-size instances.The Opt column shows the optimal value produced for the respective instance.The Time column reports the CPU time required by the respective algorithm to reach the optimal value.The GAP column shows the difference between the average value and the optimal value as a percentage.A value of GAP 0% indicates that the optimal value was reached, and GAP > 0% otherwise.The maximum time available for each instance corresponds to the time in which the MIP algorithm was able to find the optimal value for that instance.
As shown in Table 2, the BA MIP and the LS MIP were able to find the optimal solution more quickly than the MIP solver.
These results demonstrate that the robustness of our algorithm highlights the speed of our hybrid approach, which could be a determining factor for notably large-scale instances in which the MIP is not able to find the optimal value within a reasonable CPU time.
The Time column in Table 3 shows the time in seconds in which the best value was found by the respective algorithm.The GAP column shows the difference between the average value and the lower bound value as a percentage.As stated previously, the maximum time available for each instance was 2000 seconds.The excellent performance of our hybrid approach is quite obvious; it outperforms the MIP approach in almost all cases, especially those instances in which the number of warehouses is greater than or equal to 500.This situation demonstrates the robustness of our algorithm as well as its reliability.It is also important to note that, in most of the cases, the best value obtained by the MIP algorithm was improved upon by our BA-MIP algorithm before 1500 seconds had elapsed, which implies a time saving of 25% compared with the mathematical approach.Moreover, in most cases, the best value provided by our hybrid algorithm was found near the time limit, which could be considered as a quite promising feature because it means that our algorithm is able to escape from local optima.Additionally, we note that, in the few cases in which the MIP reached a better value than our hybrid approach did, the difference is not greater than 1% of the GAP.The results obtained from the simple BA algorithm are surprisingly poor.We believe this situation may be partially due to the size of the search space as well as to a lack of precision in the parameter tuning.However, even when a careful parameter tuning process was conducted, the improvement is only marginal if we consider the current GAP between the LB for each instance and the UB obtained by BA.However, the LS-MIP obtains very good results for all instances.These results confirm that even rather simple heuristic strategies can be significantly improved when combined with mathematical programming approaches.However, we note that most of the best values were obtained early in the time elapsed, which may mean that the LS-MIP algorithm is not able to explore large spaces and quickly converges to local minima.Despite this weakness, the good performance shown by the LS-MIP algorithm encourages further exploration with different hybrid approaches.
Table 4 summarizes the results GAP, as a percentage obtained by our three algorithms the BA is not considered independent of any instance-dependent effects.It can be clearly observed that despite the effects provoked by particular instances, our algorithm presents the most desirable features, that is, a reduced mean, which suggests that the result is closer to the LB on average, and the small variance can be interpreted as a measure of robustness and reliability.

Conclusions and Future Work
In this paper, we have presented a hybrid algorithm based on the BA and the MIP and used it to solve the CFLP.The BA is applied primarily for the purpose of solving the location problem, that is, finding a subset from the available set of locations that satisfy the entire demand of the system.In contrast, the MIP is applied for the purpose of finding the optimal solution by considering a specific subset of locations, that is, solving the allocation problem.Certain considerations are important and must be highlighted in this approach.First, because the algorithm is able to find the optimal solution for each subset of selected warehouses, we are able to fairly compare those subsets.At the same time, our hybrid algorithm is able to find the optimal solution to medium-large problems, which is not possible with the use of common local search strategies.Additionally, obtaining the optimal solution allows us to compare the different strategies for warehouse selection and local search.
Second, compared with the widely used local search approaches that implement a swap in the assignation vector as a neighborhood move, our approach is able to visit more solutions because the entire tree corresponding to the feasible allocations is considered for each subset of warehouses in the search for the optimal solution of the subproblem by the MIP algorithm.
Third, our BA does not require extensive computational resources because most of the calculations are handled by the MIP solver, which is by far more efficient than many common heuristics implementations.
Our hybrid BA-MIP algorithm was able to find the optimal solution for a set of well-known large-scale instances found in the literature.Moreover, it outperformed both the BA and MIP approaches as applied separately.Compared with the state-of-the-art algorithms for location problems, our BA-MIP algorithm is highly competitive and reaches an optimal solution using less CPU time than both exact and stochastic approaches.Moreover, when the instances are notably large, the algorithm is able to find a better solution than that of the MIP approach in a fraction of the time.As widely reported in the literature over the last three decades, the combination of mathematical programming with heuristics and/or metaheuristics Matheuristics 1 provides robust and straightforward approaches to solving these types of problems.
In this paper, we have developed a hybrid algorithm using BA and LP using an approach that has not yet been reported.Due to the hybrid nature of the algorithm, certain changes were made in the structure of the BA heuristic.Use of the variable neighborhood "size," the neighborhood-based initialization procedure, the intensification procedure, and the experienced foragers are the most important features of our implementation.These distinctive components allow us to improve the performance of the simple BA when combined with MIP.
The hybridization of BA with other such mathematical programming approaches as interior point methods, column generation, and gradient-based algorithms shows promise as a potentially valuable research area.Additionally, the application of similar hybrid approaches to more complex large-scale problem will be an interesting future research line.

Figure 1 :
Figure 1: A two-level supply chain network configuration.

Figure 2 :
Figure 2: Three implemented neighborhood sizes: a small-size movement, b medium-size movement, and c large-size movement.

Figure 3 :Figure 4 :
Figure 3: Comparison of the convergence between the MIP solver and the BA-MIP algorithm instance 500 3 .

Table 1 :
Results of the parameter tuning process.

Table 2 :
Optimal DND obtained for both the ILM-PR and ILM-CR models.

Table 3 :
Obtained results per algorithm for all very large instances.

Table 4 :
Summary of the results obtained by the MIP, LS-MIP, and BA-MIP.