1. INTRODUCTION
The hub-and-spoke network was first introduced by Goldman in Reference Goldman1969, followed by O'Kelly (Reference O'Kelly1986a , Reference O'Kelly1986b ). O'Kelly also proposed the first mathematical formulation for this problem as a quadratic integer programming model in Reference O'Kelly1987. Afterward, many researchers have been working on various categories of the hub location problems (HLPs), including both theoretical and practical issues. The major reason for this deep interest is the necessity and implication of modern infrastructures for transportation and telecommunication systems (Campbell et al., Reference Campbell, Ernst, Krishnamoorthy, Drezner and Hamacher2002).
The classical HLP can be described as follows: let G be a complete graph of vertices and edges; G = (V; E), where V = {v 1, v 2, … , v n } is the set of all vertices and E is the set of all edges. The elements of V represent all potential vertices for locating hubs. The distance between two vertices i and j is represented by d ij and w ij and shows the flow between these vertices. Distances satisfy the triangle inequality. The goal of the problem is assigning hubs to some vertices to minimize the total flow cost of the network.
In HLP, vertices locating hubs are called hub vertices, and others are called spoke vertices. Consequently, a hub-level network is formed by connecting pairs of hub vertices by a hub edge, Each spoke vertex will be allocated eventually to these hub vertices by spoke edge links.
Based on real-world problems, this paper considers a capacitated model of the multiple-allocation cluster hub location problem for the first time. Moreover, the capacity constraints are considered on the incoming and outgoing flow of vertices and the amount of flow traversing the transportation lines. These assumptions represent a more realistic model that is able to handle the congestion-related issues in the resulting network. Two major applications of cluster HLP are telecommunication and cargo systems. Finding the location of hub facilities is the major decision in these networks because it affects the total transportation cost and also the total service time.
The presented problem, the capacitated multiple-allocation cluster hub location problem (CMACHLP), is NP-hard since its special case is the classical HLP belonging to the class of NP-hard problems. Therefore, we propose a novel hybrid algorithm based on simulated annealing (SA) and ant colony optimization (ACO) to solve the problem. The proposed algorithm employs a new version of ACO that is combined with a powerful ternary local search procedure to improve the quality of solutions. Finally, the SA algorithm is used as a framework for the proposed hybrid heuristic to prevent falling into a local optimum.
The rest of paper is organized as follows. Section 2 provides a brief literature review on various HLP and their solution approaches. Mathematical formulation is presented in Section 3. Section 4 gives an overview of SA and ACO algorithms and discusses the principles of the proposed algorithm used to solve the CMACHLP. Section 5 shows the performance of the algorithms on a set of benchmark problems. Finally, Section 6 concludes the paper.
2. LITERATURE REVIEW
A brief description of different types of HLP is considered in this section, and then a literature survey is presented of some heuristic algorithms that are proposed in this area.
2.1. HLP
The HLP has extensive variants, including the following: HLP: minimizing of the total cost is considered; the total cost consists of the cost of establishing hub facilities (vertices and/or edges) plus transportation costs; the p-hub median problem (pHMP) in which minimizing the transportation cost is considered by locating the given number of hubs, p; the p-hub center problem (pHCP): in which given the number of hubs, p, the objective is to minimize the maximum cost for each origin–destination pair (pHCP-I), or for movement on each single link (pHCP-II) or for movements between a hub and origin/destination (pHCP-III); the hub covering problem (HCP): in which minimizing the total transportation cost is considered. Transportation costs should not exceed a certain threshold on origin–destination pairs (HCP-I), or on each link (HCP-II) or between hub and origin/destination (HCP-III); and the hub arc location problem (HALP): in which minimizing the total transportation cost is considered by locating given the number of hub edges, q (Gelareh, Reference Gelareh2008).
Regarding HLP, two major decisions are important: finding vertices to locate hubs and allocation of spoke vertices to each hub. In general there are two categories of HLPs. In the first category, a single hub vertex should receive (and deliver) total flow of an origin (or destination) vertex (hub or spoke) it is assigned to, known as single allocation (SAl). In the second category, spoke vertices can be processed by more than one hub. This category is called multiple allocation (MAl) and obviously is more applicable than SAl.
HLPs involve different varieties based on different types of capacities: capacity on the amount of flow from spoke to hub vertices, capacity on the all through traffic, and capacity on the hub edges. However, usually hub and edge capacities are considered separately.
As mentioned earlier, the idea of HLPs is attributed to Goldman in Reference Goldman1969. He considered locating n centers (processing facilities) in a network, so as to minimize the total transportation cost. He proposed two basic models. In the first one, materials are sent from origins to destinations and pass through a processing as well as a collect-and-dispatch center. The second model assumed flows can be processed in more than one center. The problem optimally finds p centers and assigns the flows to the center(s), based on the transportation cost.
O'Kelly (Reference O'Kelly1986a , Reference O'Kelly1986b ) is the pioneer of the HLP field of research, and after his scientific efforts, many researchers throughout the world worked on various issues of this field. He also made the first study on the discrete HLP in Reference O'Kelly1987, where he proposed the first formulation as a quadratic integer program for HLPs. He presented the formulation for the SAl p-hub median problem (SApHMP), which is also known as uncapacitated SApHMP.
The first study and model of the MAl problem is attributed to Campbell in Reference Campbell1992. He formulated the MAl p-hub median problem (MApHMP) as a linear integer program. In Reference Ernst and Krishnamoorthy1998, Ernst and Krishnmoorthy considered a new modeling of uncapacitated MApHMP based on their earlier work on the SAl problem. Sung and Jin (Reference Sung and Jin2001) formulated the uncapacitated cluster HLP, which can be seen as an extension to the UMApHMP. In Reference Wagner2007, Wagner improved the formulation of the uncapacitated cluster HLP proposed by Sung and Jin (Reference Sung and Jin2001). He showed that the location problem-relaxation of the proposed model provides very tight lower bounds.
Capacitated studies against uncapacitated ones are scare. Yaman in Reference Yaman2005 considered uncapacitated HLP with modular arc capacities (HLM) and proposed a model for a two-level telecommunication network with modular capacity. In his research, integer amounts of capacities are assumed on the arcs between hub pairs and terminals and their hubs. Marin (Reference Marin2005) considered the multiple-allocation HLP and described a new capacitated model in which the capacity is considered as an upper bound on the total flow coming directly from the origins. Furthermore, he assumed that the flow between a given origin–destination pair can be split into several routes. He presented tight integer linear programming formulations for the problem, along with some useful properties of the optimal solutions, which can be used to speed up the resolution. The presented approach could solve medium-size problems very efficiently and outperformed other approaches in the literature.
Correia et al. (Reference Correia, Nickel and Saldanha-da-Gama2010) revisited the classical formulation for the capacitated single-allocation HLP. They completed the formulation by the introduction of a new set of constraints, which also turned out to be computationally very effective as a cut. Contreras et al. (Reference Contreras, Díaz and Fernández2011) presented a branch-and-price algorithm for the capacitated single-allocation HLP in which Lagrangian relaxation is used to obtain tight lower bounds of the restricted master problem.
Contreras et al. (Reference Contreras, Cordeau and Laporte2012) presented an extension of the classical capacitated multiple-allocation HLP in which the amount of capacity installed at the hubs is part of the decision process. They proposed a Benders’ decomposition algorithm to solve the splittable and nonsplittable demand cases of the problem. Chen (Reference Chen, Huang, Bevilacqua, Figueroa and Premaratne2013) studied the multiple-allocation HLP in which the number of hubs and the hub capacity are decision variables. He proposed two heuristics based on the threshold accepting method to solve the problem. Correia et al. (Reference Correia, Nickel and Saldanha-da-Gama2014) extended the classical capacitated single-allocation HLP by considering that multiple products are shipped through the network. They presented several sets of valid inequalities for strengthening the linear relaxation of the proposed models.
To the best of our knowledge, there is not a significant work that addresses the capacitated multiple-allocation cluster HLP, despite its wide application in the transportation sector. The mentioned gap is the main focus of this paper. Furthermore, another contribution of this study is to present an efficient solution algorithm for this hard class of HLPs, which will be discussed in more detail later.
2.2. Solution approaches
The HLPs are known as NP-hard problems (Kara, Reference Kara1999), and this means finding optimal solution for the real-world size problems in desirable time and adjustable cost is impossible. Therefore, considering heuristics and metaheuristics is a major part of research in this area.
2.2.1. Genetic algorithm
Abdinnour-Helm and Venkataramanan proposed the first genetic algorithm (GA) for uncapacitated HLP in Reference Abdinnour-Helm and Venkataramanan1998. About one decade latter Cunha and Silva proposed another GA in Reference Cunha and Silva2007. The problem consists of determining the number of hubs, their locations, and the assignment of the spokes to the hubs, aiming to minimize the total cost. Presented GA incorporated an efficient local improvement procedure that was applied to each generated individual of the population. They embedded a fast and easy approach to allow nonimproving solutions to be accepted with a given probability. A heuristic based on GA was proposed for solving the capacitated SApHMP by Stanimirovic was Reference Stanimirovic2010. Naeem and Ombuki-Berman proposed an effective GA for solving uncapacitated single-allocation hub location problem (USAHLP) in Reference Naeem and Ombuki-Berman2010. Marić et al. in Reference Marić, Stanimirović and Stanojević2013 presented a memetic algorithm for solving the USAHLP. They proposed two efficient local search heuristics in the frame of an evolutionary algorithm in order to improve both the location and the allocation parts of the problem. Computational results showed that the proposed approach is superior over existing heuristic approaches for solving the USAHLP. Damgacioglu et al. in Reference Damgacioglu, Dinler, Ozdemirel and Iyigun2015 studied the uncapacitated single-allocation planar HLP in which capacities of hubs are infinite, they can be located anywhere on the plane, and are fully connected. Authors proposed a mathematical formulation and a GA to solve the problem in a reasonable time.
2.2.2. SA
Ernst and Krishnamoorthy proposed a SA for the capacitated SAl hub location problem (CSAHLP) in 1999. They applied the results as the upper bound in a branch-and-bound. Jabalameli et al. addressed the single-allocation maximal covering HLP in Reference Jabalameli, Barzinpour, Saboury and Ghaffari-Nasab2012. They proposed a SA-based heuristic to solve it. Rasoulinejad et al. in Reference Rasoulinejad, Bashiri and Mehrbanfar2013 applied a clustering-based SA approach for solving uncapacitated single-allocation p-hub location problem. Parvaresh et al. formulated MApHMP under intentional disruptions in Reference Parvaresh, Hashemi Golpayegany, Moattar Husseini and Karimi2013. They applied two algorithms based on SA to solve the defined problem.
2.2.3. Tabu search
Several tabu search (TS) algorithms were also proposed in the literature. The first TS was proposed for SApHMP by Klincewicz in Reference Klincewicz1992. In Reference Skorin-Kapov and Skorin-Kapov1994, Skorin-Kapov and Skorin-Kapov developed another TS based on 1-exchange for the problem of locating p interacting hub facilities among n interacting nodes in a network. Xu et al. considered a multiple hub location routing problem and designed a two-stage TS algorithm to solve it in Reference Xu, Hu, Xuan, Lin, Liu, Zhang and Guan2009. Calık et al. proposed a heuristic based on TS for solving a single-allocation hub covering problem in Reference Calık, Alumur, Kara and Karasan2009.
2.2.4. Particle swarm optimization
The first particle swarm optimization was proposed for USAHLP by Bailey et al. in Reference Bailey, Ornbuki-Berrnan and Asobiela2013. Geramianfar et al. considered a multiobjective hub covering location problem under congestion in Reference Geramianfar, Pakzad, Golhashem and Moghaddam2013. SI was used to solve the problem, and the performance of the proposed algorithm was compared against particle swarm optimization and NSGA-II. The computational results showed that SA outperforms the other two alternative methods in terms of quality.
2.2.5. ACO
ACO is a strong metaheuristic that has been successfully applied to a range of problems, especially in scheduling and vehicle routing problems. Randall in Reference Randall2008 considered the capacitated single-allocation HLP and assumed the number of hubs is not fixed and must be determined by the solver system. He outlined four possible solution construction techniques within an ACO framework and applied to the CSAHLP. Sun (Reference Sun2011) made the same research on the capacitated SAl p-hub median problem. He presented a mathematical formulation for the problem and developed a solution method based on ACO. Ma and Thing (Reference Ma and Thing2008) considered uncapacitated HMP with multiple allocation (UMApHMP). They modified the solution construction rule of the canonical ACO algorithm, which involves the locations of p hubs and the allocations of nonhub nodes to those selected hubs.
2.2.6. Other approaches
A hybrid of GA and TS is proposed by Abdinnour-Helm in Reference Abdinnour-Helm1998 for USAHLP. They showed that the hybrid algorithm provides better performance compared to applying GA alone. In 1999, Ernst and Krishnamoorthy presented an algorithm based on a random descent heuristic, which is one of the best algorithms proposed for CSAHLP up to now. Another hybrid algorithm in the literature is attributed to Chen in Reference Chen2007 for USAHLP. He proposed two approaches in order to determine the upper bound for the number of hubs along with a hybrid heuristic based on the SA method, tabu list, and improvement procedures. The proposed heuristic was capable of obtaining optimal solutions for all small-scaled problems very efficiently. This algorithm outperformed the GA of Topcuoglu et al. (Reference Topcuoglu, Corut, Ermis and Yilmaz2005).
Puerto et al. (Reference Puerto, Ramos and Rodríguez-Chía2013) presented a branch-and-bound-and-cut based algorithm to solve the single-allocation ordered median HLP. The proposed approach performs well for small- to medium-sized problems. Saboury et al. (Reference Saboury, Ghaffari-Nasab, Barzinpour and Jabalameli2013) considered a two-layered hub-and-spoke network in which both access and backbone networks are interconnected. Authors developed two hybrid heuristics based on SA and TS to solve the problem. Sender and Clausen (Reference Sender and Clausen2013) proposed different heuristic approaches based on local improvements to solve a capacitated multiple-allocation HLP derived from a practical application in the network design of German wagonload traffic.
He et al. (Reference He, Wu, Zhang and Liang2015) developed an improved heuristic combining branch-and-bound, Lagrangian relaxation, and linear programming relaxation to solve the intermodal HLP. The proposed heuristic generates a population of initial feasible solutions using the branch-and-bound and Lagrangian relaxation methods and creates a linear-relaxed solution using the linear programming relaxation method. These feasible and linear-relaxed solutions are combined to fix a number of hub location variables, which leads to a number of restricted hub location subproblems. Then, these restricted subproblems are iteratively solved by the branch-and-bound method to improve the solution quality.
For more reviews about various issues of HLPs, we refer readers to the two latest reviews by Alumur and Kara (Reference Alumur and Kara2008) and Farahani et al. (Reference Farahani, Hekmatfar, Arabani and Nikbakhsh2013), where more details about HLPs are discussed. A fast review of the literature shows some variants like uncapacitated and single allocation have received most attention by researchers so far.
3. MATHEMATICAL MODELING
In this section, we propose a new mathematical model for CMACHLP. This mathematical model minimizes two major costs: transportation cost and establishment cost. The main assumptions are the following:
-
1. The amount of flows and transportation costs are constant.
-
2. The number of clusters is predetermined and is equal to p.
-
3. Every vertex is only located in one cluster and the vertex distribution is predetermined.
-
4. Every cluster can only have one hub.
-
5. The hub's establishment cost is constant.
-
6. Vertices can directly be connected to each other (direct service).
-
7. A transportation path is a direct edge between origin and destination or an edge through two various hubs.
-
8. Each hub has a predetermined capacity, which it cannot exceed.
-
9. Each transportation path has a predetermined capacity, which it cannot exceed.
The notations that we use are the following:
-
C r is the set of vertices that belongs to cluster r;
-
C(i) is the set of vertices that is included in the cluster, which corresponds to vertex i;
-
c ij is the transportation cost per unit for the direct edge from vertex i to vertex j;
-
d ij is the amount of flows from vertex i to vertex j;
-
f i is the predetermined establishment cost of one hub in vertex i;
-
α is the discount factor related to transfer between hubs;
-
q k is the capacity of hub k; and
-
q ij is the capacity of the transportation path from vertex i to vertex j.
Our decision variables are the following:
-
y i is a binary variable that is equal to 1 if one hub is located in vertex i;
-
z ij is a binary variable that is equal to 1 if the direct edge from vertex i to vertex j is used to flow transportation; and
-
x ijkm is a binary variable that is equal to 1 if the transportation from vertex i to vertex j is flowed through hubs k and m, respectively.
Figure 1 demonstrates the parameters and decision variables of the problem in a typical example. Now we can describe the mathematical model as follows:
$$\hbox{min}\; \mathop \sum \limits_{i = 1}^n \mathop \sum \limits_{\,j = 1}^n \left( {{\tilde c}_{ij}.z_{ij} + \mathop \sum \limits_{k \in C\lpar i \rpar } \mathop \sum \limits_{m \in C\lpar j \rpar } {\tilde c}_{ijkm}.x_{ijkm}{\rm \;}} \right) + \mathop \sum \limits_{i = 1}^n \,f_i.y_i \hbox{,} $$
such that
$$\mathop \sum \limits_{k \in C_r} y_k = 1 \quad r = 1\hbox{,} \; .\,.\,. \; \hbox{,} \;p$$
$$z_{ij} + \mathop \sum \limits_{k \in C\lpar i \rpar } \mathop \sum \limits_{m \in C\lpar \,j \rpar } x_{ijkm} = 1 \quad i \hbox{,} \; j = 1 \hbox{,}\, .\,.\,. \, \hbox{,} \; n \hbox{,} \; i \ne j$$
$$y_k \ge \mathop \sum \limits_{m \in C\lpar \,j \rpar } x_{ijkm}\quad k \hbox{,} \; j = 1 \hbox{,}\; .\,.\,.\, \hbox{,} \; n \hbox{,}\; i \in C\lpar k \rpar $$
$$y_m \ge \mathop \sum \limits_{k \in C\lpar i \rpar } x_{ijkm} \quad i \hbox{,}\, m = 1 \hbox{,}\; .\,.\,. \, \hbox{,} \; n \hbox{,} \; j \in C\lpar m \rpar $$
$$\eqalign{\mathop \sum \limits_{m \in C\backslash C\lpar k \rpar }& \mathop \sum \limits_{i \in C\lpar k \rpar } \mathop \sum \limits_{\,j \in C\lpar m \rpar } d_{ij}. x_{ijkm} + \mathop \sum \limits_{m \in C\backslash C\lpar k \rpar } \mathop \sum \limits_{i \in C\lpar m \rpar } \mathop \sum \limits_{\,j \in C\lpar k \rpar } d_{ij}.\, x_{ijmk} \cr & \quad \quad \quad + \mathop \sum \limits_{i \in C\lpar k \rpar } \mathop \sum \limits_{\,j \in C\lpar k \rpar } d_{ij}.x_{ijkk} \le q_k.y_k \quad k = 1 ,\; .\,.\,.\ , \; n }$$
$$\eqalign{d_{ij}.z_{ij} & + \mathop \sum \limits_{k \in C(i)} \mathop \sum \limits_{m \in C\lpar j \rpar } \; d_{km}.x_{kmij} + \mathop \sum \limits_{k = 1}^n \mathop \sum \limits_{m \in C\lpar k \rpar } \; d_{ik}.x_{ikjm} \cr & \quad\quad\quad\quad\quad + \mathop \sum \limits_{m = 1}^n \mathop \sum \limits_{k \in C\lpar m \rpar } \; d_{mj}.x_{mjki} \le q_{ij} \quad i\hbox{,} \; j = 1 ,\, .\,.\, .\ , \; n}$$
$$z_{ij} \in \lcub {0 \hbox{,} \,1} \rcub \quad i \hbox{,} \; j = 1 \hbox{,}\; .\,.\,.\; \hbox{,} \; n\; $$
$$x_{ijkm} \in \lcub {0 \hbox{,}\, 1} \rcub \quad i \hbox{,} \;j = 1 \hbox{,} \, .\,.\,.\ \hbox{,} \; n \hbox{,} \; k \in C\lpar i \rpar \hbox{,} \; m \in C\lpar\, j \rpar $$
$$y_i \in \lcub {0 \hbox{,}\,1} \rcub \quad i = 1 ,\; .\,.\,.\; , \; n$$
where
$\tilde c_{ijkm} = d_{ij}.(c_{ik} + {\rm \alpha}. c_{km} + c_{mj})$
and
$\tilde c_{ij} = d_{ij}.c_{ij}$
.
An illustrated example of the problem parameters and decision variables.
The objective function (1) consists of two terms and minimizes the transportation and establishment costs. Constraint (2) ensures that only one hub is assigned to each cluster. Constraint (3) implies that the path between two vertices is a direct edge or indirect by an edge-connected intermediate hubs. Constraints (4) and (5) show that the transition through the hub vertex is allowable only and just only if there is a hub in the vertex. Constraint (6) is related to the capacity constraints on the amount of flow arriving to and leaving from the hub vertices. Constraint (7) indicates the capacity constraints of transportation paths. Finally, constraints (8) to (10) impose integrality restrictions on the binary variables.
It can be noted that such mathematical models of the problems usually lead to very large and hard-to-solve problems, in many cases to NP-hard problems. This is the indispensable feature of many real-life applications.
4. HYBRID ALGORITHM BASED ON SA AND ACO
As mentioned earlier, the HLP is known as NP-hard problems (Kara, Reference Kara1999), and metaheuristics are a good choice to solve this group of problems. One of the most powerful metaheuristic to solve hard optimization problems, especially scheduling and vehicle routing problems, is the ACO technique. Based on the previous research, ACO can work effectively to find near-optimal solutions in scheduling problems if it benefits from an effective initialization phase and also applies an appropriate local search. Moreover, SA is a metaheuristic, which is able to prevent falling into the local optimums. These facts have stimulated us to develop a new hybrid SA-ACO algorithm to solve CMACHLP.
SA is a stochastic local search that was presented by Kirkpatrick et al. in Reference Kirkpatrick, Gelatt and Vecchi1983. This algorithm searches solution space by using a moving procedure from the current solution to a new solution in its neighborhood. The objective function value is evaluated in each move as a criterion for accepting the new solution as the current solution. A new solution that has a bad quality (higher objective function value in comparison with the current solution in a minimization problem) may be accepted as well with a specific probability. We use the main idea of the SA algorithm as the framework for our hybrid SA-ACO.
The behavior that provided the inspiration for ACO is the foraging behavior of ants, and in particular, how ants can find the shortest paths between food sources and their nest. When searching for food, ants initially explore the area surrounding their nest in a random manner. While moving, ants leave a chemical pheromone trail on the ground. Ants can smell pheromones. When choosing their way, they tend to choose, in all probability, paths marked by strong pheromone concentrations. As soon as an ant finds a food source, it evaluates the quantity and the quality of the food and carries some of it back to the nest. During the return trip, the quantity of the pheromone that an ant leaves on the ground may depend on the quantity and quality of the food. The pheromone trails will guide other ants to the food source. It has been shown that the indirect communication between the ants via pheromone trails, known as stigmergy, enables them to find the shortest paths between their nest and food sources (Blum, Reference Blum2005). The ACO metaheuristic, as we know it today, was first formulized by Dorigo et al. in Reference Dorigo, Di Caro and Gambardella1999. The recent book by Dorigo and Stützle (Reference Dorigo and Stützle2004) gives a more comprehensive description about ACO. The definition of the ACO metaheuristic covers most (if not all) existing ACO variants for discrete optimization problems.
We use an altered version of ant colony system, one of the popular ACO's variant, which is originally proposed by Dorigo and Gambardella in Reference Dorigo and Gambardella1997. Then, it is combined with a ternary local search procedure to improve the quality of solutions. To prevent falling into a local optimum, we use the routine of the SA algorithm as a framework for the proposed SA-ACO. Figure 2 demonstrates the flowchart of the general structure of the presented algorithm. The parameters we employed are defined as follows:
-
T o is the initial temperature,
Fig. 2.General structure of the proposed algorithm.
-
T final is the final temperature,
-
T is the current temperature,
-
objcur is the objective function value of the current solution,
-
objnew is the objective function value of the new solution;
-
B is a constant,
-
CS is the cooling schedule, and
-
max_iter is the number of iterations of local search procedure in each temperature.
4.1. Initialization
In the proposed algorithm, u ants must generate all initial solutions. All ants are located in vertices at random. Then in each step, the vth ant uses transition rule for selection of the next vertex to construct the route. Therefore, the vth ant inserted in vertex i moves toward vertex j through hubs k and m, which is determined based on the following transition rule:
$$(\,j,k,m) = \left\{ {\matrix{ {\arg \mathop {{\rm \; max}}\limits_{(l,s,t) \in N_i^v} \; \lcub {{\lsqb {{\rm \tau}_{ilst}} \rsqb }^{\rm \alpha} {\lsqb {{\rm \eta}_{ilst}} \rsqb }^{\rm \beta}} \rcub \quad \hskip -1.2pt {\rm if}\; q \le q_0} \hfill \cr {\lpar {J \comma\, K \comma\, M} \rpar \quad\quad\quad\quad\quad\quad\quad\quad \hskip 2.7pt {\rm if}\;q \gt q_0} \hfill \cr}.} \right. $$
In the above equation,
${\rm \eta }_{ilst} = 1/(\tilde c_{ilst} + f_s + f_t)$
is a heuristic value; τ
ilst
is the pheromone value from vertex i to l through hubs s and t; α and β are parameters controlling the relative importance of global and local preferences for edge (i, j);
$N_i^v $
is the set of possible neighbors of the vth ant at vertex i, which means the collection of vertices at the neighborhood of vertex i that the vth ant have not met yet; q is a random number uniformly distributed in [0, 1]; and q
0 is a parameter (0 ≤ q
0 ≤ 1) that determines the relative importance of the current situation and the neighborhood search situation, which means the ant moves toward the best direction with the probability of q
0. The random variable (J, K, M) is selected according to the probability distribution given below:
$$P_{ijkm}^v = \displaystyle{{{\lsqb {{\rm \tau}_{ijkm}} \rsqb }^{\rm \alpha} {\lsqb {{\rm \eta}_{ijkm}} \rsqb }^{\rm \beta}} \over {\mathop \sum \nolimits_{(l\hbox{,}s\hbox{,}t) \in N_i^v} \; {\lsqb {{\rm \tau}_{ilst}} \rsqb }^{\rm \alpha} {\lsqb {{\rm \eta}_{ilst}} \rsqb }^{\rm \beta}}} \quad {\rm if} \; (\, j \hbox{,}\, k \hbox{,}\, m) \in \; N_i^v. $$
In the above equation,
$P_{ijkm}^v $
is the probability of moving from vertex i to j for the vth ant through hubs k and m, respectively. This amount is equal to zero for all vertices that do not belong to the neighborhood of i.
Note that we use (i, j, j, i) to show the direct service from spoke vertex i to j. By using this notation, we can unify our formulations for direct service and indirect service (through intermediate hubs) in the proposed algorithm. Moreover, the ant's moves are selecting and opening hubs along with the route construction. It implies that the ant needs to check the capacity of hubs and edges during this process. In addition, by opening a hub for a specific cluster, other paths that are connected to this cluster have to go through that hub. These issues are resolved by updating
$N_i^v $
based on the thus-far opened hubs and the remaining capacities of hubs and edges.
4.2. Pheromone updating rule
We use the two following forms to update pheromone values:
-
General update rule: After all ants complete their routes, to increase the pheromone of the edge with the highest performance, the following rule is used:
(13)where ρ is the global evaporation rate of the pheromone intensity and is chosen between 0 and 1, T bs is the best-so-far route, and
$$\eqalign{{\rm \tau} _{ijkm} \leftarrow \lpar {1 - {\rm \rho}} \rpar \times {\rm \tau} _{ijkm} + {\rm \rho} \times \Delta {\rm \tau} _{ijkm}^{{\rm bs}} \quad \forall \lpar {i \hbox{,}\; j \hbox{,}\, k \hbox{,} \,m} \rpar \in T^{{\rm bs}}\comma}$$
$\Delta {\rm \tau} _{ijkm}^{{\rm bs}} $
is the pheromone of the best-so-far route, which is calculated based on the following equation:
(14)where objbs is the cost of the best-so-far route (T bs).
$$\Delta {\rm \tau} _{ijkm}^{{\rm bs}} = \left\{ {\matrix{ {1/ \hbox{ob}\hbox{j}^{{\rm bs}}\quad \forall \lpar {i \hbox{,}\; j \hbox{,}\, k \hbox{,}\, m} \rpar \in T^{{\rm bs}}\;} \hfill \cr 0 \quad\quad\quad\quad\;\;\; {\rm otherwise} \hfill \cr}} \right. $$
Note that, as shown in Figure 2, the global pheromone update is only performed for the best-so-far solution in each temperature.
-
Local update rule: In order to prevent the solution from falling into a local optimum caused by the pheromone evaporating to zero, a local pheromone updating policy is deployed.
Therefore, when an ant passed the edge (i, j), its pheromone must be updated as follows:
(15)where ζ is the local evaporation rate of the pheromone intensity and is chosen between 0 and 1; τ0 is the initial pheromone value and is equal to
$${\rm \tau} _{ijkm} \leftarrow \lpar {1 - {\rm \zeta}} \rpar \times {\rm \tau} _{ijkm} + {\rm \zeta} \times {\rm \tau} _0 \comma$$
${1/n} \times {\rm ob}{\rm j}^{{\rm nn}}$
; n is the number of vertices; and objnn is the cost of the nearest neighbor route in which edges with minimum cost are selected between each origin–destination pair. Local update pheromone is caused when an ant uses edge (i, j); the amount of its pheromone is decreased and this leads other ants to search more unvisited edges by now. As a result, the search variety is increased and the probability of falling into a local optimum region is decreased.
To improve the quality of the solutions, we apply the local pheromone update based on a fitness function for edge (i, j). The fitness function evaluates the final value of the objective function by selecting edge (i, j). The local update is applied on the new edge if the value of the fitness function for this edge is better (lower in a minimization problem) than the fitness value of the current edge. Moreover, the local update may be applied on a new edge with a higher fitness value with a specific probability. Figure 3 depicts the proposed ACO.
Flowchart of the proposed ant colony algorithm.
Note that the ant is constructing the route (solution); therefore, we cannot directly use an objective function such as a fitness function. To overcome this issue, we use a prediction based on the relative effect of the selected edge on the final cost.
4.3. Local search
The local search in the proposed algorithm is a combination of the following three neighborhood searches:
Loc1: The hub of a cluster is reassigned to another vertex. This move is described as follows:
-
Step 1. Choose cluster r randomly and determine current hub h in this cluster.
-
Step 2. Choose randomly vertex v of cluster r as a new hub.
-
Step 3. If the new hub has enough capacity to serve all flows assigned to hub h, go to Step 4; otherwise, exit from this move and choose another move randomly.
-
Step 4. Reassign all paths to the new hub in cluster r.
Loc2: A direct service is changed to an indirect service through corresponding hubs. This move is described as follows:
-
Step 1. Choose direct service (i, j) randomly.
-
Step 2. Determine hub k, which is located in the cluster corresponding to vertex i.
-
Step 3. Determine hub m, which is located in the cluster corresponding to vertex j.
-
Step 4. If hub k, hub m, and corresponding edges have enough capacities to serve the flow from vertex i to vertex j, go to Step 5; otherwise, exit from this move and choose another move randomly.
-
Step 5. Replace direct service (i, j) with the indirect service through hub k and hub m.
-
Step 6. Update the remaining capacities of hub k, hub m, and corresponding edges.
Loc3: An indirect service through intermediate hubs is changed to a direct service. This move is described as follows:
-
Step 1. Choose indirect service (i, j, k, m) randomly.
-
Step 2. Replace indirect service (i, j, k, m) with direct service (i, j).
-
Step 3. Update the remaining capacities of hub k, hub m, and corresponding edges.
Note that the proposed algorithm randomly chooses these neighborhoods. This procedure iterates for a specific number in each temperature. Moreover, this part can be modified based on different assumptions included in the problem. The later feature makes a flexible algorithm for a variety of applications.
5. COMPUTATIONAL RESULTS
This section intends to evaluate the performance of the hybrid SA-ACO. We compare the result of the proposed algorithm with optimal solutions in small- and medium-sized instances. Then, the hybrid SA-ACO is compared with SA and the GA in large-sized instances. Finally, we compare the hybrid SA-ACO algorithm with one of the proposed heuristics in the literature. The platform of all experiments is a personal computer with a Pentium-4 2.7-GHz CPU and 2 GB RAM. The programs are coded in MatLAB. To develop case studies, we use a Turkish transportation data system that is very usual in the literature. Moreover, we employ a full factor experiment design on a variety of values to tune the parameters of the proposed algorithm based on similar studies and the problem features. Based on our experiments, because of the SA framework of the hybrid SA-ACO, the proposed algorithm is robust enough to be used with minor changes in parameters settings at least for instances applied in this paper. The parameters of all algorithms applied in this section are tuned in the ranges in Table 1.
The parameters of all algorithms tuned in ranges
5.1. Small- and medium-sized instances
In this section, we consider two classes of problem: small size and medium size. For each class, the algorithm is run 20 times, and the best results and also the average results are extracted. Finally, we obtain the optimal results using Lingo software and compare the optimal solutions and heuristic ones. All results are summarized in Tables 2 and 3. In each table the first column is the problem name, the second column is the number of vertices, the third column is the number of clusters, and the fourth column is the value of α. The CPU times are in seconds. As it is seen, the hybrid SA-ACO is not fast compared to Lingo in small-sized problems, and by increasing the size of instances, it is much faster than Lingo. In the case of 45 vertices, Lingo cannot reach the optimal solution, and therefore we report the best solution after 24 h run. For evaluating the gap between solutions we use the following formula:
$${\rm Gap}\left( {\rm \% } \right) = \displaystyle{{{\rm lingo\; objective} - {\rm best\; objective}} \over {{\rm best\; objective}}} \times 100.$$
The results of the hybrid SA-ACO and Lingo software in small sized problems
The results of the hybrid SA-ACO and Lingo software in medium sized problems
a The best result after a 24-h run.
Tables 2 and 3 demonstrate that the hybrid SA-ACO is highly competitive in medium-sized problems. Moreover, Figure 4 illustrates the performance of the hybrid SA-ACO in medium-sized problems.
Performance of hybrid simulated annealing–ant colony optimization in medium-sized problems.
5.2. Large-sized instances
To evaluate the performance of the hybrid SA-ACO in the large-sized instances, we compare it with SA and the GA. The SA algorithm that is used in this section is obtained by omitting the ACO part of the proposed algorithm in Figure 2. We replace the ACO with a random procedure to construct an initial solution. Then, this solution is improved by the ternary local searches in each temperature. Note that an initial solution is only produced in the beginning of the SA algorithm.
The GA used in this section is an altered version of prior studies (Holland, Reference Holland1992; Reeves, Reference Reeves, Glover and Kochenberger2003). The algorithm starts with a random population of individuals. Each individual has a chromosome that consists of two parts: hub indicator and direct service vector. The hub indicator is a genotype with integer values that indicates opened hubs in each cluster. Figure 5 shows the genotype of the hub indicator for an instance that has three clusters. The direct service vector is a vectorized version of the direct service matrix with binary values that indicate established direct services in each candidate solution (individual). Figure 6 shows the genotype of the direct service vector for an instance that consists of four vertices.
A hub indicator's genotype for a typical candidate solution.
A typical genotype of the direct service vector.
In each iteration of the algorithm, a new population is produced by applying genetic operators, crossover, and mutation, on current population. In crossover, two individuals (parents) are randomly selected and a mask's genotype of binary values is generated as well. Two children are produced from their parents based on the mask's genotype. If a mask's gen is 1, the first child inherits the corresponding gen from the first parent; otherwise, the first child inherits the corresponding gen from the second parent. The complement pattern of selected gens is used for generating the second child. The crossover operator is shown for typical parents in Figure 7.
A crossover operator on typical parents.
In mutation, a parent is selected randomly, and some of its gens are randomly changed in the solution domain. Figure 8 indicates this operator for a typical parent. After generating the new population, they are added to the current population. Then, all candidate solutions are evaluated by a fitness function. The fitness function is defined as the inverse of the objective function value for each individual. Finally, a roulette wheel is used to select a new generation of candidate solutions. This generation is employed as the current population for the next iteration of the algorithm. The algorithm terminates when a maximum number of generations is produced.
A mutation operator on a typical parent.
The results are given in Table 4. Figure 9 depicts the performance comparison of the hybrid SA-ACO with the above-mentioned SA and GA in large-sized instances. It can be seen that the proposed algorithm considerably outperforms the two other heuristics. Moreover, the difference between the hybrid SA-ACO and the SA, which use a same framework, highlights the role of the proposed ACO in the quality of the final solution.
Performance of the genetic algorithm, simulated annealing, and the hybrid simulated annealing–ant colony optimization in large-sized problems.
Comparison of genetic algorithm, simulated annealing and hybrid SA-ACO in large sized problems
Note: G1 (%) = (GA objective – SA objective)/SA objective × 100; G2 (%) = (SA objective – hybrid SA-ACO objective)/hybrid SA-ACO objective × 100.
5.3. Comparison with a proposed algorithm in the literature
We compare the hybrid SA-ACO algorithm with one of the proposed heuristics in the literature for a simplified version of our problem, in this section. We found that it is very hard to find some research that addresses a cluster HLP. Therefore, the random descent heuristic (RDH) proposed by Ernst and Krishnamoorthy (Reference Ernst and Krishnamoorthy1999) is selected due to its high performance and similarity to our problem. Figure 10 demonstrates the steps of RDH algorithm.
-
For A predetermined number
Fig. 10.Steps of the random descent heuristic algorithm.
-
Choose an initial solution randomly as the current solution
We use the AP data set presented in Ernst and Krishnamoorthy (Reference Ernst and Krishnamoorthy1999) to compare the hybrid SA-ACO with RDH. Table 5 reports the results for problems with 100 and 200 nodes based on the AP data set. As seen in this table, the proposed algorithm considerably outperforms the RDH algorithm adopted from the literature. This can result from the capability of the ACO algorithms for solving the graph-like problems in line with the main characteristic of the SA algorithms, which prevent trapping in local optima.
Comparison of RDH and hybrid SA-ACO in large AP problems
5.4. Convergence analysis
In this section, the convergence of the proposed algorithm is assessed. In this regard, we use an instance with 36 vertices, which includes five clusters and depicts the current status of the hybrid SA-ACO during the algorithm execution. Figure 11 shows the results. As it is seen in this figure, the procedure of the algorithm tends to be fixed in the final stage of the algorithm execution. In other words, the proposed algorithm shows a converging behavior in this diagram. This trend is very predictable because of using a SA framework in the main body of the hybrid SA-ACO. The proposed algorithm inherits the convergence among other important features of the SA algorithm. For further information about converging of the SA refer to Granville et al. (Reference Granville, Krivánek and Rasson1994) and Henderson et al. (Reference Henderson, Jacobson, Johnson, Glover and Kochenberger2003).
Variations during the proposed algorithm execution.
6. CONCLUSIONS
In this research the CMACHLP as an important branch of HLPs is considered for the first time. It is so applicable in many ranges of real-world problems, especially in telecommunication and cargo systems. The capacity constraints on the amount of the incoming and outgoing flow at the hub vertices and on the flow traversing the transportation lines are the major assumptions that are considered simultaneously to reflect the critical aspects of the real-world problems. The problem was formulated using an integer linear programming model. Then, based on the complexity of the problem that belongs to NP-hard problems, we applied a novel SA-ACO as a hybrid of the strongest metaheuristics to solve it. To enrich the SA-ACO, we applied three powerful neighborhood search strategies. The proposed method is performed extremely well, obtaining optimal solutions on nearly every run. These results were obtained with surprisingly little computational time, especially in medium-sized problems. Computational experiments demonstrate the hybrid SA-ACO performs highly competitively with respect to accuracy and efficiency of solution in comparison with the optimal solution, GA, and SA. Then, the hybrid SA-ACO is applied on the large-sized instances of the AP data set and considerably outperforms the RDH algorithm adopted from the literature. In addition, computational results show the robustness and converging behavior of the proposed algorithm as other advantages. For future research, multiobjective CMACHLP and a new efficient hybrid algorithm in that context can be interesting extensions.
Mohammad Mirabi is an Assistant Professor in the Group of Industrial Engineering at Ayatollah Haeri University of Meybod. He received his MS and PhD degrees in industrial engineering from Amirkabir University of Technology. Dr. Mirabi's research and teaching interests include project management, production management, flowshop scheduling, and economics.
Parya Seddighi received his BS in industrial engineering at Eyvanekey Institute of Higher Education.
