Optimization of a multi-modal tree hub location network with transportation energy consumption: A fuzzy approach

Abstract

A tree hub location problem (THLP) is a recently introduced extension of the classical hub location problem with an incomplete graph. The aim of this problem is to design a network more economic and practical. This paper presents a new multi-objective model to design a multi-modal tree hub location network under uncertainty. For this purpose, a fuzzy approach is applied to cope with the inherent uncertainty of input data in the THLP. One important issue, which is recently introduced in transportation network, is the amount of fuel consumption effected on economic and environmental problems, In this model, the amount of fuel consumption used in a transportation sector are accounted and the effect of road and vehicle types on consumed fuel in the THLP are studied. This model allows having different transportation modes between hubs and a set of capacity levels for each potential hub so that only one of them can be chosen. The objectives of this model include the minimization of energy consumption and the minimization of transportation costs and fixed costs of locating hubs and hub links. To solve the model, a multi-objective imperialist competitive algorithm (MOICA) is proposed to obtain the Pareto-optimal solutions of the problem. Furthermore, the performance of this algorithm is compared with non-dominated sorting genetic algorithm (NSGA-II).

Keywords

Multi-modal tree hub location problem energy consumption fuzzy approach meta-heuristics

1 Introduction

Hub location problems (HLPs) are used to design of transportation networks with several nodes that transfer passengers, cargo, express package and information between each pair of nodes. In such networks, the flows between origin-destinations are not shipped directly, but transferred through a set of nodes called hubs. Hubs are centralized facilities whose functions are to consolidate process and redistribute flows along the network. Indeed, a hub network has been operated with at least two major benefits that make it economic and practical. These benefits include decreasing the number of links between origin and destination nodes and use of discount factor for hub links because the volume of goods moving on the hub links is larger than the flows moving between non-hub and hub links.

In HLPs, the location of hub nodes and the allocation of non-hub nodes to these hubs should be determined. It is normally assumed that hubs are connected by a complete network that is to say that there is a link between each pair of hubs; however, some hub location problems have been presented the networks that do not have this restriction. In other words, two hubs are not necessarily connected by a direct link.

O’Kelly [1] formulated the HLP for the first time. According to the literature, HLPs are classified into several major types: (i) single/multiple allocation, (ii) capacitated/uncapacitated hubs, (iii) the existence or not of a pre-specified number of hubs to locate, (iv) complete/incomplete hub network, and (v) direct connection are allowed or not between spoke nodes (e.g., Alumur and Kara [2], Tan and Kara [3], Contreras et al. [4], Correia et al. [5] and Yaman [6]).

In a single allocation problem, each node or demand center can receive and send flows through exactly one hub node; whereas in a multiple allocation problem, there is not this restriction and every node can receive and send flows through more than one hub. The capacity of hubs can be limited or unlimited. In most real problems, the capacity of hubs is limited. In a p-hub problem the number of hub nodes is determined previously and in the case of a complete hub problem, every hub pair in the hub network is interconnected with a hub link; whereas on the incomplete one, there is not this restriction. The interested reader is referred to Alumur and Kara [2], Campbell et al. [7] and Zanjirani Farahani et al. [8] for surveys on HLPs.

In the real world, sometimes it is not applicable to establish a complete hub network or the cost of operating full interconnection is very high; in these cases, the application of an incomplete network arises and the connectivity must be achieved by the use of the minimum number of links. Generally, these networks are commonly used in telecommunication and transportation network.

In the case of the incomplete hub network, Campbell et al. [9, 10] proposed the hub arc location problems so that hub arcs are located with reduced unit costs. In these models, the fixed number of hub arcs is located with minimizing the total transportation costs. In this problem, hubs do not need to be connected unless forced. Alumur and Kara [11] presented the incomplete p-hub median, incomplete hub location with fixed costs, incomplete hub covering, and incomplete p-hub center network design problems. The aim of these models is to find the location of hub nodes, allocation of non-hub nodes to these hubs, and which hub links must be established between the hub nodes.

The tree of the hub location problem (THLP) is defined as a single allocation incomplete hub location problem where p hubs must be located and connected by means of a (non-directed) tree on a network.Contreras et al. [12] developed the tree hubs location problem with a single assignment. In this model, a fixed number of hubs have to be located and hubs are connected by means of a tree. The proposed problem combines the location network design and routing problems. In order to illustrate the impact of tree hub network, we consider a typical small example. Figure 1 depicts a network with 10 nodes, in which four of them are selected as hub nodes. Figure 1 (b) and (c) show the typical form of a complete and tree hub network, respectively. As it can be seen in Fig. 1, an incomplete hub network has a fewer links rather than a complete one.

Furthermore, the researches in the area of a multi-modal transportation network are increased during the last decades. It is usually assumed that there is one type of transportation mode in most of the hub location models; whereas in the real world, there are more choices for moving between origins and destinations. For example, many passengers transshipment occurs by using different transportation modes (e.g., air, ground, and water transportation systems (mainly air and ground transportation)). This assumption makes the network economic and operationally viable. There are several studies on the design of multi-modal transportation networks.

Alumur et al. [13] developed a multi-modal hub location network, in which different transportation modes between hubs and different types of service time between the origins and destinations are considered. Ishfaq and Sox [14] proposed an integrated hub-queue model and studied the effect of limited hub resources on the design of inter-modal logistics systems under service time constraints. Mohammadi et al. [15] also presented a new stochastic multi-modal hub location problem under uncertainty so that transportation time between each pair of nodes is uncertain and influenced by a risk factor. All the literature focuses on a multi-modal hub problem with a complete graph; however, the previous study shows the high quality of an incomplete hub network.

In this paper, we introduce a problem that addresses the multi-modal tree hub network under uncertainty. Moreover, one main concept which is not considered in THLPs is about the energy minimization and the amount of fuel used in transportation. The main contributions of this paper are as follows:

Designing a new multi-modal tree hub location problem under uncertainty;

Proposing multi-objective model including a new energy objective function to minimize the amount of consumed energy in the network and a cost objective function to minimize the transportation cost, fixed cost of locating hub and fixed cost of pathway construction between hub nodes;

Considering a set of the capacity level for each hub at a tree hub network, in which one of them can be chosen;

Describing a multi-objective imperialist competitive algorithm to solve the tree hub network;

Proposing a fuzzy approach to cope with the uncertainty and make the model more realistic;

This problem decides on the location and capacity of hubs, allocation of the non-hub node to these hubs, which hub links should be established between hubs and which type of transportation mode is chosen. The outline of this paper is as follows. In the next section, the literature review is described and in the third one, the mathematical formulation for multi-modal tree hub location is presented. In the fourth section, a multi-objective imperialist algorithm is described to solve the problem. The computational results are presented in the fifth section, followed by conclusions in the sixth section.

2 Literature review

Transportation sectors account for a large amount of energy consumption in many countries and have grown rapidly in recent years. It is also probable that this growth will continue in the next decades through the increase of a vehicle ownership level.Furthermore, transportation and the related amount of energy consumed in this sector are likely to grow up further with population growth, economic development, rapid industrialization, urbanization, agricultural development, rapid industrialization, and the like. The International Energy Agency (IEA) states that approximately 19% of global energy consumption and 23% of energy-related carbon dioxide (CO₂) emissions have been created by transportation sectors [16]. Hence, this section presents a brief overview of the existing literature on energy consumption in a transportation sector. Vanek and Morlok [17] presented a commodity-based approach to improve the energy efficiency of commonly used trucks and transfer more freight to energy efficient rail transportation. Tolouei and Titheridge [18] analyzed the effect of the vehicle mass on fuel consumption and secondary safety performance. The results presented that the higher vehicle mass increases the fuel consumption; furthermore, different combinations of fuel and transmission types affect the amount of consumed energy. They concluded that increasing the vehicle mass decreases the risk of driver injury in the event of a crash, so there should be a trade-off in the vehicle design between fuel economy and secondary safety performance imposed by mass. Besides their results, it should be noted that the trade-off in the vehicle design is often occurred between new technology cost, fuel consumption, and CO₂ emission and in most cases, volume capacity is paid more attention than the capacity of vehicle in term of transported weight. However, the mass criterion is less important in trade-off between technology cost and Co₂ emission. Liimatainen and Pöllänen [19] presented a new framework for modeling and assessing the trends of the energy efficiency in road freight transport. They estimated fuel consumption functions based on the vehicle and road types and predicted the trends of these factors that influence the energy efficiency until 2016.

Suzuki [20] proposed a new truck-routing approach to reduce fuel consumption and pollutant emission that minimizes the distance travelled by a delivery vehicle with a heavy payload. In this approach, heavy items deliver in the early segments of a tour while light items deliver in the later segments, such that the distance travelled by the vehicle with heavy items can be minimized. He concluded that this approach can produce fuel saving up to 6.9% over the existing methods. Mraihi et al. [21] used a decomposition analysis approach to determine the driving factors of transport energy consumption changes for the road mode. They found vehicle fuel and road vehicle intensities as principal factors andsuggested economic, fiscal and regulatory instruments used to lead road transport to an eco-friendly system by evolutiing the road transport related energy consumption and potential driving factors. Martínez and Miranda [22] studied on energy consumption and intensity of toll highway transport in Spanish roads. Wang et al. [23] surveyed on some reports included a number of theoretical and technical methods focusing on China’s transport energy consumption and saving. They also analyzed the trend of energy consumption in four current different transport sector involving road, railway, waterway and civil aviation in China. Gil et al. [24] proposed a general approach to reduce the overall energy consumption of urban rail. In this study, the main available strategies and technologies minimizing energy such as regenerative braking, traction losses reduction, comfort optimization functions, energy-efficient driving, smart power management and renewable energy micro-generation are provided. As a result, regenerative braking is illustrated as the greatest energy saving potential.

HLPs are the main branch of transportation networks, in which energy consumption and related costs need special attention. As it was mentioned, fuel consumption is a central variable in economic, environmental and engineering analyses in transportation and it also directly and indirectly determines the energy cost, driving behavior, emissions of greenhouse gases, and initial cost and performance, respectively. In recent years, energy consumption in the transport sector continues to rise, although this issue has been not seen in the HLP. Additionally, vehicle characteristics (e.g., frontal area, vehicle speed and vehicle mass and road feathers) should be considered in accounting fuel consumption. Since, there are different transportation modes in real life and these factors are different from one vehicle or road mode to other ones, the effect of each mode can be reflected in amount of fuel/energy consumption. Nevertheless, vehicle and road modes are not precisely considered in the previous studies on fuel consumption. So, there is a gap in the literature on the application of energy optimization-based models in HLPs, especially in tree hub location problem which is more applicable and economic. Indeed, most studies fail to accurately account for energy consumption factors of transportation networks. In this paper, we aim to fill this gap by taking these factors and their respective costs into account.

Furthermore, we face with uncertainty in the real world and the value of some parameters are not deterministic. On the other hand, there are many external factors affect the values of the parameters, which are not measurable; so, it is impossible to predict or detect the exact value of the all parameters. For example, the transportation cost depends on global economic, inflation rate, policy, etc. One approach for handling this uncertainty is a fuzzy approach. The uncertainty of the input data is also seen in a transportation network. In a fuzzy transportation area, Kumar and Kaur [25] proposed a new multi-modal transportation problem with a fuzzy coefficient in the objective function. The paper presented a fuzzy transportation problem so that the different cost parameters are represented by different type of LR flat fuzzy numbers and a novel method is identified for solving the same type of the fuzzy transportation problem. Vinotha et al. [26] considered a novel procedure for solving the total time minimization in a fuzzy transportation problem. They identified the transportation time, source and destination parameters as exponential fuzzy numbers, which are detected by the decision makers. Pramanik et al. [27] presented a new multi-objective solid transportation problem under a random fuzzy environment. In the paper, unit transportation costs, fixed-charge cost, demands, capacities, and resources of conveyances are illustrated as random fuzzy numbers. They also applied a reduced deterministic constraint in the proposed model and solved it using a generalized reduced gradient (GRG) method.

Bashiri et al. [28] presented a new p-hub center model and considered qualitative parameters as same as quantitative ones. The qualitative parameters (e.g., quality of service, zone traffic, environmental issues and capability for development in the future) can be considered in a transportation network. The paper also introduced a hybrid approach to the problem, in which the location of hubs is detected by both qualitative and quantitative parameters concurrently. Also, fuzzy VIKOR system is applied to cope with uncertainty. Ghodratnama et al. [29] proposed a new multi-objective uncapacitated hub location problem with fuzzy approach and solved the proposed model by five different meta-heuristic algorithms. The three objective functions of this model include of the total cost, congestion and hub installation costs. They assumed that three characteristics of transportation include costs, crowding and traffic costs, plus the costs of hub installation are not deterministic and defined the independent cost function in order to connect to the crowding rate and incurred cost in an exponential way.

3 Problem definition and formulation

In HLPs, there exists a network with n nodes includes the set of origins, destinations, and potential hub nodes. The p-hub median problem is one type of a hub location problem, which aims to locate p hubs in a network and allocate non-hub node to these hubs so that the total cost (time or distance) is minimized. The main assumptions of this paper are listed as follows:

Number of hubs is predefined and denoted by p;

Each non-hub node can be allocated to exactly one hub;

This network is an incomplete graph;

There are different transportation modes between hubs;

A number of capacity levels are considered for each hub.

3.1 Amount of fuel consumption

In this section, the amount of energy consumed in the hub network at arcs is described. When a vehicle traverses through an arc, a certain amount of energy is consumed. The amount of fuel/energy consumption depends on the travelled distances and varies from different vehicle modes. In spite of vehicles characteristics that have a significant effect on consumed energy, road type is also important because each road type has specific feathers (e.g., slope and friction) that has several effects on the amount of used energy. For example, a road with higher slope needs more energy to overcome the gravity insomuch the vehicle movement is occurred on the opposite direction of gravity.

In other aspect, driver behaviors (e.g., driving speed, intensity of braking and acceleration) have a great impact on energy efficiency. Since a journey includes periods of low-speed, urban driving and high speed, highway driving, in which highway driving has a large velocity ratio, speed cannot be constant and driver corresponds to control and constant speed driving. It is clear that high speed travelling consumes more fuel; however, speed changes throughout the journey has a more critical effect on amount of fuel consumption and rapid accelerator that is increased the amount of energy consumption.

In this paper, different vehicle types that have quite different fuel/energy consumption depending on their weight; power, its characteristics (e.g., friction of vehicle’s tires and air drag), road type and driving behavior are considered. It should be noted that the estimated energy consumption are for typical trips, which can be comprised of traveling on various roadway facility types (i.e. freeways, arterials, and local streets) and vehicle behavior effects are regarded as driving speed, intensity of braking and acceleration. The calculation of fuel/energy consumption (F) is complex and depends on several factors. A relevant study, which is considered here for energy consumption, is done by Simpson [30]. In this paper, the simplified calculation is as follows: $P_{road} = P_{aero} + P_{roll} + P_{accel} + P_{grade}$ (1) $P_{aero} = \frac{1}{2} ρ C_{d} {Av}^{3}$ (2) $P_{roll} = C_{r} M_{total} gv$ (3) $P_{accel} = k_{m} M_{total} av$ (4) $P_{grade} = M_{total} gZv$ (5)where,

P _road

Road load power

P _aero

Aerodynamic drag forces

P _roll

Rolling drag forces

P _accel

Power required for vehicle acceleration,

P _grade

Power for gradeability,

Vehicle speed (m/s),

Vehicle acceleration (m/s²),

Density of air (∼1.2 kg/m³),

C _d

Aerodynamic drag coefficient,

Frontal area (m²),

C _r

Rolling resistance coefficient,

M _total

Total vehicle mass (kg) (sum of empty and carried load),

Gravitational acceleration (9.81 m/s²),

Road gradient (%),

k _m

Factor to account for the rotational inertia of the powertrain

The road equation includes of four main components as illustrate in Equation 1. The first two components show the irrevocable power losses due to aerodynamic and rolling drag; whereas, the next two components represent the fully recoverable powers due to vehicle acceleration and hill climbing because kinetic and potential energy storage in the vehicle inertia and in theory are fully recoverable. These recoverable energies are actually recapture in the vehicle powertrain and should be attended as such. For example, the friction brake dissipates some of the kinetic and gravitational potential energy stored through the inertia of vehicle and calculates as follows:

$\begin{matrix} P_{braking} & = & (1 - k_{regen}) P_{inertia} \\ = & (1 - k_{regen}) k_{m} M_{total} av \end{matrix}$ (6)

P_inertia is the average of drivetrain power required for vehicle acceleration and k_regen is the regenerative braking fraction. Finally, the total amount of energy consumption is represented as:

$\begin{matrix} P_{total} & = & P_{road} + P_{braking} = \frac{1}{2} ρ C_{d} {Av}^{3} \\ + C_{r} M_{total} gv + k_{m} M_{total} av + M_{total} gZv \\ + (1 - k_{regen}) k_{m} M_{total} av \end{matrix}$ (7)

Some of above parameters are fixed at each arc (e.g., gravity) whereas some others are variable (e.g. load and speed). Some parameters (e.g., speed, load, acceleration and angle) change from one arc to another and some others (e.g., frontal area and speed) change from one mode to another. Let vehicle speed at an origin-destination arc (i, j) is denoted by v_ij(m/s) with distance d_ij (in meters) and road gradient Z = Z_ij. Since we suppose that there have been different types of modes between hubs, vehicle speed between hubs has been influenced by modes. So, speed between hubs (k, l) using mode m is depicted by $v_{ij}^{m}$ . In addition, M_total is divided into vehicle mass (M) depended on the mode type and load ( $f_{i j}^{'}$ ) carried at each arc. Consequently, the energy/fuel consumption on the arc can be calculated by: $\begin{matrix} P_{total} = \frac{1}{2} ρ C_{d}^{m} A^{m} (v_{ij}^{m})^{3} + C_{r}^{m} (M^{m} + f_{ij}) {gv}_{ij}^{m} \\ + k_{m}^{m} (M^{m} + f_{ij}) a_{ij}^{m} v_{ij}^{m} + (M^{m} + f_{ij}) {gZ}_{ij} v_{ij}^{m} \\ + (1 - k_{regen}^{m}) k_{m}^{m} (M^{m} + f_{ij}) a_{ij}^{m} v_{ij}^{m} \end{matrix}$ (8)

Consequently, the amount of energy consumed in a hub network is calculated by Equations (9) and (10). Equation (9) calculates energy consumption between hubs and Equation (10) calculates it through traversing from non-hub to hub nodes.

$\begin{matrix} P_{kl}^{m} (inter - hub) = \frac{1}{2} ρ C_{d}^{m} A^{m} (v_{kl}^{m})^{3} h_{kl}^{m} \\ + C_{r}^{m} [M^{m} h_{kl}^{m} + y_{klm}^{i}] {gv}_{kl}^{m} \\ + k_{m}^{m} [M^{m} h_{kl}^{m} + y_{klm}^{i}] a_{kl}^{m} v_{kl}^{m} \\ + [M^{m} h_{kl}^{m} + y_{klm}^{i}] {gZ}_{kl} v_{kl}^{m} \\ + (1 - k_{regen}^{m}) k_{m}^{m} [M^{m} h_{kl}^{m} + y_{klm}^{i}] a_{kl}^{m} v_{kl}^{m} \end{matrix}$ (9)

$\begin{matrix} P_{ik} (node - hub) = \frac{1}{2} ρ C_{d} A (v_{ik})^{3} x_{ik} \\ + C_{r} [M + (O_{i} + D_{i})] x_{ik} {gv}_{ik} \\ + k_{m} [M + (O_{i} + D_{i})] x_{ik} a_{ik} v_{ik} \\ + [M + (O_{i} + D_{i})] x_{ik} {gZ}_{ik} v_{ik} \\ + (1 - k_{regen}) k_{m} [M + (O_{i} + D_{i})] x_{ik} a_{ik} v_{ik} \end{matrix}$ (10)where O_i = ∑_jw_ij, D_i = ∑_jw_ji. In addition, w_ij and $y_{kl}^{i}$ show the flow to be sent from node i ∈ N to node j ∈ N and the flow that originates at node i and travels from hub k ∈ N to hub l ∈ N, respectively. x_ik defines an allocation variable whose value is 1 if node i allocated to hub k; otherwise, it is 0. Significantly, $h_{kl}^{m}$ shows that there is an arc between hubs k and l using mode m.

3.2 Mathematical formulation

In the proposed model, there is a set of n nodes that some of them are chosen as hubs and the remaining ones are defined as non-hub nodes, which are allocated to the located hubs. Hub nodes should be connected by means of a (non-directed) tree. In a tree hub network, the links that connect hubs are defined as a small tree and the links connecting non-hub nodes to the allocated hub nodes are defined as a large tree. Moreover, in a tree hub network with p hub nodes, there are p-1 hub links and the path between each pair of node is unique. It should be noted that the inter-hub transportation cost is different from other travelling costs; because in this network, it is possible to establish a various type of hub links and use different transportation modes that is made the hub arc transportation more efficient and economic. Thus, this parameter is depicted by its transportation mode ( ${\dot{c}}_{kl}^{m}$ ). In this model, uncertain parameters are also signed by “∼” (e.g., ${\tilde{c}}_{ik}$ ). Decision variables and parameters of the model are as follows.

Indices:

n = {1, … ,N} Set of nodes.

q = {1, … ,Q} Set of capacity levels.

m = {1, … ,M} Set of transportation modes.

Parameters: w _ij

Flow to be sent from node i to node j (i, j ∈ N).

c _ik

Unit transportation cost from node i to node k.

f_{k}^{q, m}

Fixed cost for installing a hub at node k using mode m at capacity level q.

{\dot{c}}_{kl}^{m}

Unit transportation cost from hub k to hub l using mode m.

F_{kl}^{m}

Fixed cost for establishing hub link between hub k and l using mode m.

{cap}_{k}^{q}

Capacity of hub installed at node k at capacity level q.

O _i

Total flow originating at node i ∈ N.

D _i

Total flow destined for node i ∈ N.

Number of hubs.

δ ^m

Inter-hub transportation discount factor, 0 δ^m< 1.

Variables: x _ik

1 if node i ∈ N is allocated to hub at node k ∈ N; 0, otherwise.

h_{kl}^{m}

1 if a hub link is established between hubs k ∈ N and l ∈ N using transportation mode m ∈ M; 0, otherwise.

z_{k}^{q, m}

1 if hub is located at node k ∈ N using mode m; 0, otherwise.

y_{kl}^{i}

Amount of flow from origin i that is routed through hubs k and l.

According to the above-mentioned notations, a new mathematical model for a multi-modal tree hub location problem is presented. The aim of this problem is to locate hub nodes, link them as a tree and allocate non-hub nodes to these located hubs. The first objective of this model minimizes the total amount of consumed energy at hub arcs and at arcs between non-hubs and hubs while the second one minimizes the total transportation cost. In other words, establishing a transportation mode with less consumed energy/fuel needs more cost in comparison with other modes. Respect to the above descriptions, the multi-objective problem can be developed as follows: $\begin{matrix} \min z = \sum_{k} \sum_{i} \frac{1}{2} ρ C_{d} A ({\tilde{v}}_{ik})^{3} x_{ik} \\ + \sum_{k} \sum_{i} C_{r} [M + (O_{i} + D_{i})] x_{ik} g {\tilde{v}}_{ik} \\ + \sum_{k} \sum_{i} k_{m} [M + (O_{i} + D_{i})] x_{ik} a_{ik} {\tilde{v}}_{ik} \\ + \sum_{k} \sum_{i} [M + (O_{i} + D_{i})] x_{ik} {gZ}_{ik} {\tilde{v}}_{ik} \\ + \sum_{k} \sum_{i} (1 - k_{regen}) k_{m} [M + (O_{i} + D_{i})] a_{ik} {\tilde{v}}_{ik} x_{ik} \\ + \sum_{k} \sum_{l} \sum_{m} \frac{1}{2} ρ C_{d}^{m} A^{m} ({\tilde{v}}_{kl}^{m})^{3} h_{kl}^{m} \end{matrix}$ $\begin{matrix} + \sum_{k} \sum_{l} \sum_{m} \sum_{i} C_{r}^{m} (M^{m} h_{kl}^{m} + y_{klm}^{i}) g {\tilde{v}}_{kl}^{m} \\ + \sum_{k} \sum_{l} \sum_{m} \sum_{i} k_{m}^{m} (M^{m} h_{kl}^{m} + y_{klm}^{i}) a_{kl}^{m} {\tilde{v}}_{kl}^{m} \\ + \sum_{k} \sum_{l} \sum_{m} \sum_{i} (M^{m} h_{kl}^{m} + y_{klm}^{i}) {gZ}_{ik} {\tilde{v}}_{kl}^{m} \\ + \sum_{k} \sum_{l} \sum_{m} \sum_{i} (1 - k_{regen}^{m}) k_{m}^{m} \\ (M^{m} h_{kl}^{m} + y_{klm}^{i}) a_{kl}^{m} {\tilde{v}}_{kl}^{m} \end{matrix}$ (11) $\begin{matrix} \min z = \sum_{k} \sum_{i} {\tilde{c}}_{ik} x_{ik} (O_{i} + D_{i}) \\ + \sum_{m} \sum_{i} \sum_{l} \sum_{k} δ^{m ˜ ̇} c_{kl}^{m} y_{klm}^{i} \\ + \sum_{k} \sum_{m} \sum_{q} {\tilde{f}}_{k}^{q, m} z_{k}^{q, m} \\ + \sum_{k} \sum_{l} \sum_{m} {\tilde{F}}_{kl}^{m} h_{kl}^{m} \end{matrix}$ (12) $\begin{matrix} s . t . \\ x_{ik} \leq \sum_{q, m} z_{k}^{q, m} \forall i, k \in N \end{matrix}$ (13) $\sum_{k} x_{ik} = 1 \forall i \in N$ (14) $\sum_{q} \sum_{k} \sum_{m} z_{k}^{q, m} = p$ (15) $\sum_{q} z_{k}^{q, m} \leq 1 \forall k \in N, \forall m \in M$ (16) $\sum_{m} z_{k}^{q, m} \leq 1 \forall k \in N, \forall q \in Q$ (17) $h_{kl}^{m} \leq \sum_{q} z_{l}^{q, m} \forall k, l \in N, l \geq k, \forall m \in M$ (18) $h_{kl}^{m} \leq \sum_{q} z_{k}^{q, m} \forall k, l \in N, l \geq k, m \in M$ (19)

$\begin{matrix} O_{i} x_{ik} + \sum_{m} \sum_{l} y_{lk}^{im} = \sum_{l} \sum_{m} y_{kl}^{im} \\ + \sum_{l} w_{il} x_{lk} \forall i, k \in N \end{matrix}$ (20)

$\begin{matrix} \sum_{i} (O_{i} + D_{i}) x_{ik} \leq \sum_{m} z_{k}^{q, m} \times c {\tilde{ap}}_{k}^{q} \\ \forall k \in N, q \in Q \end{matrix}$ (21) $y_{klm}^{i} + y_{lkm}^{i} \leq O_{i} h_{kl}^{m} \forall i, k, l, m \in N$ (22) $\sum_{m} \sum_{l} \sum_{k} h_{kl}^{m} = p - 1$ (23) $x_{ik} \leq 1 - \sum_{q, m} z_{i}^{q, m} \forall i, k \in N$ (24) $\begin{matrix} x_{ik}, z_{k}^{q, m}, h_{kl}^{m} \in {0, 1}, y_{klm}^{i} \geq 0, int \\ \forall i, k, l, m, q \end{matrix}$ (25)

In this model, the first objective function is directly related to minimization of energy consumption in the multi-modal tree hub network while the second one is related to the minimization of the total transportation cost of the THLP. The first two terms of the first objective function show the irrevocable power losses due to aerodynamic and rolling drag through the transshipment from non-hub to hub nodes. The next two terms depict the recoverable powers due to vehicle acceleration and hill climbing recaptured from non-hub nodes to hubs and the next one is a friction brake function that dissipates some of the kinetic and gravitational potential energy stored through the inertia of vehicle. The last fifth ones are as same as the first five terms; however, these functions compute the amount of consumed/recovered energy produced between hub transportation via different modes. In the second objective function, the sum of transportation cost, fixed cost of locating hub and establishing hub link are minimized. In the hub network, there is no any direct transshipment between non-hub nodes, and transportation between non-hub nodes should be occurred though passing hubs so each non-hub node should be allocated to the open hub for travelling to the another non-hub node. Constraints (13) ensure that non-hub nodes can only be allocated to the opened hubs. As it was mentioned, in the single allocation hub location, each non-hub node should be allocated to the exactly one hub. So, Equation (14) shows the single allocation for each non-hub node.

Constraint (15) requires that exactly p hubs are chosen. This paper assumes that each node has a potential to be selected as a hub node with different capacity levels and modes. The hub mode prepares the foundation of establishing a hub link via the selected mode. Therefore, when one node is selected as a hub, its capacity level and mode should be determined. Inequalities (16) guarantee that for each hub at most one capacity level is selected. Constraint (17) assures that each hub can choose at most one mode. It is clear that a hub link should be selected and established between hubs through designing the hub network. The link between hubs is set up via mode m when the related hubs are opened via mode m. Constrains (18) and (19) make sure that the establishing hub link using mode m occurs when the related hubs using mode m are opened. Equation (20) is the divergence equations for flow i at node k in a tree network. The left-hand-sides of equation add up the flow incoming to node k from node i, when node i is a non-hub node and it is allocated to hub k, plus the flow originated as node i and transferred from hub k to another hub l. The right-hand-side of the equation shows flow that origin in i and destination in any node allocated to the hub. Constraints (21) are the capacity constraints, which limit the amount of flows processed by each hub. As it was mentioned, the THLP requires the hubs to be connected by means of a tree. Constraints (22) and (23) show that variable h defines tree and p-1 hub links are chosen to define the tree. Constraint (24) ensures that only non-hub nodes can be allocate to the hubs. Finally, constraint (25) shows types of variables.

3.3 Fuzzy approach

In this real world, the value of some parameters, such as fixed cost of opening hubs ( $f_{k}^{q, m}$ ), fixed cost of establishing hub arcs ( $F_{kl}^{m}$ ), transportation costs between nodes (c_ik), speed of vehicle ( $v_{kl}^{m}$ ) and capacity of hub ( ${cap}_{k}^{q}$ ) are uncertain. The cost of transportation or opening hub/hub link is not deterministic because these costs are depended on the society, government policy, inflation and global economy. Moreover, other parameters include vehicle speed and hub capacity are also depended on some cases (e.g., road and weather conditions) and company policy that want to use a partial or total sector of its capacity (because of environmental or traffic issues), respectively. There are many different shapes of fuzzy numbers, in which a triangular fuzzy number (TFN) is the most popular one. Noticeably, the use of the TFN has three central advantages. Firstly, the TFN is the simplest method to implement a function with only three parameters that can be decreased to two parameters in the case of symmetric fuzzy numbers. Secondly, the unsophisticated realization of the elementary fuzzy arithmetical operations result to more use of this method. Finally, in this approach, the linearity of the model and the number of the objective functions are not changed. In this paper, a triangular fuzzy approach is also utilized to cope with the uncertainty.Jimenez et al. [31] presented a method to convert a fuzzy mixed-integer programming into an equivalent auxiliary crisp one. Some initial concepts of this method are briefly defined as follows.

Suppose that the fuzzy parameter is adapted to triangular fuzzy numbers and depicted by as $\tilde{ω} = {ω^{p}, ω^{m}, ω^{o}}$ , where ω^p, ω^m, ω^o are the pessimistic, most possible, and the optimistic value of $\tilde{ω}$ ; these values are obtained by decision makers. To show the performance of this method a simple fuzzy mathematical model is presented, in which all parameters of the model are considered as a fuzzy number: $\begin{matrix} min z = \tilde{Cx} \\ s . t . \\ \tilde{ax} \leq \tilde{b} \\ \tilde{Hx} \geq \tilde{r} \\ x \geq 0 \end{matrix}$ (26)

The equivalent crisp model is presented as follows [32]:

$\begin{matrix} min z = (\frac{C^{p} + 2 C^{m} + C^{o}}{4}) \times x \\ s . t . \\ (1 - α) (\frac{a^{p} + a^{m}}{2}) + α (\frac{a^{m} + a^{o}}{2}) x \\ \leq α (\frac{b^{p} + b^{m}}{2}) + (1 - α) (\frac{b^{m} + b^{o}}{2}) \\ α (\frac{H^{p} + H^{m}}{2}) + (1 - α) (\frac{H^{m} + H^{o}}{2}) x \\ \geq (1 - α) (\frac{r^{p} + r^{m}}{2}) + α (\frac{r^{m} + r^{o}}{2}) \\ x \geq 0 \end{matrix}$ (27)

Consequently, the equivalent auxiliary crisp model of our model is represented as follows: $\begin{matrix} \min \sum_{k} {\sum_{i} \frac{1}{2} ρ C_{d} A (\frac{v_{ik}^{p} + 2 v_{ik}^{m} + v_{ik}^{o}}{4})}^{3} x_{ik} \\ + \sum_{k} \sum_{i} C_{r} [M + (O_{i} + D_{i})] x_{ik} g \\ (\frac{v_{ik}^{p} + 2 v_{ik}^{m} + v_{ik}^{o}}{4}) \\ + \sum_{k} \sum_{i} k_{m} [M + (O_{i} + D_{i})] x_{ik} a_{ik} \\ (\frac{v_{ik}^{p} + 2 v_{ik}^{m} + v_{ik}^{o}}{4}) \\ + \sum_{k} \sum_{i} [M + (O_{i} + D_{i})] x_{ik} {gZ}_{ik} \\ (\frac{v_{ik}^{p} + 2 v_{ik}^{m} + v_{ik}^{o}}{4}) \\ + \sum_{k} \sum_{i} (1 - k_{regen}) k_{m} [M + (O_{i} + D_{i})] a_{ik} \\ (\frac{v_{ik}^{p} + 2 v_{ik}^{m} + v_{ik}^{o}}{4}) x_{ik} \\ + \sum_{k} \sum_{l} \sum_{m} \frac{1}{2} ρ C_{d}^{m} A^{m} \\ {(\frac{v_{klm}^{p} + 2 v_{klm}^{m} + v_{klm}^{o}}{4})}^{3} h_{kl}^{m} \\ + \sum_{k} \sum_{l} \sum_{m} \sum_{i} C_{r}^{m} (M^{m} h_{kl}^{m} + y_{klm}^{i}) g \\ (\frac{v_{klm}^{p} + 2 v_{klm}^{m} + v_{klm}^{o}}{4}) \\ + \sum_{k} \sum_{l} \sum_{m} \sum_{i} k_{m}^{m} (M^{m} h_{kl}^{m} + y_{klm}^{i}) a_{kl}^{m} \\ (\frac{v_{klm}^{p} + 2 v_{klm}^{m} + v_{klm}^{o}}{4}) \\ + \sum_{k} \sum_{l} \sum_{m} \sum_{i} (M^{m} h_{kl}^{m} + y_{klm}^{i}) {gZ}_{ik} \\ (\frac{v_{klm}^{p} + 2 v_{klm}^{m} + v_{klm}^{o}}{4}) \\ + \sum_{k} \sum_{l} \sum_{m} \sum_{i} (1 - k_{regen}^{m}) k_{m}^{m} (M^{m} h_{kl}^{m} + y_{klm}^{i}) a_{kl}^{m} \\ (\frac{v_{klm}^{p} + 2 v_{klm}^{m} + v_{klm}^{o}}{4}) \end{matrix}$ (28)

$\begin{matrix} \min \sum_{k} \sum_{i} (\frac{c_{ik}^{p} + 2 c_{ik}^{m} + c_{ik}^{o}}{4}) x_{ik} (O_{i} + D_{i}) \\ + \sum_{m} \sum_{i} \sum_{l} \sum_{k} δ^{m} (\frac{{\dot{c}}_{klm}^{p} + 2 {\dot{c}}_{klm}^{m} + {\dot{c}}_{klm}^{o}}{4}) y_{klm}^{i} & + \sum_{k} \sum_{m} \sum_{q} (\frac{f_{kqm}^{p} + 2 f_{kqm}^{m} + f_{kqm}^{o}}{4}) z_{k}^{q, m} \\ + \sum_{k} \sum_{l} \sum_{m} (\frac{F_{klm}^{p} + 2 F_{klm}^{m} + F_{klm}^{o}}{4}) h_{kl}^{m} \end{matrix}$ (29)

$\begin{matrix} s . t . \\ Constraints (13) - (20) and (22) - (25) \\ \sum_{i} (O_{i} + D_{i}) x_{ik} \leq \sum_{m} z_{k}^{q, m} \times a (\frac{{cap}_{kq}^{p} + {cap}_{kq}^{m}}{2}) \\ + (1 - a) (\frac{{cap}_{kq}^{m} + {cap}_{kq}^{m}}{2}) \end{matrix}$ (30)

4 Multi-objective imperialist competitive algorithm

The imperialist competitive algorithm (ICA) has been introduced as a new evolutionary algorithm for optimization problems. Most evolutionary algorithms (e.g., genetic algorithm and ant colony) simulate the natural processes, whereas the ICA simulates the social political process of imperialism and imperialistic competition [33]. In this paper, a multi-objective imperialist competitive algorithm (MOICA) is proposed to solve the model and obtained the Pareto optimal solutions. Finally, the performance of this algorithm is compared with a famous multi-objective meta-heuristic algorithm (NSGA-II).

4.1 Solution derived representation

Representing and decoding solutions in an easy and suitable way to searching the space is important. In the THLP, solution representation must be depicted the location of hub nodes and hub link and the allocation of non-hub nodes to hubs.

The presented continuous solution encoding is included the 1×(n+p) matrix [15]; where n and p are equal to the number of nodes and number of hubs, respectively. In the first step, bits 1 to n are filled with real random values and shown the allocation phase. In the next step, the biggest random value includes number 1; the next biggest one includes number 2, and so on. In this way, the real random numbers are replaced by a number of nodes. After that, bits (n+1) to the (n+p) are filled with random integer values limited 1 to n; however, it should be attended that the value of the (n+p)th bit should be equal to n. These numbers show the place of hub nodes in the previous part of the matrix. Finally, the sorted nodes in the first part of the matrix are allocated to the nearest right signed hub (Fig. 2).

In the THLP, p-1 hub links should be selected, in which p equals to the hub number. The proposed continuous solution encoding (CSE) of hub arc consists of a (p×p) matrix, which is filled with real randomvalues (elements on the main diagonal are equal to zero) and the two biggest ones are selected as the first hub arc which is connected hub nodes. In this matrix, these selected values are replaced by zero and in the second matrix, the selected hub arcs include the number equal to one and other bits are filled by zero. After that, numbers in each column of second matrix are summed; if this computed number is more than one, the rows of the related column that has a value more than zero must be found and related bits are filled by 2 indicated that there should not establish any link between these nodes. This method prevents to face with a sub-tour in the network. In the next step, we find the next biggest value as a candidate of a hub link and if this arc satisfies the previous condition, the arc is placed in the second matrix and other steps continues; else, the next arc candidates to enter in the matrix. Figure 3 represents the sample CSE with five hubs, a hub arc between hubs 1 and 3, and hubs 1 and 4 are selected as two first arcs. According to the proposed steps, we cannot have a hub link between hubs 3 and 4, and then a new hub arc selected. This procedure is shown in Fig. 3.

To depict the transportation mode, we use a n × n matrix where n equals to the number of nodes and fill this matrix with real random numbers. After selecting hubs, the transportation modes between these hubs are read in this matrix. Then, the value of bits are multiplied by the number of available modes and round up (Fig. 4).

4.2 Generating initial empires

The ICA is a robust method based on imperialism with the policy of extending the power. Each country has a social-political feature of a country (e.g., race, language, and economic feature). In this algorithm, the initial population is known as initial countries. This algorithm finds the best country with the best combination of socio-political characteristics.

Some of the best countries among the population are selected to be the imperialists and the rest of the population is divided among the mentioned imperialists as colonies [33]. The initial number of colonies of an empire is as follows: ${NC}_{n} = round {p_{n} . N_{col}}$ (31)where NC_n is the initial number of colonies of the n-th empire, N_col is the number of all colonies, and p_n is the normalized power of each imperialist defined by: $p_{n} = | \frac{C_{n}}{\sum_{i = 1}^{N_{imp}} C_{i}} |$ (32)where C_n is the cost of the n-th imperialist. Finally, the total power of an empire is calculated by: $\begin{matrix} TP of {Emp}_{n} = Total cos t ({imp}_{n}) \\ + ξ mean (total cos t (colonies of {emp}_{n})) \end{matrix}$ (33)where ξ is a position coefficient and is equal to 0.02 in our problem. As can be seen, the total power is mainly depended on the power of the imperialist country; however, the power of the colonies of an empire has small effect on the total power of that empire.

4.3 Searching the solution space

Three steps including assimilation, crossover and revolution are used to search the solution space for obtaining the better solution during each iteration. In the assimilation step, after dividing colonies between imperialists, the colonies move toward their corresponding imperialists (Fig. 5). For colony’s movement toward imperialist in the other direction, θ angle is applied with the uniform distribution between −γ and γ (where β is a number, greater than 1).

In the crossover step, 1-point, 2-point and uniform crossover operators are used for sharing colonies information with together. In the revolution step, three different procedures (i.e., swap, reversion and inversion operators) are applied for each colony. Finally, these created populations are merged and archive updating is done to select the best colonies for initial population. In archive adaption, the population are ranked and selected by the use of a non-dominance technique and crowding distance (which is equal to the archive size) and the best ones are selected as the imperialist as well as the remaining ones are considered as colonies. It is necessary to mention that the power of each country is mainly associated with its rank such that the weakest country in a higher rank is more powerful than the strongest in a lower rank and the countries with the same rank are compared with the crowding distancemetric.

4.4 Exchanging positions of the imperialist and a colony

In some decades, a colony may access to a better position than the relevant imperialist, and as a result the imperialist moves to the position of that colony and vice versa. Consequently, the total power of imperialists during the time is changed. Then the competition between imperialists start, the weakest empire loses its colonies and powerful ones try to grab it, and the empire that has lost all the colonies will collapse. Figure 6 shows the flowchart of the basic MOICA.

5 Experimental results

In this section, the performance of the proposed MOICA is compared with a well-known multi-objective algorithm called NSGA-II via four comparison metrics. In this paper, to optimize the performance of the algorithms, tuning of the parameters has been obtained throughout the response surface methodology (RSM). A regression equation analysis is applied to detect the response surface model. At the first step, parameters that statistically affect the algorithms are identified. Two different sizes including small (S) and large (L) sizes are applied to determine the parameters value. Additionally, each factor is calculated at two levels coded as –1 to identify the factor at its low level (L) and +1 to identify the factor at its high level (H). Finally, the coded variable computed by: $X_{i} = \frac{r_{i} - (\frac{h + l}{2})}{\frac{h - l}{2}}$ where X_i and r_i are shown as coded variable and real variable and h and l show the high level and low level quality, respectively.

The parameters of proposed MOICA for small (S) and large (L) sizes of the given problems have been shown in Table 1. (P_a, P_c and P_r show the probability of assimilation, crossover and revolution operators, respectively).

The NSGA-II assumptions are as follows:

The initial population is randomly generated.

Crossover operator is applied on random selected solutions using one of the operators included: 1-point crossover, 2-point crossover and uniform crossover, randomly.

Mutation operator is applied on random selected solution using one of the operators, which is included inversion, swap and reversion randomly.

Crossover and mutation ratios are set to 0.72 and 0.3, respectively.

The number of the initial population is set to 150 and 300 for small and large-sized problems.

The number of iteration as stopping criteria is set to 300 and 500 for small and large-sized problems.

5.1 Comparison metrics

The performance of the proposed MOICA is compared with NSGA-II via four comparison metrics described as follows (Mohammadi et al. 2013):

Quality metrics (QM): This index shows the contribution of each algorithm from Pareto solutions. This metric is measured by putting the non-dominated solutions together and calculating the ratios between non-dominated solutions of each algorithm. As a result, an algorithm with a higher value of the QM has a better performance.

Mean ideal distance (MID): This metric illustrates the distance between Pareto solutions and ideal point ( $f_{1}^{best}, f_{2}^{best}$ ):

$MID = \frac{\sum_{i = 1}^{n} \sqrt{(\frac{f_{1 i} - f_{1}^{best}}{f_{1, total}^{p, max} - f_{1, total}^{p, \min}})^{2} + (\frac{f_{2 i} - f_{2}^{best}}{f_{2, total}^{p, max} - f_{2, total}^{p, \min}})^{2}}}{n}$ (34)where n is a number of non-dominated solution and $f_{i, total}^{p, max}$ and $f_{i, total}^{p, min}$ are the maximum and minimum values of each functions. Consequently, an algorithm with a lower value of the MID has a better performance.

Diversification metric (DM): This index shows the spread of the Pareto solutions (Equation 35). An algorithm with a higher value has a better performance.

DM = \sqrt{(\frac{max f_{1 i} - min f_{1 i}}{f_{1, total}^{p, max} - f_{1, total}^{p, \min}})^{2} + (\frac{max f_{2 i} - min f_{2 i}}{f_{2, total}^{p, max} - f_{2, total}^{p, \min}})^{2}}

(35)

Spacing metric (SM): This metric describes the uniformity of the spread of non-dominated set solutions. An algorithm with a lower value has a better performance.

SM = \frac{\sum_{i = 1}^{n - 1} | \bar{d} - d_{i} |}{(n - 1) \bar{d}}

(36)where, d_i is the Euclidean distance between sequential solutions in the non-dominated set of solutions and

\bar{d}

is the average of these distances.

5.2 Data generation

Required data of the presented problem involves a number of nodes, number of hubs, transportation cost, fixed cost, capacity levels, speed, frontal area and flows. The value and distribution of parameters are generated randomly and illustrated as Table 2. It should be noted that random data are estimated based on the availability of sufficient data (e.g., available historical data) and relying on opinions, experiences and feelings of field experts. Some special numbers of hubs are shown for each number of nodes. Also, each problem instance is illustrated as “Number of nodes # Number of hubs” (e.g., 50#10 means 50 nodes and 10 hubs).

5.3 Computational results

In this section, some test problems are presented and solved by the proposed MOICA, whose performance is compared with NSGA-II. Tables 3 and 4 list the four comparison metrics for small and large-sized problems, respectively. A better performance of each algorithm in each metric is shown in bold.

As can be seen in Tables 3 and 4, proposed MOICA mostly has higher quality in comparison with the NSGA-II. In other words, MOICA has more contribution on the obtained Pareto optimal solutions. Furthermore, non-dominated solutions of the MOICA have less average values of the spacing metric. These data show that non-dominated solutions of the proposed MOICA are more uniformly distributed in comparison with NSGA-II.

In values of the diversification metric and mean ideal distance, in most of test problems, MOICA has a better application than NSGA-II. In a few test problems that the proposed MOICA has a less quality than the NSGA-II (e.g., 50#18, 100#5), it should be considered that the MOICA has a better performance in other metrics (e.g., spacing and diversity), meaning that this algorithm searches more accurately in space and it can be possible to obtain better solutions throughout the more number of iterations. The comparisons between runtime of these algorithms show that the NSGA-II has a higher speed in computations. As can been seen, node and hub numbers significantly affect computation time; additionally, the computation time is more sensitive to the hub number and in test problems with constant node number, the run time increases rapidly with arising hub numbers.

Figures 7 to 11 depict the performance of the proposed MOICA in comparison with NSGA-II regarding to obtain Pareto solutions for some numericalexamples.

5.4 Conflict of objective functions

The comparison results of the objective functions are depicted to show the conflict of functions. The conflict of energy consumption and transportation cost function is shown in Fig. 12. It can be justified by that establishing transportation modes with less energy/fuel consumption absolutely need more cost in comparison with other modes. For example, using an old technology needs more energy consumption and less fixed cost; however, the new technology needs less energy but it is expensive. So, if we want to consider environmental issues and use the vehicles that consumed less fuel, we should expend more money.

5.5 Sensitivity analysis

Figures 13 to 15 depict the sensitivity analysis on three main factors employed in the model (i.e., number of hub nodes, capacity level and the amount of flow), which effect on both objective function values (OFVs). Since the sensitivity analysis should be done based on the exact results till the effect of parameters is precisely illustrated, we analyze the effect of parameter changes on several small instances solved by the GAMS software, as a powerful optimization tool. Figure 13 illustrates the trend of OFVs by increasing the number of hub nodes. As can be seen, by increasing the number of hub nodes, OFV 1 is increased; however, OFV 2 is decreased and then increased. In other words, a trend of OFV 2 is reversed from hub number 3. Increasing the hub nodes leads to the use of the discount factor in transshipment at hub arcs and consequently decreasing the transportation cost. However, the cost of establishing hub nodes leads to have changes in the trend of OFV 2. Figures 14 and 15 show the changes in OFVs by an increase of the hub capacity level and the amount of the flow routed between nodes,respectively.

As can be seen, OFV 2 is more sensitive to the hub number (Fig. 13) and OFV 1 is more sensitive to the capacity level (Fig. 14); however, both OFVs follow the similar trend in Fig. 15.

Figures 16 and 17 are drawn to illustrate the effect of some main factors whose functions are only restricted to the related OFV. As illustrated in Fig. 16, when the amount of the considered parameters includes speed, vehicle mass, frontal area and acceleration increases, the first objective function values calculated consumed energy also increase. As mentioned above, these parameters do not have any impact on the cost function value. Finally, Fig. 17 depicts the effect of cost parameters on the cost objective function. It can be resulted that the objective functions are more sensitive to the speed parameter in OFV 1 and the fixed cost of the opening hub in OFV 2 when compared with the other relevant parameters.

6 Conclusion

A hub location network is one of the main groups of a transportation network, which is applied as switching, transshipment and sorting points in many-to-many transportation and distribution systems. This paper has presented a new multi-modal tree hub location problem (THLP) that minimizes the amount of energy consumption and transportation cost. Energy consumption has been divided into five different components including irrevocable power losses due to the aerodynamic and rolling drag, the fully recoverable powers due to vehicle acceleration and hill climbing (kinetic and potential energy) and the amount of energy recaptured in the vehicle powertrain. These items include some parameters depended on the type of road and vehicles. Since each transportation mode has different characteristics (e.g., cost, load, speed, accessibility and frontal area), the amount of consumed energy is different from one mode to another. Moreover, some other parameters applied in energy consumption formulation (e.g., angle and speed) depend on the type of roads. In this paper, the mode of transportation between hubs and mode of hub arcs are different. Thus, in this network, the energy consumption and transportation cost depend on the transportation mode so that transportation mode with less energy/fuel consumption absolutely needs more cost in comparison with another mode with higher energy/fuel consumption. This model has been considered these two conflicting objective functions simultaneously because the amount of consumed energy and costs especially in transportation sector have been become the main concern of most countries (especially developing countries) in recent years. Additionally, we have considered the several available capacity levels for each potential hub that make this network more flexible. This paper has proposed a fuzzy approach to handle the inherent uncertainty of input data in the given THLP. We have also applied a multi-objective imperialist competitive algorithm (MOICA) for solving this problem. Its performance has been compared with a well-known evolutionary algorithm, namely NSGA-II. The computational results have shown that the proposed MOICA has provided non-dominated solutions with higher quality and less average values of the spacing metric and in most test problems, the proposed MOICA considerably has a better performance than the NSGA-II. Finally, the sensitivity analysis on some factors employed in the proposed model includes a number of hub nodes, capacity level, the amount flow, some vehicle characteristics and costs are added to show the effect of them on both objective functions.

Footnotes

Acknowledgments

The authors would like to thank the Editor-in-Chief of the Journal of Intelligent and Fuzzy Systems and anonymous reviewers for their helpful comments and suggestions, which greatly improved the presentation of this paper. Additionally, this work has been supported financially by the Center for International Scientific Studies & Collaboration (CISSC) and the French Embassy in Tehran as well as the Partenariats Hubert Curien (PHC) program in France.

References

O’Kelly

M.E.

, The location of interacting hub facilities, Transportation Science 20(2) (1986), 92–105.

Alumur

and Kara

B.Y.

, Network hub location problems: The state of the art, European Journal of Operation Research 190 (2008), 1–21.

Tan

P.Z.

and Kara

B.Y.

, A hub covering model for cargo delivery systems, Networks 49(1) (2007), 28–39.

Contreras

, Díaz

and Fernández

, Lagrangean relaxation for the capacitated hub location problem with single assignment, OR Spectrum 31(3) (2009), 483–505.

Correia

, Nickel

and Saldanha-da-Gama

, The capacitated single-allocation hub location problem revisited: A note on a classical formulation, European Journal of Operation Research 207 (2010), 92–96.

Yaman

, The hierarchical hub median problem with single assignment, Transportation Research - Part B: Methodological 43(6) (2009), 643–658.

Campbell

J.F.

and O’Kelly

, Twenty-five years of hub location research, Transportation Science 46 (2012), 153–169.

Zanjirani Farahani

, Hekmatfar

, BolooriArabani

and Nikbakhsh

, Hub location problems: A review of models, classification, techniques and application, Computers & Industrial Engineering 64 (2013), 1096–1109.

Campbell

J.F

, Ernst

A.T.

and Krishnamoorthy

, Hub arc location problems: Part I introduction and results, Management Science 51(10) (2005), 1540–1555.

10.

Campbell

J.F.

, Ernst

A.T.

and Krishnamoorthy

, Hub arc location problems: Part II formulations and optimal algorithms, Management Science 51(10) (2005), 1556–1571.

11.

Alumur

and Kara

B.Y.

, The design of single allocation incomplete hub networks, Transportation Research -Part B: Methodological 43(10) (2009), 936–951.

12.

Contreras

, Fernandez

and Marín

, The tree of hubs location problem, European Journal of Operational Research 202 (2010), 390–400.

13.

Alumur

S.A.

, Kara

B.Y.

and Karasan

O.E.

, Multimodal hub location and hub network design, Omega 40(6) (2012), 927–939.

14.

Ishfaq

and Sox

C.R.

, Design of intermodal logistics networks with hub delays, European Journal of Operational Research 220(3) (2012), 629–641.

15.

Mohammadi

, Jolai

and Tavakkoli-Moghaddam

, Solving a new stochastic multimode p-hub covering location problem considering risk by a novel multi-objective algorithm, Applied Mathematical Modelling 37(24) (2013), 10053–10073.

16.

IEA, World Energy Outlook 2012. International Energy Agency, 2012, Paris, France.

17.

Vanek

F.M.

and Morlok

E.K.

, Improving the energy efficiency of freight in the United States through commodity-based analysis: Justification and implementation, Transportation Research Part D: Transport and Environment 5(1) (2000), 11–29.

18.

Tolouei

and Titheridge

, Vehicle mass as a determinant of fuel consumption and secondary safety performance, Transportation Research Part D: Transport and Environment 14(6) (2009), 385–399.

19.

Liimatainen

and Pöllänen

, Trends of energy efficiency in Finnish road freight transport 1995–2009 and forecast to 2016, Energy Policy 38(12) (2010), 7676–7686.

20.

Suzuki

, A new truck-routing approach for reducing fuel consumption and pollutants emission, Transportation Research Part D: Transport and Environment 16 (2011), 73–77.

21.

Mraihia

, Abdallaha

K.H.B.

and Abidb

, Road transport-related energy consumption: Analysis of driving factors in Tunisia, Energy Policy 62 (2013), 247–253.

22.

Pérez-Martíneza

P.J.

and Miranda

R.M.

, Energy consumption and intensity of toll highway transport in Spain, Transportation Research Part D, Transport and Environment 27 (2014), 1–5.

23.

Wang

Y.F.

, Lib

K.P.

, Xu

X.M.

and Zhang

Y.R.

, Transport energy consumption and saving in China, Renewable and Sustainable Energy Reviews 29 (2014), 641–655.

24.

González-Gil

, Palacin Batty

R.P.

and Powell

J.P.

, A systems approach to reduce urban rail energy consumption, Energy Conversion and Management 80 (2014), 509–524.

25.

Kumar

and Kaur

, Optimization for different types of transportation problems with fuzzy coefficients in the objective function, Journal of Intelligent & Fuzzy Systems 23(5) (2012), 237–248.

26.

Vinotha

J.M.

, Ritha

and Abraham

, Total time minimization of fuzzy transportation problem, Journal of Intelligent & Fuzzy Systems 23(2-3) (2012), 93–99.

27.

Pramanik

, Maiti

and Jana

D.K.

, A fixed charge multi-objective solid transportation problem in random fuzzy environment, Journal of Intelligent & Fuzzy Systems (2015), DOI: 10.3233/IFS-151542, Article in Press.

28.

Bashiri

, Mirzaei

and Randall

, Modeling fuzzy capacitated p-hub center problem and a genetic algorithm solution, Applied Mathematical Modelling 37(5) (2013), 3513–3525.

29.

Ghodratnama

, Tavakkoli-Moghaddam

, Kalami-Heris

S.M.

and Nagy

, Solving a novel multi-objective uncapacitated hub location problem by five meta-heuristics, Journal of Intelligent & Fuzzy Systems (2015), Article in Press, DOI: 10.3233/IFS-141525.

30.

Simpson

A.G.

, Parametric modelling of energy consumption in road vehicles. Ph.D. thesis, School of Information Technology and Electrical Engineering, The University of Queensland,Australia, 2005.

31.

Jimenez

, Arenas

, Bilbao

and Rodriguez

M.V.

, Linear programming with fuzzy parameters: An interactive method resolution, European Journal of Operation Research 177 (2007), 1599–1609.

32.

Pishvaee

M.S.

and Torabi

S.A.

, A possibilistic programming approach for closed-loop supply chain network design under uncertainty, Fuzzy Sets and Systems 161 (2010), 2668–2683.

33.

Atashpas-Gargari

and Lucas

, Imperialist competitive algorithm: An algorithm for optimization inspired by imperialist competitive, Proceeding IEEE Congress on Evolutionary Computation, (10.1109/CEC.2007.4425083), Singapore, (2007), pp. 4661–4667.