Solving a multi-objective chance-constrained hub covering location problem by discrete invasive weed optimization

This paper presents a stochastic bi-objective model for a single-allocation hub covering problem (HCP) with the variable capacity and uncertainty parameters. Locating hubs can influence the performance of hub and spoke networks, as a strategic decision. The presented model optimizes two objectives minimizing the total transportation cost and the maximum transportation time from an origin to a destination simultaneously. Then, due to the NP-hardness of the multi-objective chance-constrained HCP, the presented model is solved by a well-known metaheuristic algorithm, namely multi-objective invasive weed optimization. Additionally, the associated results are compared with a well-known multi-objective evolutionary algorithm, namely non-dominated sorting genetic algorithm. Furthermore, the computational results of the foregoing algorithms are reported in terms four well-known metrics, namely quality, spacing, diversification, and mean ideal distance. Finally, the conclusion is reported. Subjects: Evolutionary Computing; Logistics; Operations Research; Optimization; Stochastic Models & Processes; Supply Chain Management *Corresponding author: Reza TavakkoliMoghaddam, School of Industrial Engineering, College of Engineering, University of Tehran, Tehran, Iran E-mail: tavakoli@ut.ac.ir Reviewing editor: Zude Zhou, Wuhan University of Technology, China Additional information is available at the end of the article ABOUT THE AUTHORS Seyed Hossein Nikokalam-Mozafar obtained his BSc and MSc degrees in Industrial Engineering from Islamic Azad University and Tafresh University in Iran, respectively. His interest researches are in facility location, optimization and management. Behzad Ashjari is an assistant professor of Industrial Engineering at Tafresh University. He obtained his PhD degree from AmirKabir University of Technology in Iran. His research interests are in operations research, maintenance and supply chain management. Reza Tavakkoli-Moghaddam is a professor of Industrial Engineering at University of Tehran. He obtained his PhD in Industrial Engineering from Swinburne University of Technology in Australia. His research interests include facility layouts and location design, scheduling, and metaheuristics. Aida Omidvar received her BSc and MSc degrees in Applied Mathematics and Industrial Engineering from University of Isfahan and Islamic Azad University in Iran, respectively. Her interest researches are in facility location, optimization, metaheuristics and strategic planning. PUBLIC INTEREST STATEMENT Nowadays, transportation and communication networks are of a vital role in modern societies. Thus, optimization of the network design for managing better the traffics and decreasing transport time and cost in acceptable service quality is more important. In this paper, we optimize a network structure for achieving two objectives, namely minimizing the total transport cost and maximum transport time in each allocation when a demand is uncertain. Because of the complexity of the problem, we propose a meta-heuristic algorithm that hss inspired the behavior of invasive weed colonization. Received: 07 April 2014 Accepted: 13 November 2014 Published: 16 December 2014 © 2014 The Author(s). This open access article is distributed under a Creative Commons Attribution (CC-BY) 3.0 license.


PUBLIC INTEREST STATEMENT
Nowadays, transportation and communication networks are of a vital role in modern societies. Thus, optimization of the network design for managing better the traffics and decreasing transport time and cost in acceptable service quality is more important. In this paper, we optimize a network structure for achieving two objectives, namely minimizing the total transport cost and maximum transport time in each allocation when a demand is uncertain. Because of the complexity of the problem, we propose a meta-heuristic algorithm that hss inspired the behavior of invasive weed colonization.

Introduction
A hub location problem is one of the significant classes of facility location problems that after more than two decades from appearance nevertheless many papers of various journals are titled to hub networks design and applications. The background of these problems is to manage traffic better in the network that causes decrease in the transportation cost or time in an acceptable quality level of servicing. In hub and spoke networks, supply and demand points are usually connected together via special nodes, called hubs, which are equipped with some facilities. Therefore, the initiation flow from origin nodes routes destinations through hubs. Hub covering problems (HCPs) are one of the subproblems in a hub location area, which are applied in transportation systems and telecommunication networks. Campbell (1994) introduced three types of HCPs by defining three covering criteria depending on the part of flow route that is focused by the decision-maker for optimizing. Hence, according to the definition an origin-destination pair (i, j) is covered by hubs k and l if: (1) the cost from i to j via k and l (i → k → l → j) does not exceed a specified value; (2) the cost of each link of the path (i → k, k → l, l → j) does not exceed a specified value; and (3) the cost of origin-hub or hub-destination does not exceed a specified value. Also, he presented the first mixedinteger programming model for both single-and multiple-hub set covering problems, whose objective was to locate hub nodes with minimum fixed costs in the network satisfying all demand points. Kara and Tansel (2003) studied the single-allocation HCP and proposed various linear formulations. Wagner (2004) presented a model for single-and multiple-allocation HCPs. Afterward, Ernst, Jiang, and Krishnamoorthy (2005) proposed a different model for single HCPs using a new concept, radius.
Recently, most of the investigations in the optimization area are interested in studying the problems under uncertainty, and hub location problems are no exception in this matter. One of the approaches in decision-making under uncertainty is stochastic programming. Drezner and Scott (2013) presented a location-inventory model with a perishable product for a distribution center. They considered the transportation cost and time. Marianov and Serra (2003) proposed a stochastic model for hub location in an airline network by the queuing theory and chance-constraint approach. Sim, Lowe, and Thomas (2009) introduced a p-hub center problem and used a chance-constrained formulation for minimizing the service-level requirement. Lium, Crainic, and Wallace (2009) addressed the stochastic demand in the service network design. Contrerase, Cordeau, and Laporate (2011) studied stochastic uncapacitated hub location problems with uncertainty in demand, dependent, and independent transportation costs and described a Mont-Carlo simulation to solve the problem. Alumur, Nickel, and Saldanha-da-Gama (2012) addressed several aspects concerning hub location problems under uncertainty, such as setup cost and demand. Most of the optimization problems involve optimizing more than one objective function simultaneously. Hence, the decision-maker usually deals with trade-offs between some conflictive objectives to make an optimal decision. It has been another challenge to researchers on study of the problem with realistic assumptions. There are some researches about multi-objective hub location problems (Costa, Captivo, & Clímaco, 2008). Mohammadi, Tavakkoli-Moghaddam, and Rostami (2011) proposed a new weight-based multiobjective imperialist competitive algorithm (MOICA) for the HCP. The total transportation cost and the time service in the hub nodes as two objective functions were minimized. Furthermore, Mohammadi, Jolai, and Tavakkoli-Moghaddam (2013) presented a new stochastic multi-mode transportation model for a hub covering location problem under uncertainty. They minimized two objectives, the total current investment costs and the maximum transportation time between each origin-destination pair in the network simultaneously. In addition, they proposed a MOICA that compared two wellknown multi-objective evolutionary algorithms. Mohammadi, Torabi, and Tavakkoli-Moghaddam (2014) considered a new hub location problem (HLP) considering environmental aspects for air and noise pollution of vehicles. They proposed a mixed possibilistic-stochastic programming method along with two meta-heuristics, namely simulated annealing and imperialist competitive algorithm. Ghodratnama, Tavakkoli-Moghaddam, and Azaron (2013) presented a new mathematical model for a hub location with heterogeneous fleet vehicles and proposed three meta-heuristics; namely, genetic algorithm, particle swarm optimization, and simulated annealing to solve the model.  proposed a novel fuzzy bi-objective HCP for minimizing the total cost of transportation, consisting of transportation, covering, opening and reopening the facilities in hubs, activating the facilities of hubs and transporter purchasing, and the sum of times of shipping commodities by transporters from origin to destination. Table 1 is a short extract of other literature studies in HCPs (Calik, Alumar, Kara, & Karasan, 2009;Hamacher & Meyer, 2006;Karimi & Bashiri, 2011;Lowe & Sim, 2012;Mohammadi et al., 2011;Qu & Weng, 2009;Tan & Kara, 2007;Weng & Wang, 2008. Nowadays, transportation and communication networks are a vital role in modern societies. Thus, optimization of the network design for managing better the traffics and decreasing transport time and cost in acceptable service quality is more important. In this paper, we optimize a network structure for achieving two objectives, namely minimizing the total transport cost and maximum transport time in each allocation when a demand is uncertain. Because of the complexity of the problem, we propose a meta-heuristic algorithm that has inspired the behavior of invasive weed colonization. This paper surveys multi-objective single-allocation HCPs with the uncertain demand and transportation time using chance-constrained programming (CCP). CCP is one of the powerful approaches in modeling stochastic systems by replacing the deterministic constraints involving the effects of uncertainty with probabilistic ones. These new constraints hold at least the reliability level. The authors assume that the demand and transportation time between the nodes are the only uncertain parameters. The target of the problem is to find the location of hubs and allocate the non-hub nodes to hubs such that the transportation cost between each origin-hub or hub-destination links is within a considered economical bound. Moreover, the maximum transportation time between origin-destination is minimized simultaneously. The variable capacity of hub nodes is one of the different aspects of this problem with the other hub location models in the literature. The authors propose a new formulation for the problem and solve the multi-objective stochastic HCP using a novel artificial intelligent technique called as discrete invasive weed optimization (IWO) that is more efficient in comparison to the older famous The rest of this paper is organized as follows. In the next section, we first describe the deterministic model and then employ a CCP approach to model a stochastic environment. Some definitions and concepts of multi-objective optimization are presented in Section 3. Section 4 consists of an overview on the proposed multi-objective invasive weed optimization (MOIWO) algorithm and its structure. Subsequently, we present our proposed solving method. All the computational results and conclusions are brought in Sections 5 and 6, respectively.

Problem definition
This section intends to define the studied hub location problem and its assumptions. At first, the notations are introduced as follow. Coverage radius for the hub node k χ Collection cost unit flow per unit of distance α Discount factor between the hubs δ Distribution cost unit flow per unit of distance x ik 1, if node i is allocated to a hub at node k; 0, otherwise (decision variable) y i kl Amount of flow originated at node i is routed via hub nodes (decision variable) z k Variable capacity of hub node k (decision variable) Let a complete network with n nodes, some of them, have potential to being a hub. All nodes are connected together via hub nodes. Every node is only assigned one hub; usually this type of problem is known as a single-allocation hub location problem in the literature. In this model, the transportation cost consists of collection (i.e. from origin node to hub, i → k), transfer (i.e. between the hubs, k → l), and distribution (i.e. from hub to destination, l → j), which are represented by the coefficients χ, α, and δ, respectively. Our deterministic model is based on one of the effective three-index flow variables for hub location problems proposed by Ernst and Krishnamoorthy (1999) and improved by Correia, Nickel, and da-Gama (2010). They interpreted the traffic flow originating from a particular node i as a commodity by defining the variables y i kl . Often in transportation networks, opening costs of hub nodes are dependent on the amount of traffic flow in them so that the cost of equipping hub nodes with facilities increases by rise in the amount of flows. On the other hand, the hub capacity affects the import flow in hub nodes directly and, like most of researches in the literature of capacitated hub location problems, only adding a capacity constraint to the model can be restricted the minimum and maximum of input/ output flow in a hub. Hence, in this paper, we consider variable capacities for hubs and moreover define a new parameter, g k , that represents the cost for each unit of capacity development of hub k. By this assumption, the proposed model is different with the other ones in hub location literature.

Mathematical formulation
Following this, the authors firstly elaborate the formulation of the studied problem in deterministic conditions. According to the past explanations, the deterministic model is as follows.
It is a quadratic formulation with three decision variables for finding location of hub nodes, allocations, and the hub capacities represented by x ik ,y i kl and z k , respectively. The objective function (1) calculates the total transportation costs of the hub and spoke network. The objective function (2) and Constraint (8) ensure to minimize the maximum transportation time between the nodes for every origin-destination pair (i, j). Constraint (3) ensures each node is a hub or allocated to a hub node. Constraint (4) represents that non-hub nodes are only assigned to one hub. Constraint (5) indicates for every origin-destination pair (i, j), if node j is allocated to hub l then the maximum amount of flow between i and j is equal y i kl . Constraint (6) ensures that node i can only be allocated to k, if cost c ik is, at most, the radius r k . Constraint (7) is the capacity constraint for each hub node that represents the amount of hub traffics. Constraints (9) and (10) show the domain of decision variables. The main difference in this formulation and the proposed one (Ernst & Krishnamoorthy, 1999) are in Constraint (5)

Stochastic multi-objective single-allocation HCP
In the real world, some parameters of the models are uncertain. Therefore, stochastic modeling approaches may be more efficient for studying the natural environments. CCP is one of the powerful approaches in modeling stochastic systems, in which the variation of the uncertain parameter is based on the known distribution. Charnes and Cooper (1959) introduced CCP, in which the main idea of their proposed method is to relax the constraints in deterministic programming and replace them with probabilistic ones. CCP holds the constraints with at least reliability level ( ≺ 1). Indeed, an uncertain parameter has variations in an interval and the expected value of it is considered for deterministic programming. In the other words, by changing the parameter value, the solution space is different from deterministic programming. When = 1, then the probabilistic constraints are equal deterministic ones. W ij is assumed to be normally distributed with a mean (w ij ) and a standard deviation ( ij ). With this assumption, the total flow originated at node i (O i ) and destined to node i (D i ) are also uncertainty parameters with normal probability distributions.
) show the mean and standard deviation of it, respectively, that can be calculated simply as follows: In addition, T ij is assumed an uncertain parameter with a normal distribution.
According to the explanations, the stochastic model obtained by replacing Constraints (5), (7), and (8) in the deterministic model with three following probabilistic constraints are given below: Clearly, according to the statistics theories, Constraints (14-16) can be simplified to Constraints (17-19), respectively.
So, the stochastic model can be written as follows: Min Min E( )

Multi-objective optimization
Considering more than one objective function to optimize a problem leads to conceptual changes in meaning of the term "optimize" to "Pareto-optimal set." Without loss of generality, let an optimization problem with n objective functions as follows (Coello Coello, Lamont, & van Veldhuizan, 2007): Convergence to the Pareto-optimality set and maintenance of diversity in its solutions has two general goals, in which multi-objective optimization will be achieved.
All non-dominated solutions are classified in a front, which is called the first front (front1). Afterward, the remaining solutions are categorized again based on non-domination as a new front (front2) and so on. Figure 1 depicts three non-dominated sorting fronts. To estimate the density of solutions' another attribute, the crowding distance is calculated by the following formula (22) for solutions in the same front (Deb, Pratap, Agarwal, & Meyarivan, 2000). The minimum and maximum objective function values are assigned an infinite distance value.

Discrete IWO
Recently, many researchers in optimization areas tend to solve problems under uncertainty and dynamic information that cause high complexity and difficulty. However, often classical approaches cannot be more efficient in a dynamic environment. Here, alternative techniques (e.g. meta-heuristics) solve these problems effectively in acceptable time and quality by patterning out the behavior of nature to conform events. The IWO algorithm is introduced by Mehrabian and Lucas (2006) as an efficient evolutionary algorithm, especially in nonlinearity and non-convexity. They evaluated the Figure 1. Illustration of nondominated sorting fronts for a minimization problem.

Front2
Front3 capability of the algorithm in finding global extermum of optimization test functions and compared with genetic, memetic, shuffled frog leaping, particle swarm optimization, simplex, and direct search-simulated annealing algorithms. Until now, various studies used this novel algorithm for solving their problems, such as array antenna synthesis problems (Karimkashi & Kishk, 2010), timecost trade-off problems (Ramezani-Ghalenoei & Hajimirsadeghi, 2009), no wait-two-stage multiprocessor flow shop scheduling problems (Moradi-Nasab, Shafaei, Rabiee, & Mazinani, 2012) and multi-objective problems (Kundu, Suresh, Ghosh, Das, & Panigrahi, 2011). All the results show that IWO outperforms most of the oldest powerful meta-heuristics. It is a population-based algorithm inspired by natural behavior of colonizing weeds. The solving process of an optimization problem by IWO can be summarized in four steps fundamentally below: Initialization: A finite number of feasible solutions as initial population of weeds is generated and being dispread uniformly randomly in the search space. Spatial dispersal: This step provides a selection mechanism for recognizing appropriate plants. Produced seeds must be distributed normally in every iteration with different decreasingly standard deviations (sd). This ensures generated seeds be around their plants. Equation 24 calculates the sd for each iteration (sd iter ).
where pow, sd initial , sd final , and iter max show the predefined nonlinear modulation index, initial and final standard deviations, and the maximum number of iterations, respectively (sd final ≺ sd initial ).
Competitive exclusion: Newly generated seeds in the previous step are added to the colony. When the maximum number of populations is obtained in the colony, the algorithm ranks all the weeds for screening better ones. Therefore, according to the defined number of colony population, the surplus is eliminated.
Here, we illustrate the proposed multi-objective stochastic algorithm based on discrete IWO for solving multi-objective single-allocation HCPs with CCP. The generality of this algorithm is similar to the main structure since we focus on differences of implementation. The algorithm is as follows: Step 1. Initialization: The initial population of weeds is generated and is dispread uniformly randomly in the search space. Each solution or weed is represented by two separated strings of lengths n and p, respectively, for coding a hub and spoke network with n nodes and p unknown hubs. The first string of length n is formed of different n random numbers into the set N, which denotes the allocations indirectly. The second string of length p consisting of p random natural numbers concerns the locations of hub nodes in the first string. Always, two assumptions must be observed for the second string: (1) the value of last array is equal to n and (2) the random generated natural numbers sorted in to ascending order. Hub nodes are marked in the new string according to the proposed locations of hubs in the hub location string. To achieve a non-hub nodes assignment to hubs, the unmark nodes are allocated to the right-hand side distinct array or hubs. Figure 2 shows the structure of a weed for coding a hub and spoke network with seven nodes and two hubs. The advantage of this representation is to decrease the probability of generating infeasible solutions. So, it is not necessary to check or repair weeds after each operation and finally influence computational time-saving.
Step 2. Reproduction and spatial dispersal: As described in steps b and c, in the previous section, all plants are evaluated and ranked according to fitness values. Therefore, the population must be ranked using non-dominated sorting as described in the previous sections. We note all constraints instead of covering, and the probabilistic ones are satisfied by representation mechanism. Covering restriction imposes on the solution space by employing penalty function. The difference in probabilistic constraints in comparison with deterministic ones is only a constant because of considering a reliability level and simulating an uncertain condition. Thus, it is clear that this subject leads to translation of the feasible space without any fundamental changes on it. So, the proposed formulation causes to simplify chance-constrained implementation. Flowering plants reproduce seeds and distribute them normally in the search space by Equations (23) and (24), respectively. The algorithm generates a random number (num) in the interval [sd final ,sd iter ]. Then, the num random positions of allocation string are selected and swapped. Figure 3 illustrates how the reproduction mechanism of an invasive weed works when the num is equal to 3. The cells containing nodes 1, 3, and 6 are randomly selected by the algorithm. The contents of the cells are exchanged one by one, specifically that of No. 1 is exchanged for that of No. 3 and then that of No. 3 for that of No. 5.
Step 3. Mutation: The idea of using the mutation operator for the IWO algorithm was proposed by Moradi-Nasab et al. (2012) that inspired herbicide-resistant weeds. In this step, the mutation is done on the percent of the previous step population, (p m ). The algorithm generates a random number. If this value is less than or equal to the predefined probability of mutation, the invasive weed is mutated. Generally, five types of mutation operators can occur in the algorithm.
Type one: The number of hubs (p) increases or decreases randomly. Figure 4(a) shows a mutated weed. The number of hubs is increased.
Type two: The percent of non-hub nodes reallocated where is a predefined parameter. The algorithm generates n random numbers in unit interval as a new string. The maximum of the generated numbers is found and assigned to the array, in which the node 1 is located in the main string. This rule is observed for the other assignments. percent of nodes are changed by new random numbers. Hub nodes locate according to the main string and non-hub nodes are reallocated by the  above rule. The algorithm generates a new string of random numbers in initial interval [0, 1]; for example [0.9, 0.2, 0.6, 0.87, 0.49, 0.7, 0.33]. In this new coding, the maximum of generated numbers (0.9) refers to node 1 and the maximum remaining numbers refers to node 2. Consequently, the invasive weed is coded as follows: [1, 2, 3, 5, 6, 7, 4] → [0.9, 0.87, 0.7, 0.49, 0.33, 0.2, 0.6]. The algorithm randomly selected 20% of the cells assigned to the non-hub nodes and replaces them with the new random numbers. In Figure 4(b), the cell No. 2 is selected and replaced with 0.41. The objective of this mutation is to change the structure of the spokes, so there is no change in the hub nodes. The mutated weed is decoded by algorithm in the above-mentioned manner.
Type three: Selecting two non-hub nodes in allocation string and reversing all the nodes between these positions (see Figure 4(c)).
Type four: Select a hub node randomly and change with a practicable non-hub node. What restrict this substitution are two assumptions (see Figure 4(d)).
Type five: Two non-hub nodes are selected randomly. The second node is located adjacent to the first non-hub node position. Therefore, the other nodes shifted right side (see Figure 4(e)).
Step 4. Competitive exclusion: In this step, the initial weeds, produced seed, and mutated weeds are merged together. Only pop size of this huge population with the best fitness can win this competition as new population for the next iteration. The algorithm employs this evaluation using non-dominated sorting.
Step 5. Go to Step 2 and repeat the algorithm until the stopping criteria are satisfied.

Computational experiments
In this section, the performance of the proposed MOIWO algorithm is tested for different instances. The authors use the Civil Aeronautics Board data-set and a random generation data-set for a network up to 100 nodes. For the random data-set, x and y coordinates are uniformly generated from U [0, 100000]. We assume that w ij and t ij take the values given in the data-sets. All the standard variations are considered as a coefficient ( ) of means. In addition, the studied problem is solved by the nondominated sorting genetic algorithm (NSGA-II). Deb et al. (2002) proposed an improved version of the NSGA-II. It ranks the populations due to the non-domination level.
Each problem is solved 10 times for both MOEAs. The number of population for the problems with up to 50 nodes is assumed less than larger ones. The NSGA-II and MOIWO parameter values for a given run are used for the problem provided in Table 2.
Good distribution of the proposed solutions, a wide-range convergence of non-dominated solutions and less distance of a non-dominated solution set to Pareto front are measured in multi-objective problems. The computational results of the two foregoing algorithms are reported in terms of the following four well-known metrics. The quality, spacing, and diversification metrics are used for the validation of the multi-objective evolutionary algorithm (MOEA) reliability.
Quality metric (QM): The ratio between non-dominated solutions found by algorithms is employed to measure the quality (see Figure 5).

Spacing metric (SM): It measures the uniformity of Pareto solutions' distribution as follows:
where d i denotes Euclidean distance between consecutive solutions, and the mean of d i is shown by d (see Figure 6).   Figure 7). Let (max f 1,i , min f 1,i ) and (f max 1,total , f min 1,total ) represent the maximum and minimum of the first objective values for each Pareto solution i and all of them, respectively. It is similar for another objective function. Therefore, the spread of proposed Pareto solutions can be measured by: Figure 8 shows the mean of distances between ideal solution and Pareto solutions obtained from MOIWO and NSGA-II. The mean ideal distance is calculated by the following formula (27). Less MID is another MOIWO preeminence (see Figure 7). Table 3 illustrates the results obtained by the NSGA-II and proposed MOIWO in terms of the foregoing four metrics.

Mean of distances between ideal solution and Pareto solution (MID):
The concluding results are as follows: (a) Higher quality of the proposed MOIWO algorithm for small, median, and especially large networks.    (c) Higher diversification of the Pareto solutions of the proposed MOIWO algorithm for different networks.
(d) Lower mean of distances between the ideal solution and Pareto solutions of the proposed MOIWO algorithm for different networks.
The Pareto solutions distributed more uniformly and the distance between ideal solutions and Pareto ones are lower with fewer variations.

Conclusion
In this paper, the authors have considered a hub location problem by considering realistic assumptions as uncertain parameters and two conflictive goals for optimization. A quadratic formulation has been proposed for multi-objective single-allocation HCP with the variable capacity in uncertainty demand and transportation time using CCP. Additionally, the authors have developed the MOIWO algorithm to solve the given problem. The reliability of the proposed MOIWO algorithm was evaluated by four well-known metrics; namely QM, SM, DM, and MID. The associated results are in evidence, in which the proposed MOIWO performance is appreciable, compared with the NSGA-II. These algorithms rank the individuals by non-dominated sorting; however, the proposed MOIWO algorithm can propose better Pareto solution sets regarding quality and diversification.