Abstract
We present a hybridization of two different approaches applied to the well-known Capacitated Facility Location Problem (CFLP). The Artificial Bee algorithm (BA) is used to select a promising subset of locations (warehouses) which are solely included in the Mixed Integer Programming (MIP) model. Next, the algorithm solves the subproblem by considering the entire set of customers. The hybrid implementation allows us to bypass certain inherited weaknesses of each algorithm, which means that we are able to find an optimal solution in an acceptable computational time. In this paper we demonstrate that BA can be significantly improved by use of the MIP algorithm. At the same time, our hybrid implementation allows the MIP algorithm to reach the optimal solution in a considerably shorter time than is needed to solve the model using the entire dataset directly within the model. Our hybrid approach outperforms the results obtained by each technique separately. It is able to find the optimal solution in a shorter time than each technique on its own, and the results are highly competitive with the state-of-the-art in large-scale optimization. Furthermore, according to our results, combining the BA with a mathematical programming approach appears to be an interesting research area in combinatorial optimization.
1. 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 meta-heuristics 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 two-level 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.
2. 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.
2.1. 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: s.t.: 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 ( = 1). Equation (2.4) states that the summation of the demand of each customer served by a particular warehouse must be less than or equal to a threshold , which can be different for each warehouse. Finally, (2.5) states the integrality (0-1) for the variable and sets the range of the variable . 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 large-scale 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.
2.2. 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.
3. 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.
3.1. 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.
(a)
(b)
(c)
Using notation from [27], we describe our variable neighborhoods as follow: let S denote any subset of open warehouses ; the solution space may be subsequently defined as all such possible subsets. The total number of solutions in is 2m − 1. To define the neighborhoods, [27] defines a distance function as follows: let , 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 is different from the allocation of the same customer made in . An intuitive consequence of the above is that the distance between solutions , is zero if and only if = . Figure 1 shows an example of the distance between the different solutions. For instance, in Figure 2(a), ; in Figure 2(b), ; and in Figure 2(c), . 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 well-known 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 is the set of scout bees and is our objective function (2.1). The total cost of the set of scout bees is calculated (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 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 most attractive ones. In the same way, the ne (nonelite) bees are selected from the remaining scout bees. Subsets 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.
| ||||||||||||||||||||||||||||||||||||||||||||
3.2. 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.
4. 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.
4.1. 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 and a warehouse corresponds to the transportation cost . The demand is calculated using a uniform distribution [5, 35]. The is calculated using U [100, 1600], and we amplify the capacity of the warehouse to obtain harder instances. Finally, the fixed cost of warehouse is calculated by . 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 , 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 , 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.
5. 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.