1 Introduction

Hub-and-spoke networks have wide applications in many areas, such as the airline industry, postal delivery, telecommunications, and public transportation networks. If hub nodes are fully-connected, the hub network is called a complete hub network. In an urban public transportation network with at least two modes of transportation [e.g., bus rapid transit (BRT), subway, freeways and monorail], hub nodes are commonly chosen from the nodes with less distance from other fast lines. Here, fast lines are assumed to be hub links, and regular bus lines are spoke links. Since direct links between all fast lines are too expensive and constructing them is not practically feasible, therefore the hub network is not complete in practice. In most hub location models in the literature, the demand is assumed to be known and independent of the condition of the hubs. However, especially for inessential goods and services, the demand is dependent on price, distance to the facilities, etc. (Redondo et al. 2012). Hence, considering a variable demand makes the model more realistic.

A brief overview of the literature (see Sect. 2) demonstrates that, in spite of the fact that public transportation systems can have a hub-and-spoke structure, this aspect has mostly been overlooked in the previous studies of hub networks. The models presented in this field have all been single-objective and negligent of the service levels and capacity constraints and involving just one transportation mode. Also, all these models were deterministic and the demand was assumed to be static, which is far from reality. In this study, the main focus is on designing a new hub network for public transportation, which fills the gaps identified in the previous models in this area and complies more with reality. Prime contributions of this study, which distinguish it from the published papers on this subject, can be summarized as follows:

  • Designing a multi-modal incomplete hub network and considering different types of hubs.

  • Considering three potential capacity levels for each hub node from which only one is chosen as the final hob’s capacity.

  • Considering a limited capacity for hub links.

  • Assuming the demand to be elastic and dependent on the location of the hubs.

  • Incorporating a fuzzy TH method (Torabi and Hassini 2008) to cope with uncertainty for solving small-sized problems.

  • Considering a multi-objective model that maximizes the total profit and minimizes the total time of transportation such that developing a transportation mode with shorter travel time results in an increase in the total profit.

  • Calculating the number of different kinds of vehicles used in hub links.

  • Solving a small size instances considering real data of the Qom monorail project and developing a meta-heuristic algorithm, namely MOPSO, to solve the large-sized instances of the given problem in a reasonable time using the Turkish data set.

The remaining of the paper is structured as follows. Sections 2 and 3 present a literature review and the mathematical model for the bi-objective incomplete hub median location problem, respectively. The solution method is presented in Sect. 4 in order to modify the fuzzy bi-objective model into a crisp one. Our solution algorithms are presented in Sect. 5, and the computational results are presented in Sect. 6 followed by conclusions in Sect. 7.

2 Literature review

O’Kelly (1986, 1987) introduced the first mathematical formulation for the hub-and-spoke problem as a quadratic integer programming model (Zanjirani Farahani et al. 2013). One main assumption of a typical hub location problem is that the hub network is complete; however, this assumption is unrealistic, especially in public transportation networks. Nonetheless, not many studies have been conducted on incomplete hub networks. Contreras et al. (2010) proposed a tree of the hub location problem with a single assignment (see also Sitek and Wikarek 2019). Martins de Sá et al. (2013) presented a hub location problem with a tree topology and solved it by a Benders’ decomposition approach that was an extension of the work of Papadakos (2008). Chiu et al. (1995), and Wang et al. (2006) presented a ring structure in hub network problems. In this topology, each hub node only links with two other hubs so there is only one loop in the hub network. Yoon and Current (2006) proposed a general hub network topology. In this topology, the form of the hub network is defined by the model, and any structure (e.g., a tree, ring or complete structure) can be created. O’Kelly et al. (2015) formulated a model to demonstrate how the fixed costs in transportation hub location models can generate different kinds of network types.

Although in real life cases there is always more than one objective to be optimized, the hub location literature mainly emphasizes single objective problems. However, for a realistic setting, service quality should be considered while maximizing the total profits (Campbell 2009; Yaman 2009; Alumur et al. 2012a, b). Correia et al. (2010a) developed a formulation for the capacitated single-allocation hub location problem and showed that the classical formulation is not complete in some instances. Correia et al. (2010b, 2011) extended a classical capacitated single allocation hub location problem in which the capacity of each hub is determined by solving the model. Azizi (2019) considered facility disruption in hub-and-spoke networks. In his proposed model, every demand node has a backup hub in the case of disruption. Corberán et al. (2016) also considered a single allocation hub location problem with link capacities. They proposed a meta-heuristic algorithm in accordance with the strategic oscillation and demonstrated that their proposed algorithm has a better performance compared to the tabu search. Ambrosino and Sciomachen (2016) presented a multi-modal capacitated hub location problem in freight logistics context and considered capacity bounds for both hub nodes and arcs.

In most hub location studies, the demand is presumed to be static and independent from the location of the hub facility whereas this is not a sensible presumption. In real life cases, the demand is a function of the utility perceived by the hub facilities, especially for inessential goods. Khosravi and Akbari Jokar (2014) considered elastic demand in the classical hub location problem. To the best of our knowledge, this study is the only one in the literature which assumes the demand dependent on the location of hub nodes.

Strategic planning may take considerable time to be implemented and may have a long-lasting effect. On the other hand, input data are not accurately known beforehand. Consequently, uncertainty should be taken into account in these kinds of problems. Contreras et al. (2011), Alumur et al. (2012a, b), Zhou and Liu (2007), Taghipourian et al. (2012), Mohammadi et al. (2014), Sadeghi et al. (2015) and Rahimi et al. (2016) considered a number of features of hub location problems under uncertainty.

For the first time, Nickel et al. (2001) presented new mathematical models of hub location problems (HLPs) in urban public transportation networks in 2001. They modified their models for public transportation systems by relaxing some classical assumptions of the hub location problem. Their proposed model was based on a general incomplete network. Afterwards, Gelareh and Nickel (2011) proposed a new formulation for un-capacitated multiple allocation hub location problems in urban transportation and liner shipping network design. Setak et al. (2013) proposed an incomplete hub location network for an urban transportation problem. In their model, non-hub nodes can be connected directly to each other and the network topology is incomplete. Karimi and Setak (2014) introduced a new formulation for a public transportation network with a hub structure. The advantage of their model in comparison to the previous models in this area is that the number of variables and constraints are reduced so that the computational time is reduced considerably. Sadeghi et al. (2018) considered a reliable p-hub covering location problem with links capacity and stochastic degradations in a transportation network. They proposed the artificial bee colony and differential evolution algorithms to solve their model. Bashiri et al. (2018) presented a new dynamic mobile p-hub location problem, in which hub facilities with a mobility feature are transferred to other nodes to meet demands. They proposed simulated annealing and genetic algorithms to solve this model.

On the basis of the above survey, it is clear that this paper gives answer to the several gaps existing in the literature by contributing along the following directions:

  • From the modelling viewpoint, we came up with a novel holistic bi-objective mathematical programming model for the public transportation network design encompassing several real-world limitations and constraints that have mostly been neglected in the literature filling most of the gaps present in the previous studies.

  • We have applied an efficient approach for solving the fuzzy bi-objective mixed-integer linear programming to cope with uncertainty (Torabi and Hassini 2008).

  • Concerning the methodological aspect, we have developed two efficient meta-heuristics based on NSGA-II and MOPSO algorithms that allow to address large-sized problems.

  • From the computational point of view, we solved a real-life problem and we performed an intensive sensitivity analyze phase in order to validate the effectiveness and efficiency of the proposed model and approaches.

3 Problem description and proposed model

In this section, first the problem statement is presented, and then the model assumptions are explained. Finally, the parameters, variables, and the bi-objective model are described.

3.1 Problem description

In this paper, a public transportation system with different types of transportation modes is considered. Each node can be selected as a hub. In each node, there is a choice to select different types of stations (i.e., BRT and metro stations) with different cost and attraction. As defined by Berman and Krass (2002), “The attractiveness of the facility is a function of the facility characteristics (e.g., size), characteristics of the facility’s location (e.g., availability of parking and proximity to major intersections), and other variables (e.g., marketing expenditures and prices) that may affect the perceived attractiveness of the facility”. Thus, a realistic assumption consists in considering the demand at each node as a function of the attractiveness of that hub station and a decaying function of the distance between that node and the hub station. As hub stations are not established in advance, the demand between origins and destinations depends on the location of these stations. In each hub node, one or both BRT and metro stations can be established; however, in each hub link, only one of the BRT or metro lines can be established. Building a BRT line is allowed only when there is a preexisting street between the two nodes; however, there is no such restriction for building a subway. Due to economic considerations and to the necessity of connecting different parts of the city, the subway network should be fully connected.

We assume here that each hub station has three potential capacity levels of which only one is selected as the final capacity of the station. The capacities of hub links are assumed to be limited too. The required number of transportation vehicles on each link is calculated by the model. The allocation of demand nodes to hub stations is assumed to be multiple. To avoid the allocation of each demand node to all the established hubs, an allocation cost is considered. The first objective of the model is to maximize the benefits of the total transportation, construction of hub stations, and deployment of the vehicles. The second objective is to minimize the total travelling time in the network. Since we consider the real-world situation, some of the input parameters are assumed to be uncertain.

3.2 Model formulation

This section presents a new mathematical formulation for the hub location problem in public transportation utilizing the modelling structure described in Sect. 3.1. The main assumptions considered while developing our mathematical model are as follows:

  • The number of hubs to be located is predefined.

  • Each non-hub node can be allocated to one or more hubs (i.e., multiple allocations).

  • The hub network is incomplete; however, it is a connected graph.

  • Hubs and arcs capacities are limited.

  • Each hub includes three capacity levels from which only one level is selected as the optimal one considering the total cost.

  • The capacities of transportation vehicles are limited.

  • All flows travelling from (to) a node pass through one or more hubs.

  • The demand of each node is a function of utility perceived by the customers at that node.

  • There are two types of transportation modes in the hub network (i.e., BRT and subway lines).

  • Building a BRT line is only allowed when there is a preexisting street between the two nodes; however, there is no such restriction to build a subway line.

  • The metro network should be connected; but there is no such restriction for the BRT network.

  • Some input parameters of the model are assumed to be uncertain.

It should be noted that the presented model is readily extendable to more than two transportation modes. A list of notations used in the models is presented below.

Indices and Sets

N:

Set of all origin and destination nodes

H ∈ N:

Set of all potential hub nodes

i, j ∈N:

Origin and destination indices

k, l ∈ H:

Hub nodes indices

M:

Set of the transportation modes

b, t, m ∈ M:

BRT and metro indices

Q:

Set of the capacity levels

q ∈ Q:

Capacity level index

Parameters

\( W_{ij}^{\rm max } \):

Maximum demand from i to j

\( p_{ik} \):

Income per unit of flow from origin i to hub k

\( p_{kl}^{m} \):

Income per unit of flow from hub k to hub l via mode m

\( c_{ik} \):

Cost per unit of flow from origin i to hub k

\( c_{kl}^{m} \):

Cost per unit of flow for hub link k–l via mode m

\( F_{k}^{mq} \):

Fixed cost for establishing a hub station at node k with transportation mode m and capacity level q

\( I_{kl}^{m} \):

Fixed cost for establishing a hub link between hub nodes k and l with transportation mode m

\( t_{ik} \):

Transportation time from node i to hub k

\( t_{kl}^{m} \):

Transportation time from hub k to hub l via mode m

\( t_{k}^{m} \):

Waiting time in hub k with transportation mode m

\( \alpha^{m} \):

Discount factor of transportation cost using mode m

\( \beta^{m} \):

Discount factor of transportation time using mode m

\( \varGamma_{k}^{mq} \):

Capacity of hub station k with capacity level q using transportation mode m

\( R_{kl}^{m} \):

Capacity of link k–l using transportation mode m

\( V^{m} \):

Capacity of vehicle m

\( S_{kl} \):

Binary parameter showing whether a route between hubs k and l exists

\( CV^{m} \):

Cost of vehicle m

\( J_{ik} \):

Cost of allocating non-hub node i to hub node k

P:

Number of hubs to be established

Decision variables

\( z_{ik} \):

Flow from node i to hub k

\( x_{ij}^{k} \):

Flow originated from node i flowing from hub k to node j

\( y_{ikl}^{m} \):

Flow originated from node i routed between hub nodes k and l with transportation mode m

\( a_{ik} \):

1 if non-hub node i is allocated to hub node k; 0, otherwise

\( h_{k}^{mq} \):

1 if hub k with transportation mode m and capacity level q is established; 0, otherwise

\( b_{kl}^{m} \):

1 if a hub link k–l with transportation mode m is established between hub nodes k and l

\( n_{kl}^{m} \):

Number of transportation vehicle m on hub link k–l

In the sequel, first the mathematical model for the HLP in a public transportation network with a static demand is proposed. Then, the elastic demand formulation and its application in the presented model are explained.

3.2.1 Mathematical formulation with a static demand

Considering the assumptions and notations in the previous sections, the bi-objective formulation for a public transportation network with a hub structure and static demand is developed as follows:

$$ \begin{aligned} {\text{Max}}\quad z_{1} & = \sum\limits_{i \in N} \sum\limits_{k \in H} \left( {p_{ik} z_{ik} + \sum\limits_{j \in N} {p_{kj} x_{ij}^{k} } } \right) + \sum\limits_{i \in N} \sum\limits_{k \in H} \sum\limits_{l \in H} \sum\limits_{m \in M} p_{kl}^{m} y_{ikl}^{m} - \left(\sum\limits_{i \in N} {\sum\limits_{k \in H} {\left( {\tilde{c}_{ik} z_{ik} + \sum\limits_{j \in N} {\tilde{c}_{kj} x_{ij}^{k} } } \right)} } \right.\\ & \left.\quad + \sum\limits_{i \in N} {\sum\limits_{k \in H} {\sum\limits_{l \in H} {\sum\limits_{m \in M} {\alpha^{m} \tilde{c}_{kl}^{m} y_{ikl}^{m} } } } } + \sum\limits_{k \in H} {\sum\limits_{m \in M} {\sum\limits_{q \in Q} {\tilde{F}_{k}^{mq} } } } h_{k}^{mq} + \sum\limits_{k \in H} {\sum\limits_{l \in H} {\sum\limits_{m \in M} {\tilde{I}_{kl}^{m} } } } b_{kl}^{m} + \sum\limits_{i \in N} {\sum\limits_{k \in H} {\tilde{J}_{ik} a_{ik} } } \right.\\ & \left.\quad + \sum\limits_{k \in H} {\sum\limits_{l \in H} {\sum\limits_{m \in M} {CV^{m} n_{kl}^{m} } } }\right) \\ \end{aligned} $$
(1)
$$ {\text{Min}}\quad z_{2} = \sum\limits_{i \in N} {\sum\limits_{k \ne i \in H} {\tilde{t}_{ik} a_{ik} + \sum\limits_{k \in H} {\sum\limits_{j \ne k \in N} {\tilde{t}_{kj} a_{jk} + \sum\limits_{k \in H} {\sum\limits_{l \ne k \in H} {\sum\limits_{m \in M} {(\tilde{t}_{k}^{m} + \beta^{m} \tilde{t}_{kl}^{m} } + \tilde{t}_{l}^{m} )b_{kl}^{m} } } } } } } $$
(2)

s.t.

$$ \sum\limits_{k \in H} {z_{ik} = \sum\limits_{j \in N} {W_{ij}^{\rm max } } } \quad \forall i \in N $$
(3)
$$ \sum\limits_{k \in H} {x_{ij}^{k} = W_{ij}^{\rm max } } \quad \forall i,j \in N $$
(4)
$$ \sum\limits_{l \ne k \in H} {\sum\limits_{m \in M} {y_{ikl}^{m} + \sum\limits_{j \in N} {x_{ij}^{k} } } } = \sum\limits_{l \ne k \in H} {\sum\limits_{m \in M} {y_{ilk}^{m} + z_{ik} } } \, \quad \forall i \ne k,k \in H $$
(5)
$$ z_{ik} \le a_{ik} \sum\limits_{j \in N} {W_{ij}^{\rm max } } \quad \forall i \in N,k \in H $$
(6)
$$ \sum\limits_{i \in N} {x_{ij}^{k} \le a_{jk} \sum\limits_{i \in N} {W_{ij}^{\rm max } } } \quad \forall k \in H,j \in N $$
(7)
$$ y_{ikl}^{m} + y_{ilk}^{m} \le b_{kl}^{m} \sum\limits_{j \in N} {W_{ij}^{\rm max } } \quad \forall k,l \in H,i \in N,m \in M $$
(8)
$$ a_{ik} \le a_{ii} + a_{kk} \quad \forall i \ne k \in N,k \in H $$
(9)
$$ a_{ik} \le 2 - a_{ii} - a_{kk} \quad \forall i \ne k \in N,k \in H $$
(10)
$$ \sum\limits_{k \in H} {a_{kk} = P} $$
(11)
$$ b_{kl}^{m} \le \sum\limits_{q \in Q} {h_{k}^{mq} } \quad \forall k \ne l,l \in H,m \in M $$
(12)
$$ b_{kl}^{m} \le \sum\limits_{q \in Q} {h_{l}^{mq} } \quad \forall k \ne l,l \in H,m \in M $$
(13)
$$ \sum\limits_{q \in Q} {h_{k}^{mq} } \le a_{kk} \quad \forall k \in H,m \in M $$
(14)
$$ \sum\limits_{m \in M} {b_{kl}^{m} \le 1} \quad \forall k \ne l,l \in H $$
(15)
$$ b_{kl}^{b} \le S_{kl} \quad \forall k \ne l,l \in H,b \in M $$
(16)
$$ \sum\limits_{i \in N} {z_{ik} } \le \sum\limits_{q \in Q} {\sum\limits_{m \in M} {h_{k}^{mq} {\tilde{\varGamma }}_{k}^{mq} } } \quad \forall k \in H $$
(17)
$$ \sum\limits_{i \in N} {y_{ikl}^{m} + } \sum\limits_{i \in N} {y_{ilk}^{m} } \le b_{kl}^{m} \tilde{R}_{kl}^{m} \quad \forall k,l \ne k \in H,m \in M $$
(18)
$$ V^{m} n_{kl}^{m} \ge \sum\limits_{i \in N} {y_{ikl}^{m} } \quad \forall k,l \ne k \in H,m \in M $$
(19)
$$ n_{kl}^{m} = n_{lk}^{m} \quad \forall k,l \ne k \in H,m \in M $$
(20)
$$ \sum\limits_{k \in H} {\sum\limits_{l \in H} {b_{kl}^{t} \ge \sum\limits_{q \in Q} {\sum\limits_{k \in H} {h_{k}^{tq} - 1} } } } \quad \forall t \in M $$
(21)
$$ y_{ikl}^{m} \ge 0,\,x_{ij}^{k} \ge 0,\,z_{ik} \ge 0,\,a_{ik} = \{ 0,1\} ,\,h_{k}^{mq} = \{ 0,1\} ,\,b_{kl}^{m} = \{ 0,1\} ,\,n_{kl}^{m} \ge 0 $$
(22)

The first objective (1) is to maximize the total profit of the transportation network. The first and second terms in function (1) represent the income of transportation between non-hub and hub nodes and between hubs, respectively. The third term includes the transportation costs between non-hub nodes and hub nodes, transportation costs between hub nodes, construction costs of hubs and hub links, the cost of allocating non-hub nodes to the hub nodes, and the total cost of transportation vehicles. The second objective function (2) minimizes the total transportation and waiting time in the network. Equation (3) ensures that all flows leave their origin. Equation (4) guarantees that the destination node j receives its flow originated from origin node i. Equation (5) is the flow conservation constraint. Constraints (6) and (7) ensure that if a non-hub node is allocated to a hub node, it can send (receive) its flow to (from) the assigned hub node. \( y_{ikl}^{m} \) can be positive only if the hub link between k and l with transportation mode m exists; this is ensured by Constraint (8). Constraints (9) and (10) indicate that if node i is allocated to hub node k by \( a_{ik} \), only one of these nodes must be considered as a hub. Equation (11) guarantees that exactly p hubs are chosen in the network. Constraints (12) and (13) ensure that both end-points of a hub edge must be hub nodes with the same transportation mode. A hub station with transportation mode m and capacity level q can be constructed only if node k is selected as a hub. This is guaranteed by Constraint (14). Constraint (15) specifies that only one transportation mode (i.e., BRT or metro) should be considered between two hub nodes. A BRT mode can be used in a hub link only if there is a street path between those hub nodes; This restriction is ensured by constraint (16). Constraints (17) and (18) restrict the input flow to a hub node and link, respectively. Constraints (19) and (20) determine the required number of vehicles to fulfil the existing demand on each hub link. Constraint (21) ensures that the metro hub network should be connected. For the sake of creating a connected hub network, the number of hub links should be at least equal to the number of links in a tree-shaped network. Therefore, the number of hub links to be used is at least the total number of hubs minus one. Constraint (22) represents the domains of the decision variables.

It should be noted that in this project, planners intend to connect all the significant population points using a unified and balanced urban transportation system, so the importance of all the points are equal and of the same importance. For this reason, the demands of nodes are not considered in the second objective function. However, in some cases, it is important to consider the weight of demands in each node for calculating the total travelling time of the network. In this case, the second objective function should be adapted to such a scenario by substituting \( a_{ik} \) by \( z_{ik} \) in the first term, \( a_{jk} \) by \( x_{ij}^{k} \) in the second term, and \( b_{kl}^{m} \) by \( y_{ikl}^{m} \) in the third term.

3.2.2 Mathematical formulation with an elastic demand

In this model, the demand (i.e., passengers) between every two nodes is the function of the utility. So it is shown as Wij(Ui + Uj), and it means that the demand between nodes i and j depends on the utility perceived by passengers at nodes i and j. The utility is assumed to be additive for simplicity of calculation. In this formulation, Ui means the utility perceived by node i, and it is calculated by:

$$ U_{i} = \sum\limits_{{k,d_{ik} \le D}} {\sum\limits_{m \in M} {\sum\limits_{q \in Q} {u_{ik}^{mq} } h_{k}^{mq} } } \quad \forall i \in N $$
(23)

where \( u_{ik}^{mq} \) is the utility of hub k with transportation mode m and capacity level q perceived by node i. As specified in Eq. (23), the utility perceived by each node is a function of the location of hubs within distance D from that node so the demand in each node is not defined in advance. As proposed by Khosravi and Akbari Jokar (2014), the function of the demand between every two nodes is defined as follows:

$$ W_{ij} (U_{i} + U_{j} ) = (W_{ij}^{\rm max } - W_{ij}^{\rm min } ) \cdot c_{ij} \cdot (U_{i} + U_{j} ) + W_{ij}^{\rm min } $$
(24)
$$ c_{ij} = \frac{1}{{U_{i}^{\rm max } + U_{j}^{\rm max } }} $$
(25)
$$ U_{i}^{\rm max } = \sum\limits_{{k,d_{ik} \le D}} {\sum\limits_{m \in M} {Max_{q} \{ u_{ik}^{mq} \} } } $$
(26)

where \( U_{i}^{\rm max } \) is the maximum utility perceived by passengers at node i. \( W_{ij}^{\rm max } \) and \( W_{ij}^{\rm min } \) are the maximum and minimum demands between nodes i and j, respectively. This formulation means that among nodes, the ones with high attractiveness should be established as hubs. \( u_{ik}^{mq} \) is a function of the attractiveness of the facility and a decaying function of the distance:

$$ u_{ik}^{mq} = \frac{{A_{k}^{mq} }}{{d_{ik} }} $$
(27)

where \( A_{k}^{mq} \) is the attractiveness of hub k with transportation mode m and capacity level q, and dik is the distance between node i and hub k. Therefore, by using Wij(Ui + Uj) instead of \( W_{ij}^{{\rm max} } \) in the formulation (1)–(22), the model with the elastic demand can be formulated as follows:

$$ \begin{array}{*{20}l} {\text{Max}} \hfill & {z_{1} } \hfill \\ {\text{Min}} \hfill & {z_{2} } \\ {{\rm s.t.}}& \hfill \\ \end{array} $$

Constraints (5), (9)–(22)

$$ \sum\limits_{k \in H} {z_{ik} = \sum\limits_{j \in N} {\left( {(W_{ij}^{{\rm max} } - W_{ij}^{{\rm min} } )c_{ij} \left( {\sum\limits_{{k,d_{ik} \le D}} {\sum\limits_{m \in M} {\sum\limits_{q \in Q} {(u_{ik}^{mq} + u_{jk}^{mq} )} } } h_{k}^{mq} } \right) + W_{ij}^{{\rm min} } } \right)} } \quad \forall i \in N $$
(28)
$$ \sum\limits_{k \in H} {x_{ij}^{k} = (W_{ij}^{{\rm max} } - W_{ij}^{{\rm min} } )c_{ij} \left( {\sum\limits_{{k,d_{ik} \le D}} {\sum\limits_{m \in M} {\sum\limits_{q \in Q} {(u_{ik}^{mq} + u_{jk}^{mq} )h_{k}^{mq} } } } } \right) + W_{ij}^{{\rm min} } } \quad \forall i,j \in N $$
(29)
$$ z_{ik} \le \sum\limits_{j \in N} {\left( {(W_{ij}^{{\rm max} } - W_{ij}^{{\rm min} } )c_{ij} \left( {\sum\limits_{{k^{{\prime }} ,d_{ik} \le D}} {\sum\limits_{m \in M} {\sum\limits_{q \in Q} {(u_{{ik^{{\prime }} }}^{mq} + u_{{jk^{{\prime }} }}^{mq} )} } } h_{{k^{{\prime }} }}^{mq} a_{ik} } \right) + W_{ij}^{{\rm min} } a_{ik} } \right)} \quad \forall i \in N,k \in H $$
(30)
$$ \sum\limits_{i \in N} {x_{ij}^{k} \le \sum\limits_{i \in N} {\left( {(W_{ij}^{{\rm max} } - W_{ij}^{{\rm min} } )c_{ij} \left( {\sum\limits_{{k^{{\prime }} ,d_{ik} \le D}} {\sum\limits_{m \in M} {\sum\limits_{q} {(u_{{ik^{{\prime }} }}^{mq} } } } + u_{{jk^{{\prime }} }}^{mq} )h_{{k^{{\prime }} }}^{mq} a_{jk} } \right) + W_{ij}^{{\rm min} } a_{jk} } \right)} } \quad \forall k \in H,j \in N $$
(31)
$$ y_{ikl}^{m} + y_{ilk}^{m} \le \sum\limits_{j \in N} {\left( {(W_{ij}^{{\rm max} } - W_{ij}^{{\rm min} } )c_{ij} \left( {\sum\limits_{{k^{{\prime }} ,d_{ik} \le D}} {\sum\limits_{{m^{{\prime }} }} {\sum\limits_{q} {(u_{{ik^{{\prime }} }}^{{m^{{\prime }} q}} + u_{{jk^{{\prime }} }}^{{m^{{\prime }} q}} )} } } h_{{k^{{\prime }} }}^{{m^{{\prime }} q}} b_{kl}^{m} } \right) + W_{ij}^{{\rm min} } b_{kl}^{m} } \right)} \quad \forall k,l \in H,i \in N,m \in M. $$
(32)

3.2.3 Linearization of the model

According to the definition of the demand function, the right hand side of constraints (30)–(32) are nonlinear. To avoid such complexities, the model is linearized on the basis of the method proposed in Mohammadi et al. (2011).

4 Solution methodology

There are fundamental differences between single and multi-objective optimization problems. In a single-objective one, the aim is to determine a unique optimal solution whereas in multi-objective optimization problems, a group of solutions according to the non-dominance criterion, named the Pareto set, is calculated. In this section, first, the uncertain input parameters of the presented model are proposed. Then, a hybrid two-stage approach is presented to modify the proposed bi-objective fuzzy model into the corresponding auxiliary crisp one.

4.1 Uncertain parameters

In HLPs, some parameters (e.g., costs, transportation time and demand) can change as the result of the uncertain environment. In this paper, the uncertain parameters include:

  • transportation costs (\( c_{ik} ,c_{kl}^{m} \)), fixed cost of establishing a hub station (\( F_{k}^{mq} \)), and hub and spoke links (\( I_{kl}^{m} ,J_{ik} \)) because these costs may be subject to governmental policies;

  • capacity of hub stations (\( \Gamma _{k}^{mq} \)) because it is affected by company policies;

  • capacity of hub links (\( R_{kl}^{m} \)) because it may be affected by the vehicles traffic;

  • transportation time and waiting time (\( t_{ik} ,t_{k}^{m} \)) because these parameters depend on the weather conditions or vehicles characteristics.

The above parameters are assumed to be fuzzy numbers with a triangular structure as shown in Table 1. All other parameters are assumed to be deterministic because they usually do not change a lot during a long period of time.

Table 1 Uncertain parameters
Full size table

4.2 Proposed fuzzy mathematical programming

In this paper, a hybrid two-stage method is presented for solving the proposed fuzzy bi-objective mixed-integer linear programming (FBMILP) model for small-sized problems. In the first stage, the original problem is converted into a crisp counterpart multi-objective function by means of an effective approach proposed by Jimenez et al. (2007). In the second stage, the TH method (Torabi and Hassini 2008) is used to convert the bi-objective functions into a mono-objective problem.

4.3 Equivalent auxiliary crisp model

To deal with the uncertain parameters, Jimenez et al. (2007) developed a method to convert a fuzzy MIP model into a crisp counterpart. The Jimenez et al.’s method allows to make a decision interactively with the decision maker (DM). Through the idea of a feasible optimal solution in degree \( \alpha \), the DM has enough information to fix an aspiration level. The DMs can also choose the degrees of feasibility that they are willing to admit depending on the context. It is important to highlight that the acceptable optimal solutions in degree \( \alpha \) are not fuzzy quantities, which makes it easier to take a decision in a simple way by solving a crisp parametric linear program. The DM also has additional information about the risk of violation of the constraints, and about the compatibility of the cost of the solution with his/her preferences for the objective function values.

This method is effective and capable of solving linear programming models because it does not change the linearity of the model and the number of objective functions. It can be applied to different kinds of membership functions. As the triangular fuzzy number (TFN) is the simplest and most popular one among various shapes of fuzzy numbers, it is assumed that fuzzy parameters are adapted to be TFNs and shown as \( \widetilde{c} = (c_{1} ,{\text{ c}}_{2} ,{\text{ c}}_{3} ) \) where c1, c2 and c3 are pessimistic, possible and optimistic values of the fuzzy parameters, respectively and can be obtained by the decision maker (see Fig. 1). The expected interval of a fuzzy number \( \tilde{c} \) (\( {\text{EI}}({\tilde{\text{c}}})) \), can be calculated by:

$$ {\text{EI}}({\tilde{\text{c}}}) = \left[ {{\text{E}}_{1}^{\text{c}} ,{\text{E}}_{2}^{\text{c}} } \right] = \left[ {\int_{0}^{1} {f_{A}^{ - 1} \left( \alpha \right)d\alpha } - \int_{0}^{1} {g_{A}^{ - 1} \left( \alpha \right)d\alpha } } \right] $$
(33)
$$ {\text{EV}}({\tilde{\text{c}}}) = \frac{{{\text{E}}_{1}^{\text{c}} + {\text{E}}_{2}^{\text{c}} }}{2} = \frac{{{\text{c}}_{1} + 2{\text{c}}_{2} + {\text{c}}_{3} }}{4} $$
(34)

where the expected value of a triangle fuzzy number (\( EV\left( {\tilde{c}} \right) \)) is the half point of its expected interval as per Eq. (34). In order to cope with fuzzy parameters in the constraints, the following constraints are presented:

$$ \left[ {\left( {1 - \alpha } \right)E_{1}^{{a_{i} }} + \alpha E_{2}^{{a_{i} }} } \right]x \le \alpha E_{1}^{{b_{i} }} + \left( {1 - \alpha } \right)E_{2}^{{b_{i} }} \quad \forall i = 1, \ldots ,m;\quad x \ge 0;\quad \alpha \in \left[ {0,1} \right]. $$
(35)
Fig. 1
figure 1

Triangular fuzzy number \( \tilde{c} \)

Full size image

Interested readers can refer to Pishvaee and Torabi (2010), Rabbani et al. (2018) and Zhalechian et al. (2018) for further explanations. Finally, the auxiliary crisp model counterpart can be formulated as:

$$ \begin{aligned} {\text{Min}}\quad OFV_{1} & = \sum\limits_{i \in N} {\sum\limits_{k \in H} {\left( {p_{ik} z_{ik} + \sum\limits_{j \in N} {p_{kj} x_{ij}^{k} } } \right) + \sum\limits_{i \in N} {\sum\limits_{k \in H} {\sum\limits_{l \in H} {p_{kl}^{m} y_{ikl}^{m} } } } } } - \left(\sum\limits_{i \in N} \sum\limits_{k \in H} \left(\left( {\frac{{c_{ik(1)} + 2c_{ik(2)} + c_{ik(3)} }}{4}} \right)z_{ik} \right.\right.\\ & \left.\left.\quad + \sum\limits_{j \in N} \left( {\frac{{c_{kj(1)} + 2c_{kj(2)} + c_{kj(3)} }}{4}} \right)x_{ij}^{k} \right) + \sum\limits_{i \in N} {\sum\limits_{k \in H} {\sum\limits_{l \in H} {\alpha^{m} \left( {\frac{{c_{kl(1)}^{m} + 2c_{kl(2)}^{m} + c_{kl(3)}^{m} }}{4}} \right)y_{ikl}^{m} } } } \right.\\ & \left.\quad + \sum\limits_{k \in H} {\sum\limits_{m \in M} {\sum\limits_{q \in Q} {\left( {\frac{{F_{km(1)}^{q} + 2F_{km(2)}^{q} + F_{km(3)}^{q} }}{4}} \right)} } } h_{k}^{mq} + \sum\limits_{k \in H} {\sum\limits_{l \in H} {\sum\limits_{m \in M} {\left( {\frac{{I_{kl(1)}^{m} + 2I_{kl(2)}^{m} + I_{kl(3)}^{m} }}{4}} \right)} } } b_{kl}^{m} \right.\\ & \left.\quad + \sum\limits_{i \in N} {\sum\limits_{k \in H} {\left( {\frac{{J_{ik(1)} + 2J_{ik(2)} + J_{ik(3)} }}{4}} \right)a_{ik} } } + \sum\limits_{k \in H} {\sum\limits_{l \in H} {\sum\limits_{m \in M} {CV^{m} n_{kl}^{m} } } }\right) \\ \end{aligned} $$
(36)
$$ \begin{aligned} {\text{Min}}\quad OFV_{2} & = \sum\limits_{i \in N} {\sum\limits_{k \ne i \in H} {\left( {\frac{{t_{ik(1)} + 2t_{ik(2)} + t_{ik(3)} }}{4}} \right)a_{ik} } } + \sum\limits_{k \in H} {\sum\limits_{j \ne k \in N} {\left( {\frac{{t_{kj(1)} + 2t_{kj(2)} + t_{kj(3)} }}{4}} \right)a_{jk} } } \\ & \quad + \sum\limits_{k \in H} {\sum\limits_{l \ne k \in H} {\sum\limits_{m \in M} {\left( {\left( {\frac{{t_{k(1)}^{m} + 2t_{k(2)}^{m} + t_{k(3)}^{m} }}{4}} \right) + \beta^{m} \left( {\frac{{t_{kl(1)}^{m} + 2t_{kl(2)}^{m} + t_{kl(3)}^{m} }}{4}} \right) + \left( {\frac{{t_{l(1)}^{m} + 2t_{l(2)}^{m} + t_{l(3)}^{m} }}{4}} \right)} \right)b_{kl}^{m} } } } \\ \end{aligned} $$
(37)

s.t.

$$ \sum\limits_{i \in N} {z_{ik} } \le \sum\limits_{q \in Q} {\sum\limits_{m \in M} {\left( {\alpha \left( {\frac{{\Gamma _{k(1)}^{mq} +\Gamma _{k(2)}^{mq} }}{2}} \right) + (1 - \alpha )\left( {\frac{{\Gamma _{k(2)}^{mq} +\Gamma _{k(3)}^{mq} }}{2}} \right)} \right)h_{k}^{mq} } } \quad \forall k \in H $$
(38)
$$ \sum\limits_{i \in N} {y_{ikl}^{m} } + \sum\limits_{i \in N} {y_{ilk}^{m} } \le \left( {\alpha \left( {\frac{{R_{kl(1)}^{m} + R_{kl(2)}^{m} }}{2}} \right) + (1 - \alpha )\left( {\frac{{R_{kl(2)}^{m} + R_{kl(3)}^{m} }}{2}} \right)} \right)b_{kl}^{m} \quad \forall k,l \ne k \in H,m \in M $$
(39)

Constraints (5), (28)–(32), (9)–(16), (19)–(22).

4.4 TH method

Various applicable methods have been proposed in the literature to solve crisp multi-objective programming models, such as a priori, interactive and posteriori methods (Vahdani et al. 2012; Zhalechian et al. 2017a, b; Akbari-Jafarabadi et al. 2017). One of the attractive methods is the fuzzy interactive method considering its capability in calculating the satisfaction level of each objective function regarding the decision maker’s priorities in an interactive manner (Zahiri et al. 2014). In this subsection, the two-phase interactive TH method is implemented to solve the multi-objective linear programming problem. Unlike classical multi-objective programming methods that might lead to inefficient solutions, the TH method promises to find efficient solutions. For further explanations, interested readers can refer to Torabi and Hassini (2008).

5 Solution algorithms

In this section, a meta-heuristic algorithm, namely multi-objective particle swarm optimization (MOPSO) is proposed to find near-optimal Pareto solutions of the given problem. Another multi-objective meta-heuristic algorithm, namely the non-dominated sorting genetic algorithm (NSGA-II) is also developed for the purpose of measuring the performance of the proposed MOPSO algorithm.

5.1 Solution representation

To illustrate the location of hubs and allocations of nodes in the network, a 2 × (n + p) matrix is formed, where n and p are the numbers of nodes and hubs, respectively. In the first row, the first n elements are random numbers between 0 and 1 with 2 decimal places, and the remaining p elements are random integers between 1 and n, which show the bits that should be chosen as hubs. For the second row, the first n elements are numbered in an ascending order according to random numbers in the same column above them. Each number corresponds to a node in the network, so every column represents a node in the network. The rest of the elements in the second row are filled in based on the number in the first row and represent the node number of the hubs. To allocate nodes to hubs, [the second digit of the random number in the first row modulo p] + 1 determines the number of allocations to that column (node); starting from the hubs positioned on the right of the node, hubs are allocated to the node, and if the number of hubs on the right is not sufficient, hubs on the left are allocated. For connecting hubs together, the same procedure is repeated this time detracting 1 from the number of allocations if it is equal to p (see Fig. 2 for an example). Additionally, a symmetric N × N matrix is formed, in which the (i − j)th element is chosen randomly between 1 or 2 for a hub link between i and j (satisfying the connectivity of the metro link [i.e., mode 2) and preexistence of streets for the BRT link (i.e., mode 1)], 3 is assigned to spoke links, and 0 means no link. Furthermore, the hub capacity levels are chosen after forming the network and calculating the flow each hub must support. To minimize the cost, the lowest capacity level supporting the flow is chosen as the capacity of hubs.

Fig. 2
figure 2

Continuous solution encoding for the proposed problem

Full size image

In Fig. 2, n = 10 and p = 3. For example, see column 5 corresponding to node 2 (determined by the number in the second row, column 5), the second digit of 0.72 is 2, 2 mod 3 is 2, so the number of allocations is 3. Hubs on the right are 6 and 7 (highlighted by blue color) and not enough, so hub 1 from the left is also assigned to node 2.

5.2 Particle swarm optimization

Particle swarm optimization (PSO) has been one of the popular algorithms, which has an efficient performance in a broad variety of optimization problems in comparison to lots of other meta-heuristics. This algorithm was first introduced by Eberhart and Kennedy (1995). It has a population, named particles, to pass through the search space. Each particle has a speed that works as an operator to create a new set of individuals. In real life problems, there are usually more than one objective functions to be optimized. Therefore, the single-objective algorithms cannot be used to solve the presented model and must be changed based on multi-objective assumptions. Details of the MOPSO algorithm are explained in Coello Coello et al. (2007).

5.3 Non-dominated sorting genetic algorithm

The non-dominated sorting genetic algorithm (NSGA-II) was first presented by Deb et al. (2002). The idea behind the algorithm is based on sorting and choosing the population fronts by non-dominance and crowding distance methods. In order to produce new solutions as the next offspring, the algorithm uses crossover and mutation operators. The next generation is created by combining the current population and generated offspring based on non-dominance and crowding distance method.

6 Computational experiments

In this section, we validate the proposed model through a case study based on the Qom Monorail Project (QMP) data. The conflict between two objective functions is shown based upon the results of the case study. To solve the bi-objective model for this case study, the TH aggregation function (see Sect. 4.4) is used for 15 nodes of the case study by GAMS software (24.0.1 version) using CPLEX as solver for the mixed-integer linear programs (MILP). Afterwards, various sensitivity analyses are carried out to analyze the sensitivity of the objective functions to the input parameters. In order to examine the model with large-sized problems, a performance comparison is carried out between the proposed meta-heuristics using the Turkish Data Set (TDS).

6.1 Case study

To validate the accuracy of the presented mathematical model, the model is tested on several numerical examples of the case study and the results are shown in this section. Additionally, the results of the elastic demand model are compared with the results of the static demand. The Qom city is located in the center of Iran with the area of 11,240 km2. However, this city is very populated because the presence of a holy shrine. Thousands of tourists travel to this city every day. Therefore, designing efficient public transportation for this population is of utmost importance. The data of this paper is taken from a real project executed in 2008 and evaluated by the experts. This city is divided into 15 traffic areas and the maximum demand of each area is estimated in the QMP. The centers of each area are selected as demand nodes.

The unit transportation costs (\( c_{ik} \)) and prices for spoke links are calculated based on the distances between the traffic areas. The transportation time between nodes are also calculated by assuming that the regular buses (used in spoke links) have a range of 40 km/h speed and BRT buses and metro trains have 50 km/h and 70 km/h speed, respectively.

Table 2 shows the input parameters used in this paper. In this table, parameters with the QMP indices are the parameters given from the QMP data and parameters with the TDS indices are the parameters available in the Turkish data set. \( \Gamma _{k}^{{mq_{1} }} \), \( \Gamma _{k}^{{mq_{2} }} \) and \( \Gamma _{k}^{{mq_{3} }} \)(N#P) show the capacity of the first, second and third levels of hubs with mode m (m = {t, b}) for a problem with N nodes and P hubs. \( R_{kl}^{m} \) (N#P) shows the capacity of hub link k–l with transportation mode m for a problem with N nodes and P hubs. Skl(0) means the elements in Skl matrix that have 0 value, in other words, the nodes with no streets connecting them. The results of the case study with three hubs are shown in Table 3 by both elastic and static demands and the same results are also shown in Figs. 3 and 4.

Table 2 Data used to solve the model in small and large-sized problems
Full size table
Table 3 Results of the QMP case study
Full size table
Fig. 3
figure 3

Hub topography for the case study with the elastic demand

Full size image
Fig. 4
figure 4

Hub topology for the case study with the static demand

Full size image

6.2 Sensitivity analysis

In this section, several problems are solved by changing the value of some important parameters to analyze the sensitivity of the model. γ and ϴh show the importance of the hth objective function and the coefficient of recompense \( (\sum\nolimits_{h} {\theta_{h} } = 1;\theta_{h} \succ 0) \). As shown in Table 4, changing the value of γ alters the values of both objective functions and the satisfaction degrees of each objective function correspondingly. Furthermore, greater value for ϴh indicate the importance of the hth objective function, so a higher weight should be assigned to it. Therefore, the value of that objective function should be improved. For example in this table, for a problem with 5 nodes and 2 hubs with elastic demand, by increasing ϴ1 from 0.4 to 0.5, the value of the first objective function increases from 3.13E+08 to 4.52E+08 (as it is a maximization objective function) and the value of the second objective function increases from 42.97 to 57.17 and becomes worse (as it is a minimization objective function). Consequently, greater values for ϴh indicates that higher weights are allocated to make a higher lower bound for the satisfaction degree of objectives (λ0) and subsequently more balanced compromise of solutions. The results show that using the elastic demand decreases the total profit and increases the total time of the network in comparison to the static demand. This happens because using the elastic demand model decreases the flow between every two nodes. This means that some of the potential demands of the nodes may be lost because the demand is a function of the utility perceived by each node. Thus, in the static demand model, the total profit is higher because all the passengers (\( W_{ij}^{{\rm max} } \)) of all the nodes use the hub network and increase the profit. However, it is well known that such static demand may not happen in practice. Moreover, in the elastic demand model, the location of hubs depend on the demand as the model tries to make a network that persuades passengers to use it and increase the total profit so the total time of the network increases in comparison to the static demand model.

Table 4 Summary of the results of the TH method using the case study data
Full size table

Figures 5 and 6 show the Pareto frontier for both static and elastic demands. Table 5 shows the effect of the capacity on both objective functions by changing the capacity from the value explained before for a problem with 7 nodes, 4 hubs, γ = 0.4, and ϴh = 0.4 and 0.6 in the model with the elastic demand. For each number of nodes, a special number of hubs is considered and represented as “number of nodes # number of hubs”; for example, 15#3 means 15 nodes and 3 hubs.

Fig. 5
figure 5

Pareto-optimal solutions for problem 10#4

Full size image
Fig. 6
figure 6

Pareto-optimal solutions for problem 7#4

Full size image
Table 5 Results of changing the capacity of hubs
Full size table

In Table 5, the Cap column shows the percentage of increasing or decreasing all the capacity level with respect to the initial value. The hub nodes column shows in each scenario which of the nodes are selected as hubs. Table 5, Figs. 7 and 8 show that an increase in the capacity level causes an increase in the total profit while the total time decreases. The reason is that by increasing the hub capacity, the amount of flow that can be assigned to the hub increases and it is not necessary to allocate demand nodes to further hubs because of capacity level restrictions. Therefore, the transportation cost decreases because of a decrease in the total distances and the result is that the profit increases. The total transportation time also decreases because of the allocation of demand nodes to closer hubs. As shown in Figs. 7 and 8, the total profit and time become constant from scenario 15. This means that by a 30% increase in the capacity level from the base value, the capacity constraint is not any more binding.

Fig. 7
figure 7

Effects of the capacity level on the total profit

Full size image
Fig. 8
figure 8

Effects of the capacity level on the total time in a network

Full size image

Table 6 shows the effect of passengers’ flow variation in solving the model with the elastic demand on the objective functions, hub nodes and hub links. The results suggest that by increasing the number of passengers of each node, both objective functions increase. That is because by increasing the passengers, the total income of the network increases, so the total profit increases. However, using more passengers from the hub network results in an increase in the total time of transportation. That is because having more passengers may exceed the capacity of closer hubs and the allocation of demand points to further hubs because of capacity limitation results in an increase in the total time of transportation. In Table 6, the values of objective functions, hub nodes and hub links in each scenario are shown for a problem with 7 nodes, 4 hubs, γ = 0.4 and ϴh = 0.4 and 0.6. Figures 9 and 10 show the effect of the demand changes in the total profit and total time.

Table 6 Objective functions variations according to demand changes
Full size table
Fig. 9
figure 9

Profit variations according to the demand changes

Full size image
Fig. 10
figure 10

Total time variations according to the demand changes

Full size image

6.3 Managerial insights

It is worthwhile to analyze the solutions in different scenarios based on the urban structure of Qom which is divided into four main areas. Some general yet useful information about these areas can be found in Table 7. For the static demand model, in the solutions that have higher profit than the average profit of all the Pareto solutions of the same problem the hubs are chosen to be inside areas 1 and 2 in 82% of the times, which seem to have the higher percentage of residential utilization, while if we analyze this for the solutions with shorter travel time, in 74% of the times, the hubs are situated in areas 3 and 4, which show a higher percentage of roads utilization. Also, for the model with elastic demand, for the solutions with higher profits, in 70% of the times, hubs are chosen to be in areas 3 and 2, which seem to show a good balance between percentage of commercial and residential utilizations, while for the solutions with better time, the hubs are still around 82% of the time inside areas 3 and 4 because of their high percentage of roads utilization. These findings can be interpreted below.

Table 7 General characteristics of four main areas of Qom
Full size table

Parts of the city with the highest flows happen to be the ones that include or are very close to dense residential areas, so in the static demand solutions, profit is increased if the hubs are constructed in those areas (considering that costs do not increase significantly if we choose such areas for hubs.) In the elastic model, commercial areas are also considered attractions in addition to residential areas, so this affects the utility of the hubs and consequently the demand for those hubs, which results in an increase in profit.

In both models, the solutions with shorter total travel time choose their hubs mostly in the areas with a higher percentage of roads utilizations because, in those areas, faster roads and streets and shorter paths help the network experience lower total transportation time.

One of the main challenges in studying and solving public transportation problems is the reliance on the parameters from the results of the previously conducted city planning studies. The ever-changing outlines of cities, their population and points of demand warrants a careful examination of the input parameters and analysis of sensitivity to these changes. Furthermore, the special dynamics of the city under study have to be fully understood to account for potentially significant changes that affect the optimality of the model’s solutions. Another challenge is the need for extra caution to avoid the definition of solutions that have negative effect on the population and result in under-served communities.

The results of the proposed model can serve as a good guideline when planning the development of new areas of the city. For example, if we consider an area that was not developed when the model was solved but that has similar characteristics of another area that was included into the model. The decision makers can make use of our results in order to plan ahead for the future development of the transportation for the new developing area (e.g., reserve an appropriate ground for possibly constructing one or more metro stations). Relating to the challenges, there is an opportunity to include other goals as well such as serving under-served communities and decide between optimal Pareto solutions from the lens of that goal.

6.4 Performance of the meta-heuristics

In this section, the performance of the proposed meta-heuristic algorithms is measured with some comparison metrics using the TDS. This data set includes 81 nodes. The other data, which are not in the TDS, are generated as shown in Table 2. The required parameters in the NSGA-II and MOPSO algorithms are tuned by the response surface methodology (RSM), as will be explained in the next subsection.

6.4.1 Parameters setting

The values of the algorithm’s parameters have a considerable effect on the performance of the algorithm. In this section, an appropriate tuning phase of the algorithm’s parameters is carried out in order to optimize their behavior. Therefore, the well-known RSP approach (i.e., Response Surface Methodology) is applied (Mohammadi et al. 2013). Table 8 shows the values of the parameters used while implementing NSGA-II and MOPSO methods.

Table 8 NSGA-II and MOPSO parameters settings
Full size table

6.4.2 Comparison metrics

To validate the performance of the presented MOPSO, four comparison metrics are employed. The performance of each algorithm can be calculated by the quality metric (QM). By using the mean ideal distance (MID) metric, the closeness of the solutions to the ideal point can be calculated. The diversification metric (DM) represents the spread of the solutions set. So, higher values of the DM mean better performance. The spacing metric (SM) estimates the uniformity of the spread and spacing of the solutions. It is easy to observe that lower values of this metric mean better performance. For further explanations, interested readers can refer to Mohammadi et al. (2013).

The proposed MOPSO algorithm is used to solve a number of test problems and its performance is compared against NSGA-II. Table 9 shows the QM and SM for small- and medium-sized problems. Likewise, Table 10 lists the DM and MID for small and medium-sized problems. Tables 11 and 12 show the four comparison metrics for large-sized problems. The values in the “Problem No.” column show the number of nodes versus the number of hubs. These tables show that the proposed MOPSO performs better in terms of the QM, DM, and MID metrics and NSGA-II has a better performance in the SM metric in most of the test problems. It means that MOPSO solutions have better quality compared to the NSGA-II (i.e., better QM). The spread of solutions in MOPSO is more than NSGA-II. Therefore, MOPSO Pareto solutions cover a wider range (i.e., better DM). Moreover, the MOPSO solutions are closer to the ideal solution (i.e., better MID). However, NSGA-II produces more uniform Pareto frontiers (i.e., better SM). Figures 11 and 12 show the higher performance of MOPSO in comparison with NSGA-II in obtaining better Pareto solutions.

Table 9 QM and SM for small and medium-sized problems
Full size table
Table 10 DM and MID for small and medium-sized problems
Full size table
Table 11 QM and SM for large-sized problems
Full size table
Table 12 DM and MID for large-sized problems
Full size table
Fig. 11
figure 11

Pareto-optimal solutions for problem 30#7

Full size image
Fig. 12
figure 12

Pareto-optimal solutions for problem 81#12

Full size image

7 Conclusions

In this study, a novel bi-objective mathematical programming model for the public transportation network design was proposed encompassing several real-world limitations and constraints that were mostly neglected in the literature. Some parameters were assumed to be fuzzy numbers because of the environmental impacts. An efficient approach for solving the fuzzy bi-objective mixed-integer linear programming, as proposed by Torabi and Hassini (2008), was implemented to cope with uncertainty. In addition, the presented model was tested using the real transportation data of the holy city of Qom. Several sensitivity analyses were applied to validate the effectiveness and efficiency of the proposed model, and managerial insights were provided. Moreover, two efficient meta-heuristics approaches based on multi-objective particle swarm optimization (MOPSO) and non-dominated sorting genetic algorithm (NSGA-II) are developed to address large-sized problems. The performance of the two algorithms has been compared through the use of four different metrics. The performance comparison has shown that the proposed MOPSO algorithm performs better than NSGA-II with respect to three metrics out of four.

For future work, the model can be extended using different elasticity functions. Furthermore, other meta-heuristics can be implemented and compared and more realistic constraints related, for example, to the roads congestion can be incorporated.