A Multilayer Network Approach for the Bimodal Bus–Pedestrian Line Planning Problem

Abstract

In this paper, we formulate and solve the urban line planning problem considering a multilayer representation of a bimodal transportation network. Classical formulations are usually constructed over a planar network, which implies the need to introduce several strong non-linearities in terms of frequencies when modeling transfer times. In the proposed network representation, each candidate line is stored in a specific layer and the passengers’ movements for each origin–destination pair are modelled considering a strategy subgraph, contributing to a sparse model formulation that guarantees feasibility and simplifies the assignment process. The methodology is first tested using the Mandl network, obtaining results that are comparable in terms of quality with the best metaheuristic approaches proposed in the scientific literature. With the aim of testing its applicability to large scenarios, the proposed approach is then used to design the main urban transit network of Seville, a large scenario with 141 nodes and 454 links, considering artificial unfavorable demand data. The reasonable computation time required to exactly solve the problem to optimality confirms the possibility of using the multilayer approach to deal with multimodal network design strategic problems.

1. Introduction

The line planning problem (LPP) is a critical strategic stage in the urban transit transportation planning process. Also known as the urban transit network design problem (UTNDP), it is concerned with the design of a set of transit routes and the determination of their frequencies with the aim of matching, in the most convenient way, the passenger demand and the supply; see, for instance, Goossens et al. [1], Guan et al. [2]. It is worth mentioning that several authors working in the field of passenger railway transportation have used the term line planning to specifically refer to the problem of determining the stopping pattern, frequency and capacity of different service types moving along a corridor (Wang et al. [3], Zhou et al. [4], Zhao et al. [5]). In this work, we focus on the conventional UTNDP, an NP-hard problem (Magnanti and Wong [6]), which, given its practical importance, has attracted the interest of researchers in the transportation field since the 1960s. For a comprehensive review of the different approaches used to solve this problem, we recommend that readers consult the chronologically ordered works of Ceder [7], Wirasinghe [8], Guihaire and Hao [9] and Kepaptsoglou and Karlaftis [10]. More recent reviews by Shöbel [11], Farahani et al. [12] and Ibarra-Rojas et al. [13] present a classification of the different variants of the problem, objective functions and solving procedures.

The objective functions generally considered in urban transit network design problems can be passenger- and/or operator-oriented (Guihaire and Hao [9]). From the point of view of the passengers, the designed network should cover the largest possible area of the scenario under study; the set of final lines should allow the maximum number of direct trips (that is, a lower number of transfers) and ensure the minimum travel time (Jiang et al. [14], Feng et al. [15]). From the operator point of view, the main goal is to reduce as much as possible the operating or even the construction costs (typically when dealing with urban railways), an objective that is essentially attained by minimizing the number of lines, their lengths and their frequencies. These aspects are usually considered in models where the demand is supposed to be inelastic. Otherwise, a pure cost reduction policy could imply an important loss of passengers (Robenek et al. [16], Canca et al. [17]). Multi-objective formulations consider both the user and operator points of view (Iliopoulou et al. [18]). The most common approach is to construct a weighted objective function trying to simultaneously minimize the user and operator costs. For example, Mauttone and Urquhart [19] proposed a multi-objective model considering the total travel time of users and the fleet size. A second type of multi-objective model contains bilevel formulations, as in the works by Gao et al. [20], Szeto and Jiang [21], Kim et al. [22] and Goerigk and Schmidt [23], where, usually, the upper-level problem is used to determine the lines and the frequencies, and the lower level is responsible for modeling the passenger route choice behavior for a given set of transit lines proposed by the upper level.

In general, the output of a transit network design problem is a line plan (also called a line concept) that should offer good connectivity (the possibility of traveling among major trip generator locations) and spatial coverage (geographic accessibility). The designed network has a direct influence on users and operators. On the one hand, users will expect a transit network with as many lines as possible to facilitate direct trips between selected points, minimizing then the need to transfer between lines, as well as high frequencies, minimizing the waiting times at stops (or platforms in the case of railways) and the transfer times between lines. If these goals are not achieved, the perceived quality of the service will be low. In this case, if possible, potential passengers will opt for alternative transportation modes, thus diminishing the revenue obtained by the operator company. On the other hand, more lines and higher frequencies imply higher operating costs. These opposite goals require service operators to achieve a trade-off between the quality of the service offered and the operating cost of the network.

The remainder of this work is organized as follows. Section 2 presents a detailed review of line planning problems, classifying the contributions into three groups according to the use or not of a pool of candidate lines and the type of algorithm considered to solve the problem. The section ends by positioning this research within the above scope. Section 3 contains the description of the problem addressed and explains the construction of the bimodal multilayer network that we will use to model the LPP. Section 4 presents the mechanism used to build the pool of candidate lines. Section 5 presents the motivation and the process adopted to construct the strategy subgraph for each OD pair. Using the multilayer network structure, Section 6 presents the proposed arc- and path-based formulations for the LPP. Section 7 contains the computational experiments, starting from a comparison with other approaches by using the Mandl network as a benchmark scenario. Later, the largest network is used with the aim of measuring the performance of exact procedures over our multilayer structure. Finally, Section 8 contains the main conclusions.

2. Literature on the Line Planning Problem

The literature on the LPP is rich. The line planning problem has attracted the attention of many researchers since the end of the 1960s. At an early stage in the research of the line planning problem, researchers proposed practical heuristic procedures to generate routes and determine line frequencies (see Lampkin and Saalmans [24], Bel et al. [25], Mandl [26], Furth and Wilson [27], Ceder [28], Ceder and Wilson [29], Van Nes et al. [30], Shih and Mahmassani [31]). For instance, Lampkin and Saalmans [24] proposed a sequential approach, consisting of a first phase where a heuristic algorithm iteratively added routes with the objective of maximizing the number of direct passenger trips. Then, in a second phase, a random greedy-based search procedure was used to obtain line frequencies for a given fleet without capacity limitations. Bel et al. [25] proposed a sequential approach composed of a heuristic procedure to select good candidate streets. Later, in a second stage, starting from the subset of streets, a maximal set of lines was generated. Finally, the frequencies were obtained using a gradient-based search heuristic with the objective of minimizing the passenger waiting time. Van Nes et al. [30] presented a heuristic method that started from a set of lines (the set of lines that allowed the highest number of direct trips) and later increased the frequency of each line while satisfying given budget and fleet size limits.

Since then, various approaches have followed in the literature to address the LPP and can be classified into three different groups. The first group is formed by works that model the problem as a mathematical programming problem starting from an a priori defined pool of candidate lines; see Bussieck et al. [32], Bussieck et al. [33], Goossens et al. [1] and Guan et al. [2]. The second group focuses its attention on the formulation of models that design the network from scratch. Starting from an underlying network, the models construct lines by combining edges. To this end, binary variables are used to determine whether edges belong to lines. Moreover, several constraints are needed to guarantee that the edges that are selected to define a line maintain the appropriate structure, i.e., connectivity between consecutive edges, no loops and total length bounds. This approach gives rise to very complex formulations that are difficult to solve using exact procedures; see, for instance, the works by Szeto and Wu [34], Szeto and Jiang [21] and Canca et al. [17,35]. At the beginning of the 1990s, the first applications of metaheuristic techniques to the LPP appeared—for instance, the works by Baaj and Mahmassani [36], Chakroborty et al. [37], Chakroborty et al. [38] and Chakroborty and Wivedi [39]. The contributions of this group usually address the line planning problem by starting from a small set of initial candidate lines and take advantage of the capability of metaheuristics to generate new candidate lines from the lines initially defined and to handle more complex objectives, which are difficult to explicitly formulate.

2.1. Formulations Based on a Pool of Candidate Lines

Concerning the mathematical formulations of the LPP, as mentioned before, two different approaches have been followed: models that consider an a priori defined candidate line pool and formulations that compose lines by using links from an underlying network. Among the works of the first group, Wan and Lo [40] proposed a mixed integer model to design a transit network with multiple routes while minimizing the sum of the operational costs of all the transit lines. The model does not consider passengers’ preferences, which affect the determination of line frequencies. The authors illustrate the proposed formulation on a small network with 10 nodes, 19 edges and 9 origin–destination (OD) pairs, which is solved using CPLEX. Borndörfer et al. [41] presented a multicommodity path-based network flow formulation for the line planning and frequency setting problem. The model aims to minimize a weighted combination of the total travel time of passengers and operating costs. In this formulation, passengers can be freely routed through a set of paths connecting the origin and destination of each OD pair, but no transfers between lines are considered. Using a quite similar formulation, Guan et al. [2] considered a line planning problem to simultaneously determine the line configuration (from a given predefined line pool that connects all stations) and the passengers’ assignment, while minimizing the total length of all lines and the total number of passengers transferring. Goossens et al. [1] proposed several operator-oriented models to deal with line planning problems where the lines can have different halting designs. The formulations were solved for small scenarios using CPLEX.

Cancela et al. [42] studied the bus transport network design problem incorporating a passenger assignment process based on the use of a trajectory graph (see Spiess and Florian [43]). The problem addressed consisted of determining a set of lines and their corresponding frequencies while minimizing the total travel time of all passengers subject to, among others, an upper limit on the total fleet size. Three sets of experiments were presented using different test instances (a small artificial instance, benchmark instances (Wan and Lo [40], Bagloee and Ceder [44]) and the Rivera network, Uruguay (Mauttone and Urquhart [19])), determining the initial line pool in a different way for each one of these experiments. Specifically, for the Rivera experiment, a network with 84 vertices, 143 edges and 378 OD pairs, the pool of candidate routes was generated by the Pair Insertion Algorithm described in Mauttone and Urquhart [19] and the model was solved using CPLEX. In this case, the authors achieved a relative MIP gap of 12% after 4 h of computation. Lee and Nair [45] proposed an optimization algorithm for the line planning problem based on bilevel programming that takes advantage of the structure of the problem together with the estimation of the demand range. The problem of conservativeness due to the estimation of the range of demands is mitigated by adopting a robust optimization approach.

Several authors, such as Baaj and Mahmassani [36], Israeli and Ceder [46], Pattnaik et al. [47], Bielli et al. [48], Ngamchai and Lovell [49], Chakroborty [50], Tom and Mohan [51], Lee and Vuchic [52], Mauttone and Urquhart [19], Zhao and Zeng [53], Cipriani et al. [54] and Gattermann et al. [55], have specifically addressed the construction of an initial pool of candidate lines, proposing different types of heuristics. Considering that the structure of lines is known, several authors have focused their efforts on solving the frequency setting problem under different perspectives, usually adding new characteristics. Hadas and Shnaiderman [56] presented a new approach to frequency setting that allows the use of stochastic information of the data and its corresponding costs. Canca et al. [57] proposed a mixed integer non-linear programming model (MINLP) to determine the optimal frequency of lines and their capacities in a railway rapid transit (RRT) network in which different lines could operate on common tracks. Gkiotsalitis and Cats [58] addressed the problem of determining the optimal headways with the inclusion of operational deviations. Canca et al. [59] dealt with the problem of determining the optimal set of frequencies in an RRT network considering that the operating costs are convex and non-linear as a consequence of the energy consumption of trains. In their approach, the operation cost at each track is a function of the train model characteristics, the total passenger load and the average train speed. Sun and Szeto [60] proposed a bilevel model to simultaneously determine fares and frequencies while maximizing the profit of the transit operator. The authors extended the stochastic user equilibrium assignment model to the elastic demand case. Gutiérrez-Jarpa et al. [61] proposed a model for the design of a rapid transit system. The main idea of this work consists of building segments within broad corridors to connect some vertex sets, and, at a second stage, assembling the segments into lines. The objective of the model is to maximize the captured demand while minimizing travel times. Herbon and Hadas [62] proposed a new approach that combined the interests of passengers and operators into a generalized newsvendor model. Waiting times and overcrowding costs were used to define the passenger point of view, whereas the operator point of view was modeled considering the vehicle size, lost sales and occupancy.

2.2. Constructive Approaches

The second group of modeling approaches contains formulations that do not use an a priori pool of candidate lines. In these models, the set of lines is constructed from scratch, starting from a base network that contains the set of links that can be part of the lines. In general, these models are complex in comparison with models that incorporate an initial set of possible lines, including more types of variables and constraints. Normally, the structure of these models is composed of two related parts, the first one responsible for modeling the structure of lines (topological constraints) and the second one dealing with the passenger assignment and the frequency setting problems. Szeto and Wu [34] proposed a non-linear integer formulation to minimize the weighted sum of the number of transfers and the total travel time, considering in-vehicle and waiting times. The effort made by the authors to model the waiting times in the change between lines deserves special attention. The authors applied the methodology to design the trunk bus network of the Tin Shui Wai area, Hong Kong, a network with 23 nodes and 41 edges.

Szeto and Jiang [21] proposed a bilevel transit network design problem to simultaneously determine the transit routes and their frequencies. The upper-level problem was formulated as a MINLP with the goal of minimizing the number transfers, while the lower-level problem tackled the assignment problem considering capacity constraints. The authors proposed a hybrid artificial bee colony (ABC) to solve the problem. Canca et al. [63] proposed a mathematical programming model integrating the design of the infrastructure network, the line planning problem, the determination of train capacities for each line and the fleet investment and personnel planning. The demand was considered elastic, and, therefore, a modal split model was used to determine the mode choice between the transit network and a competing mode. The objective consisted of maximizing the total profit, achieving a balance between the captured demand and the cost associated with the network construction and operation.

Canca et al. [35] presented an adaptive large neighborhood search (ALNS) algorithm to jointly handle the network design and the line planning problem, but including rolling stock and personnel planning issues. For a similar problem, Canca et al. [17] proposed an iterative procedure governed by an ALNS that, at each iteration, runs a branch-and-cut algorithm implemented in Gurobi, which solves the assignment and network operation problems. De-Los-Santos et al. [64] addressed the problem of designing a bimodal pedestrian–bus transit network, locating the bus stops and determining the set of bus lines. This work presented two formulations, which were compared using the $ϵ$-constraint method. Recently, from a theoretical point of view, Heinrich et al. [65] investigated the parameterized complexity of line planning problems when all simple paths in the network can be used as potential lines, focusing not only on the network size but also on additional parameters such as the tree width combined with the maximum degree and maximum frequency.

2.3. Metaheuristic Approaches

As mentioned in the first section, the complexity of the LPP problem has caused numerous researchers to try to solve it using metaheuristic techniques. Most of these approaches do not propose an explicit mathematical formulation of the treated problem and focus on explaining the objective function considered, the coding of the solutions and the structure of the algorithm to be used. In the case of small networks, such as the Mandl network (Mandl [26]), some of these implementations have obtained the best known results. Moreover, the fitness function can manage more complex objectives than in the case of exact methods, incorporating transfers in an easy way and allowing one to solve, in general, large scenarios. For a detailed review of the use of metaheuristics to solve the problem of transit route network design, we refer readers to the work of Iliopoulou et al. [18].

Chakroborty and Wivedi [39] used a genetic algorithm (GA) to find near-optimal solutions to the line planning problem and illustrated the algorithm using the Mandl network. Bielli et al. [48] also proposed a GA where each chromosome is evaluated by the calculation of performance indicators, which are obtained after analyzing the results of a transit assignment. Given information on passenger demand, fleet size and an underlying network, Zhao [66] proposed a methodology aimed at minimizing transfers and total users and maximizing service coverage. The authors used a stochastic global search procedure based on the combination of a genetic algorithm and simulated annealing (SA). Zhao and Zeng [67] applied a search scheme based on an integrated SA and GA method previously developed by the authors to find optimal transit network route layouts to minimize passenger transfers; see also Zhao and Zeng [68]. Later, Zhao and Zeng [69] presented a metaheuristic method for the design of the set of routes, determining the headway of each line and the timetable. The solving procedure is based on an iteratively local search procedure combined with SA, tabu search (TS) and greedy search. The authors illustrate their algorithm using first the Mandl network and comparing the results against those obtained by Mandl [26], Shih and Mahmassani [31] and Baaj and Mahmassani [36]. Fan and Mumford [70] solved the transit network design problem using simple hill climbing and simulated annealing algorithms. The authors obtained good results in comparison with previous methods for Mandl’s benchmark problems. Szeto and Wu [34] proposed a GA with a new solution representation scheme and specific genetic operators to improve the search capabilities. The algorithm incorporates a diversity control mechanism based on a new definition of the Hamming distance to avoid premature convergence. Bagloee and Ceder [44] proposed a heuristic methodology considering several of the major concerns of transit authorities, such as budget constraints, level of service and the attractiveness of the transit lines. The first step was the construction of a set of potential stops using a clustering approach. As a second step, the authors used a special shortest-path procedure to build a set of candidate routes, categorized by hierarchy (mass, feeder, local routes). Finally, a GA was launched over the set of candidate routes, incorporating budgetary constraints, until a good solution was found. Cipriani et al. [54] first proposed a heuristic route generation algorithm and then a GA to find a suboptimal set of routes and their frequencies. Nikolić and Teodorović [71] tried to maximize the number of passengers attended while minimizing the total number of transfers and the total travel times of all passengers served. The authors used a bee colony optimization (BCO) metaheuristic to solve the problem. Chew et al. [72] described a bi-objective model for the LPP with limitations on the number of lines, nodes per line and number of passengers transferring more than once. The authors proposed a GA that was tested on benchmark data. Jiang et al. [14] minimized the weighted sum of the number of transfers and the total travel times of passengers. The problem was solved by means of a hybrid improved artificial BCO algorithm. Nayeem et al. [73] dealt with the maximization of the number of attended passengers, minimizing the total number of transfers and the total travel times of all passengers served. The authors proposed a GA to solve the problem. Nikolić and Teodorović [74] considered a transit network design and frequency setting problem where the set of routes is determined by the selection and assembly of links selected from an underlying network. To solve the problem, the authors designed a BCO metaheuristic. Zhao et al. [75] used a memetic algorithm (MA) to determine the optimal route configuration and line frequencies for a urban transit network. As the objective, they proposed the minimization of the passenger cost and the reduction of the unsatisfied passenger demand. Ref. [76] applied a GA to the line planning problem, achieving the best reported results for the Mandl network.

Kim et al. [22] proposed an SA algorithm to solve a bilevel model where the upper-level problem controls the decisions corresponding to operators and the lower-level problem simulates the users’ trip assignment. The addressed problem considers two transportation modes, transit and private cars. The mode choice is selected in the lower level using a logit modal split set of constraints. The methodology is applied to the Mandl network. Feng et al. [15] proposed a GA to solve a transit route design problem with the objective of minimizing the total travel time of passengers. In their paper, the authors pay attention to the analysis of the time composition involved in transfers. Buba and Lee [77] proposed a differential evolution approach to address the network design problem. The objective consisted of minimizing passenger costs and unattended demand. Kim et al. [78] studied a bilevel transit route network design considering mobility and accessibility measures. To determine the frequency of each line, the authors used a GA in the upper level, whereas the lower level addressed the modal split and the traffic assignment. Katsaragakis et al. [79] proposed a cat swarm optimization (CSO) algorithm to design near-optimal routes for public transportation networks. Vlachopanagiotis et al. [80] presented a new approach based on an alternating-objective GA to achieve Pareto optimality between user and operator interests. Durán-Micco et al. [81] presented a transit design and frequency setting problem that incorporates additional aspects, such as discrete frequencies, the selection of terminal nodes from a specific subset and the construction of circular lines. The authors proposed a bi-objective MA that minimizes the average travel time of passengers and the fleet size. Ahern et al. [82] developed a multi-objective simulated annealing (MOSA) algorithm to solve the transit network design problem. The algorithm follows a three-phase search procedure. The first two phases, using a GA, deal independently with the passenger and operator objectives. The last phase is a multi-objective search that simultaneously works with the two objectives to improve the quality of the solutions. Iliopoulou et al. [83] proposed a VNS-based algorithm to solve the LPP. The performance of the algorithm was tested using the benchmark Mandl network and compared with several methods from the literature. Sunhyung et al. [84] proposed an algorithm that uses reinforcement learning (RL) for the simultaneous optimization of the number of bus routes, the route design and service frequencies. The algorithm was tested on the benchmark Mandl Swiss network.

Table 1. Chronological summary of scientific contributions for the LPP.

Title	Context	From Scratch	Line Pool	Single/ Multilayer	Routes	Frequencies	Mathematical Formulation	User Costs	Operation/ Fleet Costs	Construction Costs	Solving Procedure	Side Constraints	Application
Chakroborty and Wivedi [39], (2002)	Bus	X (starting from a line pool and recombining)		Single	X	X	No	Average in-vehicle travel time Percentage of users who can go directly from their origin to their destination Percentage of users making a single transfer Percentage of users transferring twice Percentage of users who cannot use the transit network to go from their origin to their destination	_	_	Genetic algorithm	_	Mandl’s network (15 nodes, 21 arcs)
Guan et al. [2], (2006)	Bus		X	Single	X	X	Yes	Minimizing the total length of all transit lines, total passenger in-vehicle travel time and total number of passenger transfers	_	_	Standard branch and bound	_	Simplified version of Hong Kong mass transit Railway (36 OD pairs, 9 nodes, 10 arcs)
Goossens et al. [1], (2006)	Railways		X	Single	X	X	Yes		Operation costs	_	CPLEX 7.5	Different halting patterns	NSRandstad network (122 nodes, 138 arcs)
Zhao and Zeng [69], (2008)	Bus	X (starting from a line pool and improving)		Single	X	X	No	Total user cost	_	Budget constraints	Integrated simulated annealing, tabu and greedy search algorithm	_	Mandl and Miami Dade transit-based networks (replanning)
Fan and Mumford [70], (2010)	Bus		X	Single	X	No	No	Weighted sum of the total travel distance accumulated over all passengers and the total number of transfers	_	_	Hill climbing and simulated annealing	_	Mandl’s network (15 nodes, 21 arcs)
Bagloee and Ceder [44], (2011)	Several transit modes		X	Single	X	X	No	Total saved generalized time with respect to no-transit-plan scenario	_	_	Three phase heuristic: location of stops + route generation + GA	_	Winnipeg and Chicago scenarios
Szeto and Wu [34], (2011)	Bus	X (starting from a line pool and recombining)		Single	X	X	Yes	Weighted sum of the number of transfers and total passengers’ travel time	_	_	Genetic algorithm + frequency setting heuristic based on neighborhood search	_	Trunk bus network of Tin Shui Wai residential area, Hong Kong (23 nodes, 41 arcs)
Hadas and Shnaiderman [56], (2012)	Bus			Single		X	Yes (Analytical solution)	Overload and un-served demand	Empty seats	_	Analytical	Stochastic demand capacity	A line with 5 stops
Nikolić and Teodorović [71], (2013)	Bus	X (starting from a line pool and recombining)		Single	X		No	Number of satisfied passengers Total number of transfers Total travel time of all served passengers	_	_	Bee colony	_	(110 nodes, 275 arcs)
Nayeem et al. [73], (2014)	Bus	X (starting from a line pool and recombining)		Single	X		No	In-vehicle time of all served passengers Total number of transfers Total number of unsatisfied passengers	_	_	Genetic algorithm with elitism	_	Mandl’s network (15 nodes, 21 arcs) Yubei (70 nodes, 210 arcs) Brighton (110 nodes, 385 arcs) Cardiff (127 nodes, 425 arcs)
Nikolić and Teodorović [74], (2014)	Bus	X (starting from a line pool and recombining)		Single	X	X	No	In-vehicle time of all served passengers Total number of transfers Waiting time of passengers	_	_	Bee colony	_	Mandl’s network (15 nodes, 21 arcs)
Zhao et al. [75], (2015)	Bus	X (starting from a line pool and applying several operators)		Single	X	X	Yes, but not in a closed formulation	Minimizing the passenger cost Reducing the unsatisfied passenger demand	_	_	Memetic algorithm	_	Mandl’s network (15 nodes, 21 arcs)
Canca et al. [57], (2016)	Railway Rapid Transit		Fixed	Single		X	Yes	Travel time per transit user Waiting time Transfer time Fare	Energy consumption and maintenance Fixed operation cost Variable operation cost Fleet acquisition cost	_	Extended cutting plane (ECP) method	Specific shared segment constraints Capacity	Madrid metropolitan railway (33 stations, 76 arcs, 800 OD pairs)
Canca et al. [63], (2016)	Railway Rapid Transit	X		Single	X	X	Yes	Total user travel time (indirectly through profit maximization)	Fixed operation cost Variable operation cost Fleet acquisition cost Revenue	Network construction cost	Branch and bound	Competing transportation mode	Small artificial instance
Kim et al. [22], (2016)	New line (different mode alternatives) into an existing bus network		X	Single	One	X (only for one line)	Yes	Travel time cost of auto and transit users	Operating costs of auto and transit users	Cost of constructing a new transit line	Simulated annealing	_	Mandl’s network (15 nodes, 21 arcs)
Goerigk and Schmidt [23], (2017)	Railways		X	Single	X	X	Yes (Two formulations)	Total travel time	_	_	Exact for artificial scenarios (CPLEX) + Genetic algorithm for the real case	_	Long-distance railway network of Germany (250 stations, 652 arcs)
Canca et al. [35], (2017)	Railway Rapid Transit	X		Single	X	X	Yes	Travel time per transit user Waiting time Transfer time	Fixed operation cost Variable operation cost Fleet acquisition cost Revenue	Network construction cost	ALNS	Alternative mode	Seville metropolitan area (49 nodes, 135 arcs, 2352 OD pairs)
Canca et al. [59], (2018)	Railway Rapid Transit		Fixed	Single		X	Yes	Travel time per transit user Waiting time Transfer time	Total variable cost due to energy consumption Fleet operation costs		Sequential optimization	Vehicle selection	Metro network of Lisbon (50 stations, 54 arcs, 2450 OD pairs)
Buba and Lee [77], (2018)	Bus	X (Generating the line pool as in Mundford 2013)		Single	X	X	Yes	Total travel time Unsatisfied demands	Bounding fleet size	_	Differential evolution	Assuming maximum two transfers	Mandl’s network (15 nodes, 21 arcs)
Canca et al. [59], (2018)	Railway Rapid Transit	X		Single	X	X	Yes	Travel time per transit user Waiting time Transfer time Fare	Fixed operation cost Variable operation cost Fleet acquisition cost Revenue	Network construction cost	ALNS	Alternative mode	Artificial network (100 nodes, 275 arcs, 9900 OD pairs)
Kim et al. [78], (2019)	Mixed (one fixed rail line, variable bus lines)		X	Single	Only bus	X	Yes	User cost (Total trip time)	Operator costs	_	Genetic algorithm	Network re-adaptation with equity aspects	Adapted Mandl network
Canca et al. [17], (2019)	Railway Rapid Transit	X		Single	X	X	Yes	Travel time per transit user Waiting time Transfer time Fare	Fixed operation cost Variable operation cost Fleet acquisition cost Revenue	Network construction cost	ALNS Metaheuristic	Alternative mode	Seville metropolitan area (49 nodes, 135 edges, 2352 OD pairs)
Katsaragakis et al. [79], (2020)	Bus		X (up to 8 routes)		X	X		Average travel time per transit user Percentage of satisfied demand without any transfers Percentage of unsatisfied demand	_	_	Cat swarm optimization based algorithm	_	Mandl’s network (15 nodes, 21 arcs)
De-Los-Santos et al. [64], (2021)	Bus	X		Two modes into a single graph	X	X	Yes	In-bus travel time Waiting time Walking time Transfer time	Constraints on the length of lines and number of stops per line	_	GAMS, using CPLEX 12.5	Circular lines Lazy constraints (sub-tour elimination constraints)	Interurban eastern bus network of Seville (43 nodes, 44 arcs, 1806 OD pairs)
Lee and Nair [45], (2021)	Bus		X	Two modes into a single graph	X		Yes	Total system travel time Total travel time and penalty for over-capacity on the link	_	_	Column generation using CPLEX	Robustness	Abidjan, Côte d’Ivoire (15,000 OD pairs)
Vlachopanagiotis et al. [80], (2021)	Bus		X	Single	X	X	Schematic	Average travel time per transit user Percentage of demand satisfied without any transfers Percentage of unsatisfied demand	Fleet size	_	Alternating-objective genetic algorithm	Capacity	Mandl’s network (15 nodes, 21 arcs)
Durán-Micco et al. [81], (2022)	Bus	X (starting from a line pool and recombining)		Single	X	X	No	Average travel time	Fleet size	_	Bi-objective Memetic algorithm	Circular lines	271 nodes, 470 arcs, 16,823 OD-pairs
This paper	Bus		X	Many layers (Pedestrian and bus superlayer with as many layers as candidate lines)	X	X	Yes	In-bus travel time Waiting time Walking time Transfer time	Constraints on the length of lines Number of lines	_	Branch and cut	Strategy subgraph for each OD pair Specific connections	Seville urban network (140 nodes, 454 links, 19,440 OD pairs)

considers several of the reviewed approaches in chronological order, describing their main characteristics.

2.4. Discussion and Contributions

From all these contributions, regarding the ability to address large-sized instances, we highlight the works by Zhao and Zeng [69], Bagloee and Ceder [44], Cipriani et al. [54] and Durán-Micco et al. [81], where the authors illustrate the applicability of their algorithms using real-life scenarios. In general, with the exception of the comparisons to benchmark networks, such as the Mandl one, the results of the application of metaheuristics to the line planning problem do not establish a measure of the quality of the obtained solutions, mainly as a consequence of the lack of comparisons with the exact solutions of the addressed problems. Most of the contributions proposing a metaheuristic algorithm to solve the problem do not incorporate an explicit mathematical formulation. Moreover, since metaheuristics are nondeterministic approaches, several runs are needed to analyze the average effectiveness of the approach, which results in different final solutions, and imply a large increment in the computation time, which, in real networks, can be extended to several hours; see, for instance, the four papers mentioned above.

In general, the majority of the referred works address the LPP considering a single-layer network, which leads to complex models if an explicit formulation of transfer times is pursued. The need to explicitly consider transfers in modeling LPPs was previously discussed by Cancela et al. [42], who stated that if the passenger flow assignment is performed ignoring the choice of individual routes, the solution can include a large number of transfers for certain OD pairs. The difficulty in modeling waiting and transfer times lies in the non-linear dependency with respect to the line frequencies. Although the passenger assignment can be formulated as a linear problem (Spiess and Florian [43]), when this problem is part of a transit network design model, non-linearities emerge because the frequencies are then decision variables.

According to Camporeale et al. [85], there is still an important gap in modeling and analyzing strategic and tactical transportation problems when more than one transportation mode is considered. In this paper, we reformulate the classical LPP (Guan et al. [2], Szeto and Wu [34]) considering an initial pool of candidate lines but using a multilayer network structure, which allows us to work with more than one transportation mode and simplify the formulation of transfer times. In fact, we are dealing with a bimodal LPP where the pedestrian mode is considered together with the bus mode.

Our network construction process allows for the inclusion of mode-dependent links, as, for instance, rail tracks cannot be used by buses or pedestrians. The inclusion of a pedestrian layer ensures feasibility even if the capacity (both in terms of the number of vehicles and the maximum frequency allowed) is not sufficient, and makes the network slightly different from the Change and Go network representation proposed by Goerigk and Schmidt [23]. Moreover, differently from the previous reference, we consider that not all the links in the network must be covered by lines (which is the usual situation when the number of lines to be constructed is limited by a given budget), thus incorporating the possibility of not attending to the whole demand and including also the usual condition that bounds the line frequencies with maximum acceptable values (Guan et al. [2]), which is of particular relevance in the case of designing railway lines. In our approach, transfers are not modeled as simple penalties but ultimately depend on the frequencies of lines.

As a summary, the proposed methodology consists of formulating a line planning problem over a multilayer bimodal network where the pedestrian network is supported by one layer and the set of candidate lines is stored in a superlayer—each candidate line in a different layer—which is appropriately connected to represent feasible transfers in the network, so that the passenger demand can move freely, according to the travel, waiting and transfer times, between the origin and the destination of each OD pair. In these movements, the use of arcs is possible only if they belong to lines that are activated. Moreover, unlike other formulations, the inclusion of the pedestrian layer ensures feasibility even if the number of allowed lines or their frequencies and capacities are not sufficient to meet the entire demand of passengers. Regarding application to large scenarios, to reduce the size of the mathematical models, for each OD pair, the movement of passengers is reduced to a set of reasonable paths. We propose two sparse formulations, using the strictly necessary variables. The resultant models are formulated using a Python object-oriented multilayer network library specifically developed to this end, and the application program interface of the branch and cut algorithm included in the Gurobi solver.

3. Problem Description and Bimodal Multilayer Network Construction

Given an OD demand matrix and a graph $G=\{N,A\}$ representing the set of intersections and streets of an urban scenario, where N is the set of nodes and A is the set of directed arcs connecting the nodes in N, we want to determine the set of loopless lines that allow the highest possible number of passengers to move through the network with the minimum possible travel time. We assume that all the streets of the scenario can be traveled by pedestrians and that only a subset of them (which could coincide with the totality) can be used for the movement of buses; then, in the next stage, we use the pedestrian network as the base network to construct the bimodal multilayer network structure. We are also interested in simultaneously determining the necessary frequencies of lines and incorporate several service operator constraints. Let $G_P=\{N_P,A_P\}$ be a directed graph defining the pedestrian network of an urban scenario, where $N_P$ is the set of nodes, and $A_P$ is the set of directed pedestrian links connecting nodes. Note that the elements in $N_P$ could represent intersections among city streets or even simply demand points connected to the real street network by means of connectors. The pedestrian graph defines the first layer in a multilayer network that contains the complete bimodal transportation network. Let $N_B$ be a subset of elements in $N_P$ that define possible bus stops and $A_B$ a set of directed links connecting nodes in $A_P$. Then, $G_B=\{N_B,A_B\}$ is the graph defining the set of streets that can support the movements of buses. In most cases, arcs in $A_B$ are also included in $A_P$ since, in general, the streets allowing the allocation of a bus line also allow pedestrian movements.

Let $L=\{1,2,3,\dots ,\mathcal{l},\dots ,L\}$ be a set of potential lines defined in $G_B$. The construction process of L will be detailed in the next Section. A generic line $\mathcal{l}\in L$ is defined by two sets $N_{\mathcal{l}}=\{i:i\in \mathcal{l}\}$ and $A_{\mathcal{l}}=\{(i,j):(i,j)\in \mathcal{l}\}$ that contain, respectively, the nodes and arcs belonging to the line. Specifically, the set $A_{\mathcal{l}}$ contains arcs in both directions connecting each pair of nodes in $N_{\mathcal{l}}$.

Each of the candidate lines will be stored in a different layer. The union of all the bus layers will result in a superlayer that will store the complete line pool of candidate lines. Note that this structure could be extended with new superlayers for new transportation modes, such as, for instance, a metro or a tram system.

To access the elements, nodes and arcs of each of the lines in the bus superlayer, consider a code function ${\varphi }_{\mathcal{l}}$ that, applied to each node $i\in \mathcal{l}$, returns a unique node code that univocally represents a copy of i in the corresponding line, i.e., ${\varphi }_{\mathcal{l}}\left(i\right)$ is the copy of node $i\in \mathcal{l}$ in the layer of line ℓ in the bus superlayer. Let ${\mathcal{N}}_{\mathcal{l}}$ be the set of nodes of line ℓ in the layer number ℓ of the bus superlayer, ${\mathcal{N}}_{\mathcal{l}}=\{{\varphi }_{\mathcal{l}}\left(i\right):i\in \mathcal{l}\}$, and let ${\mathcal{A}}_{\mathcal{l}}$ be the set of arcs in line ℓ corresponding to the $\mathcal{l}-th$ layer in the bus superlayer, ${\mathcal{A}}_{\mathcal{l}}=\{({\varphi }_{\mathcal{l}}\left(i\right),{\varphi }_{\mathcal{l}}\left(j\right)):(i,j)\in \mathcal{l}\}$. In this way, ${\mathcal{N}}_{\mathcal{l}}$ and ${\mathcal{A}}_{\mathcal{l}}$ contain a coded copy of the nodes and arcs belonging to the sets $N_{\mathcal{l}}$ and $A_{\mathcal{l}}$, respectively. To easily manage the correspondence among nodes and arcs in the bus superlayer and the graph $G_B$, it is convenient to consider an operator ${\varphi }^{-1}$ that, applied to a node ${\varphi }_{\mathcal{l}}\left(i\right)\in {\mathcal{N}}_{\mathcal{l}}$, returns the original node in $G_B$, ${\varphi }^{-1}\left({\varphi }_{\mathcal{l}}\left(i\right)\right)=i$. In a similar manner, ${\varphi }^{-1}({\varphi }_{\mathcal{l}}\left(i\right),{\varphi }_{\mathcal{l}}\left(j\right))=(i,j)$ returns the projection of a link in ${\mathcal{A}}_{\mathcal{l}}$ to the corresponding arc in $A_P$ or $A_B$, as appropriate.

Each node i in $N_B$ and each arc in $A_B$ have a set of coded copies in the bus superlayer, as many as the number of lines to which it belongs. We can then define the sets $\psi \left(i\right)=\{{\varphi }_{\mathcal{l}}\left(i\right):\mathcal{l}\in L\}$ and $\psi (i,j)=\{({\varphi }_{\mathcal{l}}\left(i\right),{\varphi }_{\mathcal{l}}\left(j\right)):(i,j)\in \mathcal{l}\}$ containing all the coded replications of node i and arc $(i,j)$, respectively, in the bus superlayer. Figure 1 summarizes the notation in a simple case with only two candidate lines. The pedestrian layer contains nodes in $N_P$ and $N_B$; bus stops are represented with a filled circle. The set $\psi \left(10\right)$ is represented in more detail in the right-hand part of the figure.

To model transfers between lines, two additional sets are built for each node $i\in N_B$. The set ${\mathcal{T}}_b\left(i\right)=\{(i,{\varphi }_{\mathcal{l}}\left(i\right)):\mathcal{l}\in L\}$ contains boarding links, whereas ${\mathcal{T}}_a\left(i\right)=\{({\varphi }_{\mathcal{l}}\left(i\right),i):\mathcal{l}\in L\}$ is composed of alighting links. These two sets are represented in the right-hand part of Figure 1 for the node 10 using solid and dashed lines, respectively. Note that, alternatively, the sets ${\mathcal{T}}_b$ and ${\mathcal{T}}_a$ can also be defined starting at nodes j belonging to ${\mathcal{N}}_{\mathcal{l}}$, concretely, ${\mathcal{T}}_b\left(i\right)=\{({\varphi }^{-1}\left({\varphi }_{\mathcal{l}}\left(j\right)\right),j):\mathcal{l}\in L\}$ and ${\mathcal{T}}_a\left(i\right)=\{(j,{\varphi }^{-1}\left({\varphi }_{\mathcal{l}}\left(i\right)\right)):\mathcal{l}\in L\}$. The union of sets ${\mathcal{T}}_b\left(i\right)$ and ${\mathcal{T}}_b\left(i\right)$ for all nodes i in $N_B$, ${\mathcal{T}}_b={\cup }_{i\in N_B}{\mathcal{T}}_b\left(i\right)$ and ${\mathcal{T}}_a={\cup }_{i\in N_B}{\mathcal{T}}_b\left(i\right)$, contains, respectively, all the transfer arcs for boarding and alighting in the bimodal multilayer network.

To complete the construction of the bimodal multilayer network, it is convenient to store the set of lines that traverse a specific node in $N_B$. Since the counterpart nodes of node $i\in N_B$ in the superlayer bus network are given by $\psi \left(i\right)$, and each line is stored in an individual layer inside the bus superlayer, the number of lines traversing a bus node can be defined by the set $\mathcal{L}\left(i\right)=\{\mathcal{l}\in L:{\varphi }_{\mathcal{l}}\left(i\right)\in \psi \left(i\right)\}$. To easily recover the line traversing a link or a node in the bus superlayer, we also define an operator ${\mathcal{L}}^{-1}$ that, applied to a specific arc $({\varphi }_{\mathcal{l}}\left(i\right),{\varphi }_{\mathcal{l}}\left(j\right))\in {\mathcal{A}}_{\mathcal{l}}$, returns the line that is traversing the link, i.e., ${\mathcal{L}}^{-1}({\varphi }_{\mathcal{l}}\left(i\right),{\varphi }_{\mathcal{l}}\left(j\right))=\mathcal{l}$ and, applied to a node $({\varphi }_{\mathcal{l}}\left(i\right)$, returns the appropriate line, ${\mathcal{L}}^{-1}\left({\varphi }_{\mathcal{l}}\left(i\right)\right)=\mathcal{l}$).

We can now finally define the full bimodal multilayer network as a directed graph $G_{ML}=\{N_{ML},A_{ML}\}$, composed of the set of nodes $N_{ML}=\{N_P\cup \left({\cup }_{\mathcal{l}\in L}{\mathcal{N}}_{\mathcal{l}}\right)\}$ and the set of arcs $A_{ML}=\{A_P\cup \left({\cup }_{\mathcal{l}\in L}{\mathcal{A}}_{\mathcal{l}}\right)\cup {\mathcal{T}}_b\cup {\mathcal{T}}_a\}$.

4. Line Pool Construction

In the previous section, we considered the existence of a pool of candidate lines L, which was used to build the bimodal multilayer network. This section details the procedure to obtain the line pool. To start the construction of L, we will use the graph $G_B=\{N_B,A_B\}$. First, for every pair of nodes $i,j\in N_B:i\ne j$, we will compute the $k_L$ shortest paths, $k_L$ being a design parameter. Every shortest path with a length between two parameters $L_{min}$ and $L_{max}$ is stored in an intermediate line pool $\widehat{L}$. By default, $L_{min}$ is set to half the diameter of the graph $G_B$, allowing relatively short lines, whereas $L_{max}$ is set to two times the diameter, thus generating lines that can completely traverse the scenario. Note that these thresholds may be set to different values depending on the knowledge and preferences of the service operator company. This process generates a large number of candidate lines. For instance, taking as a reference the Mandl network (which is formed of 15 nodes and 21 edges) and fixing $k_L$ to a value of 10, the intermediate pool contains 876 lines. This number can increase significantly for larger scenarios, such as in the real-sized illustration used in the final part of this paper. At the same time that a candidate line is stored in the intermediate line pool $\widehat{L}$, we accumulate the demand between each pair of nodes belonging to the line considering both directions. Lines with a high accumulated demand are interesting candidate lines, since they can move directly a high number of passengers without the need of transferring. Simultaneously, we count and store the number of lines traversing each node in the network. The intermediate pool of lines $\widehat{L}$ is then ordered decreasingly, according to the accumulated demand. To complete the construction of the pool of candidate lines L, we remove the worst lines from the intermediate pool of candidate lines (those whose accumulated demand is below $\vartheta \%$ of the maximum accumulated demand). In this procedure, we take care not to remove a line if the lowest value of the counter of lines traversing the nodes of the line falls below a certain predefined value. In this way, the candidate lines of the final pool will cross all nodes of the network. The final number of lines in the pool will be $N_B\times (N_B-1)\times k_L\times (1-\vartheta )$, being $\vartheta$ the parts per unit of candidate lines to be removed from $\widehat{L}$, as mentioned above. Following, for example, the illustration of the Mandl network, considering $\vartheta =0.5$, the line pool will contain 351 candidate lines. Once the line pool is defined, the construction of the Mandl bimodal multilayer network produces a graph with $N_{ML}=2516$ nodes and $A_{ML}=9344$ arcs.

5. The Strategy Subgraph for Each OD Pair

Before proceeding to the model formulation for the bimodal pedestrian–bus line planning problem and following the goals of guaranteeing full feasibility and problem solvability, it is convenient to first introduce the concept of the strategy subgraph for each origin–destination (OD) network pair.

Assume (as usual in the formulation of the line planning problem) the existence of an origin–destination demand matrix $\mathcal{D}$, which normally corresponds to the peak hour interval of a design day. Each element $d_{ij}$ in the matrix represents the number of passengers willing to travel between nodes i and j in $G_P$. Recall that $N_P$ contains all the nodes acting as origins and destinations of trips (regardless of whether all of them are origins and destinations). Passengers will start their trips at a certain node i in $N_P$ and will follow a sequence of movements over pedestrian or bus arcs using $G_{ML}$, probably using intermediate transfer arcs to finish their journey at a different node j in $N_P$. For each OD pair, the sequence of traversed arcs will depend on the passenger’s choice, which is influenced by the final selected lines from the candidate line pool, as well as by their frequencies.

Let $\mathcal{W}$ be the set of OD pairs in the network, indexed by w, so that $w_o$ and $w_d$ represent the origin and destination of pair $w\in \mathcal{W}$, respectively. However, in real circumstances, it is unlikely that a passenger corresponding to a certain OD pair uses arcs located outside a set of “reasonable” paths connecting the origin and the destination of the OD pair, as in Canca et al. [57]. In order to construct a model that is as sparse as possible, for each OD pair w, we will define a directed subgraph $G_w$ containing the possible subset of nodes and arcs used by passengers to perform their trips. To this end, for each OD pair w, we first compute the $k_P$ shortest path between $w_o$ and $w_d$ using the graph $G_{ML}$ that defines the bimodal multilayer network. The higher the value of $k_P$, the higher the possibility of moving from $w_o$ to $w_d$ and the higher the number of variables and constraints in the line planning formulation. To compute the $k_P$ shortest paths, arcs in $A_P$ are weighted by their corresponding travel times, considering a pedestrian speed $v_P$; arcs in $A_{\mathcal{l}}$ are weighted by their travel times assuming a certain bus speed $v_B$; travel times in arcs belonging to ${\mathcal{T}}_a$ are set to 0; and travel times in boarding arcs belonging to ${\mathcal{T}}_a$ are set to half of the headway corresponding to the lowest allowed frequency value.

We will denote by $N_w$ and $A_w$ the set of nodes and arcs in $G_w$, respectively. Furthermore, to facilitate the formulation, the arcs in $A_w$ are subdivided into four groups, pedestrian arcs $A_w^P$, boarding arcs $A_w^b$, alighting arcs $A_w^a$ and bus arcs $A_w^B$, so that $A_w=\{A_w^P\cup A_w^b\cup A_w^a\cup A_w^B$}. A given arc $(i,j)$ in $G_w$ could belong to the strategy subgraph of a different OD pair $w^{\prime }\ne w$; then, when constructing the set of strategy subgraphs, we also define the sets $\mathsf{\Omega }(i,j)=\{w\in \mathcal{W}:(i,j)\in A_w\}$. Figure 2 represents a simplified strategy subgraph for a generic OD pair w. Note that (a) a transfer between lines at certain nodes implies a movement from the first line to the corresponding pedestrian node in the pedestrian layer using an alighting arc with no cost, plus a second movement from the pedestrian node to the destination line whose real cost will depend on the final line frequency value. (b) On the trip from $w_o$ to $w_d$, a passenger may use a succession of pedestrian arcs, regardless of their position along the trip path, or even more than one disjoint pedestrian sub-path if required. (c) In the trip from $w_o$ to $w_d$, the passenger may use a succession of bus arcs using different lines combined with transfer alighting and boarding arcs and pedestrian arcs if required. (d) The union of subgraphs $G_w,\forall w\in \mathcal{W}$ determines a graph $G_{\mathcal{W}}$ included in $G_{ML}$ that will be used to formulate the bimodal multilayer line planning problem.

6. The Line Planning Formulations

In this section, we propose two formulations for the bimodal line planning problem. The first one is an arc-based formulation intended to minimize the total travel time while determining the best line configuration and the frequency of lines subject to operation constraints in terms of minimizing the frequency of lines for cost saving purposes. The second one is a path-based formulation that takes advantage of the specific path sequence in each strategy subgraph, which additionally aims at explicitly restricting the number of transfers allowed for each OD pair in the network. In both cases, the models will determine the best line configuration, the lines’ frequency, as well as the passenger assignment to the network. To facilitate the understanding of the formulations,

Table 2. Notation used for the construction of the bimodal multilayer network.

Symbol	Description	Symbol	Description
$G_P$	The graph defining the pedestrian network, $G_P=\{N_P,A_P\}$	$A_{ML}$	Set of arcs in the bimodal multilayer directed graph
$N_P$	The set of nodes in $G_P$	$k_L$	Number of shortest paths used to compute the pool of candidate lines
$A_P$	The set of directed arcs in $G_P$	$\widehat{L}$	The intermediate line pool before pruning candidate lines
$G_B$	The graph defining the bus network, $G_B=\{N_B,A_B\}$	$L_{min}$	Minimum length of lines in the pool of candidate lines
$N_B$	The set of nodes in $G_B$	$L_{max}$	Maximum length of lines in the pool of candidate lines
$A_B$	The set of directed arcs in $G_B$	$\mathcal{D}$	The OD demand matrix
$i,j$	Generic nodes in the multilayer network	$d_{ij}$	The demand between nodes i and j in $N_P$
L	Line pool of candidate lines	$\mathcal{W}$	The set of OD pairs between nodes of $N_P$
ℓ	Generic line in L	w	A generic OD pair
$N_{\mathcal{l}}$	Nodes belonging to line ℓ in L	$w_o$	The origin node of OD pair w
$A_{\mathcal{l}}$	Directed arcs belonging to line ℓ in L	$w_d$	The destination node of OD pair w
${\varphi }_{\mathcal{l}}\left(i\right)$	Function that returns the code of node i in $N_B$ corresponding to line ℓ	$G_w$	The strategy subgraph corresponding to the OD pair w
${\varphi }^{-1}$	Operator that returns the original code in $N_B$ of a node in line ℓ	$k_P$	The number of shortest paths used to compute the strategy subgraph of OD pair w
${\mathcal{N}}_{\mathcal{l}}$	The set of coded nodes belonging to the layer of line ℓ	$v_P$	The pedestrian speed
${\mathcal{A}}_{\mathcal{l}}$	The set of coded arcs belonging to the layer of line ℓ	$v_B$	The bus speed
$\psi \left(i\right)$	The set of counterpart nodes in the bus superlayer for the node $i\in N_B$	$N_w$	The set of nodes in the strategy subgraph of OD pair w
$\psi (i,j)$	The set of counterpart arcs in the bus superlayer bus for the arc $(i,j)\in A_B$	$A_w$	The set of directed arcs in the strategy subgraph of OD pair w
${\mathcal{T}}_b\left(i\right)$	The set of boarding arcs connecting node $i\in N_B$ with its counterpart nodes in the bus superlayer	$A_w^P$	The subset of pedestrian arcs in $A_w$
${\mathcal{T}}_a\left(i\right)$	The set of alighting arcs connecting nodes in the bus superlayer with its corresponding node $i\in N_B$	$A_w^b$	The subset of boarding arcs in $A_w$
${\mathcal{T}}_b$	All the boarding arcs between nodes in $N_B$ and its counterpart nodes in the bus superlayer	$A_w^a$	The subset of alighting arcs in $A_w$
${\mathcal{T}}_a$	All the alighting arcs connecting nodes in the bus superlayer with the corresponding nodes $i\in N_B$	$A_w^B$	The subset of bus arcs in $A_w$
$\mathcal{L}\left(i\right)$	The set of lines traversing the counterpart nodes of node $i\in N_B$ in the bus superlayer	$G_{\mathcal{W}}$	The union of the strategy subgraph for all the OD pairs in $\mathcal{W}$
${\mathcal{L}}^{-1}$	Operator that returns the line traversing a coded bus node belonging to ${\mathcal{N}}_{\mathcal{l}}$	$\mathsf{\Omega }(i,j)$	Set of OD pairs w in $\mathcal{W}$ traversing arc $(i,j)\in {\cup }_{w\in \mathcal{W}}{A}_w$
$G_{ML}$	Bimodal multilayer directed graph	${\delta }^{+}\left(i\right)$	Set of successor nodes of node i in $G_w$
$N_{ML}$	Set of arcs in the bimodal multilayer directed graph	${\delta }^{-}\left(i\right)$	Set of predecessor nodes of node i in $G_w$

summarizes the notation used so far.

6.1. The Arc-Based Formulation

To construct the arc-based formulation, we consider the following set of variables.

• $x_{ij}^w$ Positive real variables representing the flow on arc $(i,j)\in A_w$ corresponding to the OD pair w.
• $y_{\mathcal{l}}$ Binary variable taking value 1 if the candidate line $\mathcal{l}\in L$ is activated in the line concept, or 0 otherwise.
• $f_{\mathcal{l}}$ Integer variable representing the frequency of the candidate line $\mathcal{l}\in L$ measured in number of buses/hour.
• $h_{\mathcal{l}}$ Integer variable representing the headway of the candidate line $\mathcal{l}\in L$ measured in minutes.

6.1.1. Objective Function

Different objective functions may be formulated according to the review of the literature presented in Section 2. The most common is the minimization of the total passenger travel time. Usually, service operators are also interested in operating as few lines as possible to reduce the fleet and operation costs, or even operating the lines with the lowest possible frequency, since the operation costs are proportional to the lines’ frequencies. Low frequencies also mean lower fleet sizes, reducing the necessary investment in vehicles. We propose here a combined weighted objective function with the goal of minimizing the total passenger travel time, the number of lines and the frequency of lines.

$$\mathbf{Min}\sum _{w\in \mathcal{W}}\sum _{(i,j)\in A_w}{t}_{ij}{x}_{ij}^w+\alpha \sum _{\mathcal{l}\in L}{y}_{\mathcal{l}}+\beta \sum _{\mathcal{l}\in L}{f}_{\mathcal{l}},$$

(1)

being $t_{ij}$ the travel time over the arc $(i,j)$ for the different types of arcs in $G_w$, and $\alpha$ and $\beta$ are parameters to weight the importance of the subsidiary terms in the objective function.

Since the speed on pedestrian arcs is low compared to the speed of the bus, passengers will take as many bus lines as possible to reach their destinations. A good line configuration should increase the number of direct trips and minimize the use of pedestrian arcs. As previously reported, if, for instance, the number of lines is severely restricted, passengers always have the possibility of carrying out their trips by using the shortest pedestrian path connecting their origin and destination. The minimization of passenger flows over arcs in the pedestrian path of every OD pair w together with a relaxation in the constraints limiting the number of lines will guarantee the maximum demand coverage result.

6.1.2. Constraints

The first sets of constraints ensures that the demand of each OD pair is fulfilled; then, for each OD pair w, the number of passengers starting their trips at $w_o$ must be precisely $d_{w_o{w}_d}$, which must also be the number of passengers arriving at the destination $w_d$. To facilitate the formulation process, we define the sets ${\delta }^{+}\left(i\right)$ and ${\delta }^{-}\left(i\right)$ to, respectively, denote the set of successor and predecessor nodes of node i in each strategy subgraph $G_w$.

$$\sum _{j\in {\delta }^{+}\left(w_o\right)}{x}_{w_{o}j}^w=d_{w_o{w}_d}w\in \mathcal{W}$$

(2)

$$\sum _{j\in {\delta }^{-}\left(w_d\right)}{x}_{j{w}_d}^w=d_{w_o{w}_d}w\in \mathcal{W}$$

(3)

The following set of constraints enforces a balance for all nodes i belonging to $G_w$.

$$\sum _{j\in {\delta }^{+}\left(i\right)}{x}_{ij}^w-\sum _{j\in {\delta }^{-}\left(i\right)}{x}_{ji}^w=0i\in N_w\setminus \{w_o,w_d\},w\in \mathcal{W}$$

(4)

For each candidate line, its frequency and headway are related by the following equation (note that the frequency is measured in buses/hour and the headway in min). Furthermore, the frequency must be 0 if the candidate line is not activated and should also be bounded by two reasonable values $f_{min}$ and $f_{max}$. The first is related to the quality of service perceived by users, who probably will expect a service at least every 20 or 30 min. The second one is usually set by the operator depending on the available fleet and the observed demand. Reasonable values could vary between 20 and 30 buses/hour.

$$f_{\mathcal{l}}·h_{\mathcal{l}}=60·y_{\mathcal{l}}\mathcal{l}\in L$$

(5)

$$f_{min}·y_{\mathcal{l}}\le f_{\mathcal{l}}\le f_{max}·y_{\mathcal{l}}\mathcal{l}\in L$$

(6)

Note that if $y_{\mathcal{l}}=0$, the variable $f_{\mathcal{l}}$ is also 0 and the constraints (5) hold, even if the value of $h_{\mathcal{l}}$ is not determined. This is irrelevant, since line ℓ will not be selected in the final line plan. The set of constraints (5) could be linearized by including new binary variables ${\lambda }_{\mathcal{l}}^{\sigma }$ taking the value 1 if the frequency of the line ℓ takes a value $\sigma$ from a subset of feasible values—for instance, $\mathsf{\Lambda }=\{2,3,4,5,6,10,12,15,20\}$, the set of integer divisors of 60. Then,

$$f_{\mathcal{l}}=\sum _{\sigma \in \mathsf{\Lambda }}\sigma {\lambda }_{\mathcal{l}}^{\sigma }\mathcal{l}\in L$$

$$h_{\mathcal{l}}=\sum _{\sigma \in \mathsf{\Lambda }}\frac{60}{\sigma }{\lambda }_{\mathcal{l}}^{\sigma }\mathcal{l}\in L$$

$$\sum _{\sigma \in \mathsf{\Lambda }}{\lambda }_{\mathcal{l}}^{\sigma }=1\mathcal{l}\in L$$

The linearization process is completed by removing the constraint set (5).

The service operator could be interested in fixing a maximum number of lines for the network; for instance, the line planning could be simply an update of the existing line planning, and the service operator could fix the maximum number of lines to maintain approximately the same operating cost. In this case, the following set of constraints should be included:

$$\sum _{\mathcal{l}\in L}{y}_{\mathcal{l}}\le M,$$

(7)

being M the maximum allowed number of lines. Otherwise, M can be set to a large value or simply the constraint set (7) can be removed.

If the flow on a bus arc is positive, the corresponding line traversing this arc must be active. In this case, the frequency of this line must be sufficient to accommodate the total flow on the link. Recall that each line in the candidate line pool is stored in a unique layer and, as a consequence, each arc in this layer corresponds to only one line, which can be accessed using the operator ${\mathcal{L}}^{-1}$. To formulate this constraint, it is necessary to incorporate a parameter that defines the bus capacity, which will be denoted as $Cap$.

$$\sum _{w^{\prime }\in \mathsf{\Omega }(i,j)}{x}_{ij}^{w^{\prime }}\le Cap·f_{{\mathcal{L}}^{-1}(i,j)}(i,j)\in A_w^B,w\in \mathcal{W}$$

(8)

The following set of constraints specifies the travel time $t_{ij}$ in the boarding arcs as a function of the headway of the destination lines.

$$t_{ij}=\frac{h_{{\mathcal{L}}^{-1}\left(j\right)}}{2}(i,j)\in {\mathcal{T}}_w^b,w\in \mathcal{W}$$

(9)

Note that the travel times for the rest of the arcs in $G_{\mathcal{W}}$ are constant:

• $t_{ij}=0\forall (i,j)\in A_w^a$,
• $t_{ij}=l_{ij}/v_B\forall (i,j)\in A_w^B$,
• $t_{ij}=l_{ij}/v_P\forall (i,j)\in A_w^P$,

being $l_{ij}$ the length of link $(i,j)$.

Although, in the case of a transit system, the total number of buses of different lines traversing a certain arc could be quite sufficient, at some specific locations (for instance, a street with only a lane per direction or in general links that for any reason can suffer from congestion), it could be necessary to bound the total number of vehicles traversing the link during the planning horizon. If such a case is necessary, for specific arcs $(i,j)$ in a certain subset $\widehat{A}\in A_B$, the following constraints could be used:

$$\sum _{(i^{\prime },j^{\prime })\in \psi (i,j)}{f}_{{\mathcal{L}}^{-1}(i^{\prime },j^{\prime })}\le \mathcal{F}(i,j)\in \widehat{A},$$

(10)

being $\mathcal{F}$ a global limit to the sum of the affected frequencies, measured in buses/hour.

Sometimes, the service operator may be also interested in guaranteeing that a subset of nodes in $N_B$ act as interchange stations in the network, e.g., a node near a faculty, a hospital or a commercial or leisure center. Suppose that there is a subset of nodes ${\widehat{N}}_B\subset N_B$ with given minimum requirements regarding the number of lines that must allow direct access. Let ${\pi }_i$ be the minimum number of lines required at node $i\in {\widehat{N}}_B$. Then, the following set of constraints must be added to the model.

$$\sum _{\mathcal{l}\in \mathcal{L}\left(i\right)}{y}_{\mathcal{l}}\ge {\pi }_{i}i\in {\widehat{N}}_B$$

(11)

Additionally, the service operator could be also interested in limiting the maximum number of vehicles required to operate the final line concept. In this case, the following constraint should be imposed.

$$\sum _{\mathcal{l}\in L}\frac{2}{v_B}\sum _{(i,j)\in A_{\mathcal{l}}}{l}_{(i,j)}\frac{f_{\mathcal{l}}}{60}\le FS,$$

(12)

where $l_{(i,j)}$ represents the length of link $(i,j)\in A_{\mathcal{l}}$ and $FS$ is the maximum allowed number of vehicles. The left-hand side of (12) is the summation of the number of vehicles needed to operate each line. For a specific line ℓ, the required fleet size is obtained by dividing the cycle time (the amount of time needed to complete a round trip along the line), ${\sum }_{(i,j)\in A_{\mathcal{l}}}\frac{2l_{(i,j)}}{v_B}$, by the headway of the line, in this case represented in terms of the frequency as $\frac{60}{f_{\mathcal{l}}}$.

In summary, the arc-based formulation consists of minimizing (1) subject to (2)–(9), and optionally to (10)–(12), and fulfilling

$$\begin{array}{ccc}\hfill & x_{i,j}^w\ge 0,w\in \mathcal{W},(i,j)\in {\cup }_{w\in \mathcal{W}}{A}_w& \\ & y_{\mathcal{l}}\in \{0,1\},\mathcal{l}\in L& \\ & f_{\mathcal{l}},h_{\mathcal{l}}\ge 0\mathrm{integers},\mathcal{l}\in L& \\ & t_{ij}\ge 0,(i,j)\in A_w^b\end{array}$$

6.1.3. Comments Regarding the Use of a Multilayer-Type Structure

The multilayer structure used to formulate the problem drastically simplifies the formulations of constraint sets (8) and (9) (as well as the optional sets (10) and (12)) when compared to a classical single-layer graph. For instance, in a single-layer model, the travel time of an incoming arc at a transfer node must be formulated while considering the frequencies of all the possible lines traversing this node, which results in a non-linear constraint of type

$$t_{ij}=\frac{1}{2}\frac{60}{{\sum }_{\mathcal{l}\in \mathsf{\Delta }\left(j\right)}{f}_{\mathcal{l}}},$$

where, for simplicity, the set of lines traversing node j is supposed to belong to a certain set $\mathsf{\Delta }\left(j\right)$. Even when using headway variables instead of frequencies, the constraint remains non-linear.

In the classical formulation of constraints (9), since each arc supports many possible lines, the use of arc–line binary variables $y_{i,j}^{\mathcal{l}}$ is required. This implies the inclusion of an important number of binary variables, at least $L\times \left({\cup }_{w\in \mathcal{W}}{A}_w\right)$, and a similar number of constraints relating $y_{i,j}^{\mathcal{l}}$ and $y_{\mathcal{l}}$.

At the cost of building a huge multilayer graph structure, which would require more memory and the more complex preparation of auxiliary sets, the resultant line planning models are simpler than those obtained using a simple graph and, in most cases, can be solved in relatively short computation times, as reported later in the computational section.

6.2. The Path-Based Formulation

Starting from the definition of the strategy subgraphs $G_w,w\in \mathcal{W}$, it is possible to work with specific path variables instead of using variables $x_{ij}^w$, reducing, in general, even further the number of variables in the line planning formulation. In this way, for each subgraph $G_w$, we can compute all the paths $\rho$ connecting $w_o$ and $w_d$. Let ${\mathsf{\Gamma }}_w$ be the set of paths between the origin and the destination of the OD pair w. Thus, $x_w^{\rho }$ measures the number of passengers following path $\rho$ between $w_o$ and $w_d$ for each $w\in \mathcal{W}$. Note that ${\mathsf{\Gamma }}_w\ge k_{\mathcal{W}}+1$, the number of paths used to calculate $G_w$ in Section 5, since, as explained, several other trip combinations emerge when using subpath combinations of the $k_{\mathcal{W}}+1$ paths.

Let $G_{{\mathsf{\Gamma }}_w}=\{N_{{\mathsf{\Gamma }}_w},A_{{\mathsf{\Gamma }}_w}\}$ be the resulting graph that includes all nodes and arcs in paths $\rho \in {\mathsf{\Gamma }}_w$. As in the arc-based formulation, it is convenient to divide the arc set into four subsets, $A_{{\mathsf{\Gamma }}_w}=\{A_{{\mathsf{\Gamma }}_w}^P\cup A_{{\mathsf{\Gamma }}_w}^B\cup A_{{\mathsf{\Gamma }}_w}^b\cup A_{{\mathsf{\Gamma }}_w}^a\}$, containing, respectively, pedestrian, bus, boarding and alighting arcs belonging to $G_{{\mathsf{\Gamma }}_w}$. Let $N_w^{\rho }$ and $A_w^{\rho }$ be the set of nodes and arcs defining a path $\rho$ in ${\mathsf{\Gamma }}_w$, respectively. As in the arc-based formulation, a specific arc $(i,j)\in A_{{\mathsf{\Gamma }}_w}$ can also be included in the arc sets $A_{{\mathsf{\Gamma }}_w^{\prime }}$ of other pairs $w^{\prime }\ne w$; thus, it is convenient to define the set of paths belonging to all the pairs traversing a given arc, ${\mathsf{\Theta }}_w(i,j)=\{\rho \in {\mathsf{\Gamma }}_w:(i,j)\in A_{{\mathsf{\Gamma }}_w}\}$.

6.2.1. Basic Constraints

Note that the constraints (5)–(7) and optional constraints (10) and (12) of the arc-based formulation are not affected when changing the arc flow variables to path-based flow variables, and they consequently remain unchanged. The demand constraints (2)–(4) of the arc-based formulation are now replaced by

$$\sum _{\rho \in {\mathsf{\Gamma }}_w}{x}_w^{\rho }=d_{w_o{w}_d}w\in \mathcal{W},$$

(13)

and the coupling constraints (8) of the arc-based formulation are now written as

$$\sum _{w^{\prime }\in \mathcal{W}}\sum _{{\rho }^{\prime }\in {\mathsf{\Theta }}_{w^{\prime }}(i,j)}{x}_{w^{\prime }}^{{\rho }^{\prime }}\le Cap·f_{{\mathcal{L}}^{-1}(i,j)}(i,j)\in A_w^{\rho }\cap A_{{\mathsf{\Gamma }}_w}^B,\rho \in {\mathsf{\Gamma }}_w,w\in \mathcal{W},$$

(14)

where, as usual, ${\mathcal{L}}^{-1}(i,j)$ returns the unique line traversing arc $(i,j)$ in the bus superlayer. Finally, the set of constraints (9) that define the travel times at boarding arcs in the arc-based formulation are replaced by

$$t_{ij}=\frac{h_{{\mathcal{L}}^{-1}\left(j\right)}}{2}(i,j)\in A_w^{\rho }\cap A_{{\mathsf{\Gamma }}_w}^b,\rho \in {\mathsf{\Gamma }}_w,w\in \mathcal{W},$$

(15)

while the travel times for the rest of the arcs in ${\cup }_{w\in \mathcal{W}}{G}_{{\mathsf{\Gamma }}_w}$ are constant:

• $t_{ij}=0,\forall (i,j)\in A_w^{\rho }\cap A_{{\mathsf{\Gamma }}_w}^a,\rho \in {\mathsf{\Gamma }}_w,w\in \mathcal{W}$,
• $t_{ij}=l_{ij}/v_B,\forall (i,j)\in A_w^{\rho }\cap A_{{\mathsf{\Gamma }}_w}^B,\rho \in {\mathsf{\Gamma }}_w,w\in \mathcal{W}$,
• $t_{ij}=l_{ij}/v_P,\forall (i,j)\in A_w^{\rho }\cap A_{{\mathsf{\Gamma }}_w}^P,\rho \in {\mathsf{\Gamma }}_w,w\in \mathcal{W}$,

being $l_{ij}$ the length of link $(i,j)$.

6.2.2. Objective Function and Basic Model

Now, the arc-based objective function (1) can be rewritten as

$$\mathbf{Min}\sum _{w\in \mathcal{W}}\sum _{\rho \in {\mathsf{\Gamma }}_w}{x}_w^{\rho }\sum _{(i,j)\in A_{{\mathsf{\Gamma }}_w}}{t}_{ij}+\alpha \sum _{\mathcal{l}\in L}{y}_{\mathcal{l}}+\beta \sum _{\mathcal{l}\in L}{f}_{\mathcal{l}}$$

(16)

and the full path-based formulation consists of minimizing (16) subject to (5)–(7), (13)–(15) and, optionally, (10)–(12), and satisfying the following domain constraints:

$$\begin{array}{ccc}& x_w^{\rho }\ge 0,\rho \in {\mathsf{\Gamma }}_w,w\in \mathcal{W}& \\ & y_{\mathcal{l}}\in \{0,1\},\mathcal{l}\in L& \\ & f_{\mathcal{l}},h_{\mathcal{l}}\ge 0\mathrm{integers},\mathcal{l}\in L& \\ & t_{ij}\ge 0,(i,j)\in A_w^{\rho }\cap A_{{\mathsf{\Gamma }}_w}^b\end{array}$$

6.2.3. Path-Based Formulation: Additional Features

Taking advantage of the path structure of this formulation, the paths $\rho$ in ${\mathsf{\Gamma }}_w$ can be classified according to the number of transfers, denoted as $n_{\rho }$, that passengers experiment with when traversing them. In this way, let ${\xi }_w^1=\{\rho \in {\mathsf{\Gamma }}_w:n_{\rho }=1\}$ be the subset of paths with 1 transfer in ${\mathsf{\Gamma }}_w$, and ${\xi }_w^2=\{\rho \in {\mathsf{\Gamma }}_w:n_{\rho }=2\}$, ${\xi }_w^3=\{\rho \in {\mathsf{\Gamma }}_w:n_{\rho }\ge 2\}$ the subsets of paths with 2 and more than 2 transfers. Now, it is possible to define new terms in the model objective function to specifically minimize the flow traversing paths in ${\xi }_w^1$, ${\xi }_w^2$ and ${\xi }_w^3$, weighted with its corresponding importance:

$$\begin{array}{ccc}\hfill & \mathbf{Min}\sum _{w\in \mathcal{W}}\sum _{\rho \in {\mathsf{\Gamma }}_w}{x}_w^{\rho }\sum _{(i,j)\in A_{{\mathsf{\Gamma }}_w}}{t}_{ij}+\alpha \sum _{\mathcal{l}\in L}{y}_{\mathcal{l}}+\beta \sum _{\mathcal{l}\in L}{f}_{\mathcal{l}}+\hfill & \hfill \\ \\ \hfill & {\gamma }_1\sum _{w\in \mathcal{W}}\sum _{\rho \in {\xi }_w^1}{x}_w^{\rho }+{\gamma }_2\sum _{w\in \mathcal{W}}\sum _{\rho \in {\xi }_w^2}{x}_w^{\rho }+{\gamma }_3\sum _{w\in \mathcal{W}}\sum _{\rho \in {\xi }_w^3}{x}_w^{\rho }\hfill & \hfill \end{array}$$

(17)

The path-based formulation also allows us to work with average travel times instead of using the total travel time considered in the arc-based formulation. In several cases, the total waiting time could hide some poor results for some OD pairs for the benefit of others. Moreover, the average travel time may be included as a constraint, instead of being incorporated into the objective function, although, in this situation, it is possible to reach infeasible situations, which can be avoided by penalizing any excess over the specified target value. Denoting as $t_{\rho }^w$ the time needed to traverse the path $\rho$ in ${\mathsf{\Gamma }}_w$, given by,

$$t_{\rho }^w=\sum _{(i,j)\in A_{{\mathsf{\Gamma }}_w}^P}\frac{l_{ij}}{v_P}+\sum _{(i,j)\in A_{{\mathsf{\Gamma }}_w}^B}\frac{l_{ij}}{v_B}+\sum _{(i,j)\in A_{{\mathsf{\Gamma }}_w}^b}\frac{1}{2}{h}_{{\mathcal{L}}^{-1}(i,j)},$$

(18)

the average travel time experienced by users traveling from $w_o$ to $w_d$ is

$$t^w=\frac{1}{d_{w_o{w}_d}}\sum _{\rho \in {\mathsf{\Gamma }}_w}{t}_{\rho }^w{x}_{\rho }^w,w\in \mathcal{W},$$

(19)

and we can write an alternative objective function as

$$\begin{array}{ccc}\hfill & \mathbf{Min}\sum _{w\in \mathcal{W}}\frac{1}{d_{w_o{w}_d}}\sum _{\rho \in {\mathsf{\Gamma }}_w}{t}_{\rho }^w{x}_{\rho }^w+\alpha \sum _{\mathcal{l}\in L}{y}_{\mathcal{l}}+\beta \sum _{\mathcal{l}\in L}{f}_{\mathcal{l}}+\hfill & \hfill \\ \\ \hfill & {\gamma }_1\sum _{w\in \mathcal{W}}\sum _{\rho \in {\xi }_w^1}{x}_w^{\rho }+{\gamma }_2\sum _{w\in \mathcal{W}}\sum _{\rho \in {\xi }_w^2}{x}_w^{\rho }+{\gamma }_3\sum _{w\in \mathcal{W}}\sum _{\rho \in {\xi }_w^3}{x}_w^{\rho },\hfill & \hfill \end{array}$$

(20)

incorporating linear constraints (18) into the formulation.

Additionally, the designer could be interested in guaranteeing that, for each OD pair w in the line concept solution, the average travel time be as close as possible to a design value ${\tau }_w$, which can be established as the travel between $w_o$ and $w_d$ without waiting. Then, new constraints imposing the proximity to this value could be added to the model:

$$t^w\le {\tau }_w+{\chi }_w,w\in \mathcal{W},$$

where ${\chi }_w$ represents deviations with respect to the target values, which can be added to the objective function with positive weights.

7. Computational Experiments

7.1. The Mandl Network

In this section, we apply both the arc-based and the path-based formulations to the network analyzed by Mandl and present the results by varying the maximum number of allowed lines, as in other research included in the revised literature. Before starting this comparison, several aspects of the Mandl scenario must be commented on in relation to our bimodal approach.

(a)

The available information of the Mandl network is limited to the code of nodes and edges, the OD matrix and the travel time on edges. There is no information about the coordinates of the nodes or the real distances between the nodes.

(b)

No information is known about the real scale of the network. In fact, if we suppose that buses run at a constant speed, the lengths of edges in the classical representation of this network (the graph reported in many research papers) are not real. Then, it is very difficult to determine real walking times over pedestrian links to be considered as part of a bimodal approach.

(c)

Consequently, to produce a scenario that is as similar as possible to the one analyzed in the literature, we will use the information on bus travel times on links in the bus superlayer and impose large costs on pedestrian links in the pedestrian layer, trying to ensure that they are never chosen by passengers.

(d)

Since the original scenario is a unimodal network, passengers who do not use the bus (passengers who start their trips at nodes that do not belong to any line) are considered as unattended demand in the previous works. To fairly compare with them, we will consider that passengers using pedestrian links (even if their costs are high) in our bimodal network are part of the unattended demand. The pedestrian layer is set identical to the bus layer, i.e., every link in the network can be used as a pedestrian or as a bus link.

Although many researchers have based their algorithms on the Mandl network, it must be noted that the results have been presented in different formats and with different amplitudes. Therefore, only the following commonly reported indicators will be considered here for comparison purposes.

(1)

$d_0$, Percentage of demand that is served directly.

(2)

$d_1$, Percentage of demand that is served with one transfer.

(3)

$d_2$, Percentage of demand that is served with two transfers.

(4)

$d_u$, Percentage of unattended demand.

(5)

$FS$, Fleet size—usually the most significant measure from an operator’s point of view.

Surprisingly, in most of the previous metaheuristic applications, there is no information on the computational effort required to achieve good solutions or on the number of replications used to compute average results, perhaps due to the strategic nature of this problem. In addition, there are important differences in the assumptions considered by different researchers to illustrate the results of their work. For instance, some papers consider unlimited frequencies while, at the same time, fixing a constant value for the boarding and transfer times. There is also a discrepancy in several indicators concerning travel, waiting and transfer times among different works; see, for instance, the works of Bagloee and Ceder [44] or Buba and Lee [77], reporting times consistent with the structure and length of the lines in the final achieved line concepts, and those of Fan and Mumford [70] or Vlachopanagiotis et al. [80], where the authors report average in-vehicle waiting times.

Figure 3 shows the OD matrix and a representation of the desired movements on the network. As shown, the matrix is symmetric and contains a total number of 15,570 trips. As in several previous works, we consider buses with a capacity of 40 pax. and maximum headway values of 20 min.

Since we can expect a direct relationship between the quality of the solutions obtained and the number and quality of the routes in the line pool, we choose a large number $k_L=12$ of alternative paths when computing the initial pool of candidate lines, as described in Section 4. We initially select candidate routes whose length (measured in minutes) is greater than half the diameter of the network and less than 2 times the diameter. Since, as previously mentioned, there is no information about the network scale, we measure the diameter of the network by considering the total travel time corresponding to the shortest path between nodes 1 and 13. We fix $\vartheta =0.5$, removing from the intermediate line pool $\widehat{L}$ the candidate lines whose accumulated demand is below $50\%$ of the best accumulated demand. From the 876 initial candidate routes (holding length conditions), the cleaning method reduces the line pool size to 351 candidate lines.

We also consider $k_P=12$ when computing the number of alternative paths followed by passengers to travel between the origin and destination of each OD pair, i.e., when computing the strategy subgraph for each OD pair as described in Section 5. The higher the number of alternative paths, the greater the freedom to move over the network, and the greater the size of the resulting optimization model.

The preparation, including data reading, line pool generation (intermediate line pool construction and line pool cleaning), multilayer network construction (including the bus superlayer composed of as many layers as candidate lines) and the construction of the OD pair strategy subgraphs, consumes approximately 3 min of time on a personal computer with an Intel I7-1165G7 processor running at 2.80 GH with 16 GB of RAM. All these methods have been programmed in the Python programming language by developing a specialized object-oriented library to deal with multimodal transportation problems. The full multimodal network contains 2516 nodes and 9344 directed links. The number of nodes and links used to formulate the models (obtained as the union of nodes and links of the OD subgraphs) is 1322 nodes and 4064 directed links, which contributes to reducing the number of necessary variables and constraints. Note that the model preparation is performed once since the same pool of candidate lines is used to solve the arc-based and the path-based formulations, independently of the maximum allowed number of lines in the final line concept.

We solve the two basic models by using a standard branch and cut algorithm, imposing a maximum number of $4,6$ and 8 lines. To give an idea of the size and difficulty of both formulations, the arc-based formulation contains 13,377 variables (12,324 continuous, 702 integer, 351 binary) and 14,997 constraints, while the path-based formulation contains 9724 variables (8671 continuous, 702 integer, 351 binary) and 9661 constraints. We use weights $\alpha =1$ and $\beta =50$ in the objective function to penalize high frequencies and impose the maximum number of lines $(4,6,8)$ by including the appropriate constraint in the model. The arc-based formulation is solved in 28 s, and the path-based one in 34 s.

Our results in terms of comparison indicators are shown in Table 3. We also highlight in red the best previous results for different maximum numbers of lines. AIVTT and ATT represent the average in-vehicle travel time and the average travel time, respectively. As shown, the values of the demand indicators are competitive in comparison with those obtained in previous works using metaheuristic approaches. As shown, in the experiment with four lines, our approach generates a solution with a value of the attended demand that is $3\%$ worse than the best reported solution in the literature, Arbex and da Cunha [76], and ranges from 5 to 18 of the reported approaches. The score obtained for the case of six lanes is slightly better, only $1\%$ worse than the best reported solution—see Zhao and Zeng [69]—achieving third place among the works analyzed. In the case of eight lines, our model achieves again fifth position according to the percentage of directly served demand. In this case, the best reported solution corresponds to the work of Nayeem et al. [73]. As expected, from a user perspective, the general results improve as the number of lines increases. Recall that, unlike the metaheuristic approaches, in the mathematical formulations, the number of lines is not fixed but bounded with a maximum value. In each of the performed experiments, the best solution in terms of the percentage of demand that is served directly is attained for a number of lines equal to the maximum allowed value. This result seems to be logical, since a higher number of lines covers more nodes in the network, thus increasing the possibilities of taking a bus. However, this behavior is not general, as can be seen, for instance, in the case of four and six lines for the works of Chakroborty and Wivedi [39], Zhao and Zeng [68] and Fan and Mumford [70], among others, where more lines give rise to a lower value of $d_0$. The reason is clear: the quality of the solution depends on the number of lines but also on the shape of the lines and the nodes that each line visits. A line traversing nodes with low demand and few connections contributes poorly to the movement of passengers, and the addition of a new line sharing all its stops with nodes belonging to the existing lines does not improve the total captured demand but, contrarily, increases the operating costs.

On the one hand, the results in Table 3 confirm the effectiveness of several previously used metaheuristics, such as, for instance, the ones proposed by Nayeem et al. [73] and Buba and Lee [77], which were not compared with the exact solution procedures. On the other hand, our approach has generated good-quality solutions even when working with a fixed pool of candidate lines, which possibly could be improved by using a more sophisticated line pool generation method, such as the one proposed by Mauttone and Urquhart [19] or Gattermann et al. [55]. In any case, as the latter authors noted, the results are highly sensitive to the size of the pool and the quality of the lines contained therein, which, in our case, required the realization of several tests on the parameters $k_L$ (the number of shortest paths used in the generation of lines for each pair of nodes), $k_P$ (the number of shortest paths used in the generation of the strategy subgraph for each OD pair) and $\vartheta$ (the percentage of lines to be removed from the intermediate line pool). In this process, we also corroborated the experiences of Kechagiopoulos and Beligiannis [86] and Buba and Lee [77] regarding the effect of the length of the candidate routes on the quality of the solutions obtained. It is worth mentioning that our interest with this experiment is only in validating its applicability using a benchmark instance to assess the methodology’s ability to address, in the future, more complicated problems where several modes are simultaneously operating on the same network. In these complicated scenarios, we will likely need to consider a hybridized approach by combining the mathematical models with metaheuristics algorithms or applying the metaheuristics in a second phase starting from one or a set of good achieved solutions, depending on the nature of the algorithm.

Table 3. Comparison of results with respect to previous works.

	4 Lines						6 Lines						8 Lines
	do	d1	d2	du	AIVTT (ATT)	FS	do	d1	d2	du	AIVTT (ATT)	FS	do	d1	d2	du	AIVTT (ATT)	FS
Mandl (1980) [26]	69.94	29.03	0.13	0	12.9
Baaj and Mah. (1991) [36]							78.61	21.39	0	0	11.86	89	79.96	20.04	0	0	11.86	77
Kidwai (1998) [87]	72.95	26.91	0.13	0	12.72		77.92	19.62	2.4	0	11.87		84.73	15.27	0	0	11.22
Chakroborty and Dwivedi (2002) [39]	86.86	12	1.14	0	11.9		86.04	13.96	0	0	10.3		90.38	9.68	0	0	10.46
Zhao et al. (2006) [68]	95.31	4.69	0	0	(11.89)	99	95.18	4.82	0	0	(12.26)	89	95.44	4.56	0	0	(12.55)	77
Zhao et al. (2008) [69]	96.66	3.34	0	0	(11.68)	99	98.39	1.61	0	0	(11.7)	89	95.83	4.17	0	0	(12.07)	77
Fan and Machemehl (2008) [88]	93.26	6.74	0	0	11.37		91.52	8.48	0	0	10.48		94.54	5.46	0	0	10.36
Fan and Mundford (2010) HC Avg. [70]	91.83	8.17	0	0	11.69		90.23	9.26	0.51	0	10.78		93.23	6.18	0.59	0	10.69
Fan and Mundford(2010) SA Avg. [70]	92.48	7.52	0.51	0	11.55		90.87	8.74	0.39	0	10.65		93.65	5.88	0.47	0	10.58
Nikolic and Teodorovic (2013) [71]	92.1	7.19	0.71	0	10.51		95.63	4.37	0	0	10.23		98.97	1.03	0	0	10.09
Nayeem et al. (GAWIP) (2014) Avg. [73]	93.76	5.34	0.9	0	10.45		98.08	1.92	0	0	10.14		99.54	0.46	0	0	10.05
Nayeem et al. (GAWE) (2014) Avg. [73]	93.39	5.55	1.06	0	10.5		97.5	2.49	0.01	0	10.17		99.28	0.72	0	0	10.07
Nikolic andTeodorovic (User) (2014) [74]	95.05	4.95	0	0	10.36 (11.96)	94	94.34	5.65	0	0	10.21 (11.86)	99	96.4	3.6	0	0	10.15 (11.91)	99
Kechagiopoulos and Beligiannis (2014) [86]	90.52	8.75	0.73	0	10.71		95.62	4.28	0.1	0	10.28		97.47	2.53	0	0	10.17
Arbex and da Cunha (2015) [76]	98.27	1.73	0	0	11.13 (14.35)	79	98.2	1.8	0	0	11.55 (13.86)	77	98.95	1.35	0	0	11.24 (13.72)	74
Zhao et al. (2015) [75]	92.95	7.05	0	0	(13.39)
Buba and Lee (2018) Avg. [77]	90.43	9.57	0	0	11.39 (13.63)	98	95.65	4.35	0	0	10.79 (12.49)	95	95.74	4.24	0	0	10.70 (12.94)	98
Katsaragakis et al. (2020) Avg. [79]	89.60	9.985	0.414	0	10.67		95.25	4.012	0.038	0	10.26		98.47	1.523	0	0	10.125
Vlachopanagiotis et al. (2021) [80]	95.5	4.5	0	0	10.9 (14.54)	78	96.5	3.5	0	0	10.5 (13.62)	77	99.1	0.9	0	0	10.3 (12.95)	77
This research	95.12	4.88	0	0	10.5 (12.96)		97.39	2.61	0	0	10.15 (12.5)		98.65	1.35	0	0	10.1 (12.48)

Table 3. Comparison of results with respect to previous works.

	4 Lines						6 Lines						8 Lines
	do	d1	d2	du	AIVTT (ATT)	FS	do	d1	d2	du	AIVTT (ATT)	FS	do	d1	d2	du	AIVTT (ATT)	FS
Mandl (1980) [26]	69.94	29.03	0.13	0	12.9
Baaj and Mah. (1991) [36]							78.61	21.39	0	0	11.86	89	79.96	20.04	0	0	11.86	77
Kidwai (1998) [87]	72.95	26.91	0.13	0	12.72		77.92	19.62	2.4	0	11.87		84.73	15.27	0	0	11.22
Chakroborty and Dwivedi (2002) [39]	86.86	12	1.14	0	11.9		86.04	13.96	0	0	10.3		90.38	9.68	0	0	10.46
Zhao et al. (2006) [68]	95.31	4.69	0	0	(11.89)	99	95.18	4.82	0	0	(12.26)	89	95.44	4.56	0	0	(12.55)	77
Zhao et al. (2008) [69]	96.66	3.34	0	0	(11.68)	99	98.39	1.61	0	0	(11.7)	89	95.83	4.17	0	0	(12.07)	77
Fan and Machemehl (2008) [88]	93.26	6.74	0	0	11.37		91.52	8.48	0	0	10.48		94.54	5.46	0	0	10.36
Fan and Mundford (2010) HC Avg. [70]	91.83	8.17	0	0	11.69		90.23	9.26	0.51	0	10.78		93.23	6.18	0.59	0	10.69
Fan and Mundford(2010) SA Avg. [70]	92.48	7.52	0.51	0	11.55		90.87	8.74	0.39	0	10.65		93.65	5.88	0.47	0	10.58
Nikolic and Teodorovic (2013) [71]	92.1	7.19	0.71	0	10.51		95.63	4.37	0	0	10.23		98.97	1.03	0	0	10.09
Nayeem et al. (GAWIP) (2014) Avg. [73]	93.76	5.34	0.9	0	10.45		98.08	1.92	0	0	10.14		99.54	0.46	0	0	10.05
Nayeem et al. (GAWE) (2014) Avg. [73]	93.39	5.55	1.06	0	10.5		97.5	2.49	0.01	0	10.17		99.28	0.72	0	0	10.07
Nikolic andTeodorovic (User) (2014) [74]	95.05	4.95	0	0	10.36 (11.96)	94	94.34	5.65	0	0	10.21 (11.86)	99	96.4	3.6	0	0	10.15 (11.91)	99
Kechagiopoulos and Beligiannis (2014) [86]	90.52	8.75	0.73	0	10.71		95.62	4.28	0.1	0	10.28		97.47	2.53	0	0	10.17
Arbex and da Cunha (2015) [76]	98.27	1.73	0	0	11.13 (14.35)	79	98.2	1.8	0	0	11.55 (13.86)	77	98.95	1.35	0	0	11.24 (13.72)	74
Zhao et al. (2015) [75]	92.95	7.05	0	0	(13.39)
Buba and Lee (2018) Avg. [77]	90.43	9.57	0	0	11.39 (13.63)	98	95.65	4.35	0	0	10.79 (12.49)	95	95.74	4.24	0	0	10.70 (12.94)	98
Katsaragakis et al. (2020) Avg. [79]	89.60	9.985	0.414	0	10.67		95.25	4.012	0.038	0	10.26		98.47	1.523	0	0	10.125
Vlachopanagiotis et al. (2021) [80]	95.5	4.5	0	0	10.9 (14.54)	78	96.5	3.5	0	0	10.5 (13.62)	77	99.1	0.9	0	0	10.3 (12.95)	77
This research	95.12	4.88	0	0	10.5 (12.96)		97.39	2.61	0	0	10.15 (12.5)		98.65	1.35	0	0	10.1 (12.48)

7.2. A Larger Case Illustration

In this section, with the aim of examining the behavior of the multilayer formulation in a larger scenario, we will use the main network of Seville. Figure 4 shows the bus layer, containing 141 nodes, 454 directed links and 19,440 OD pairs. To test the formulation in adverse conditions, we consider a populated OD matrix with 107,780 trips, a number slightly higher than the real peak-hour bus demand in the city, which have been intentionally assigned in a random fashion to the OD pairs, thus giving rise to an unfavorable “non-structured” set of movements, which clearly complicates the determination of the line plan and which is depicted, according to the origin of the demand, in four graphs in Figure 5. For simplicity, as in the Mandl network case, we will consider a pedestrian layer identical to the graph supporting the bus mode.

To generate the lines, we consider minimum and maximum feasible length values of 4 and 10 km, respectively. These values are similar to those used by the transportation company currently operating the bus service. Moreover, we impose a maximum number of 45 lines, as in the real case. In this illustration, we reduce both $k_L$ and $k_P$ to a value of 8 and maintain $\vartheta =0.5$ of the maximum accumulated demand. The intermediate line pool contains 9867 candidate lines. After removing the less promising ones, the final pool of candidate lines is composed of 1638 lines. The multilayer network contains 31,342 nodes and 121,986 links. The arc-based model contains 1638 binary variables, 3276 integer variables, 1,309,781 variables measuring the passenger flow over links, 1,182,238 balance constraints and 46,794 constraints corresponding to transfer links. We first solve the model to optimality in the same processor as in the Mandl experiment, reaching a computation time of 3400 s. We have repeated the experiment using an Intel I9-12900 processor running at 3.5 GHz with 64 GB of RAM; in this case, the computation time decreases substantially to 813 s. The general results obtained for the Seville network are as follows (Figure 6 shows the optimal line plan).

• The set of lines practically covers the totality of nodes. This is a good indicator of the accessibility of the bus network, since practically all destinations are reachable by bus.
• The few uncovered nodes are very close (less than 400 m) to the nodes included in the line plan. This implies that all the passengers are near a bus stop in the obtained solution, which is an important condition in providing a public service.
• The directed demand (percentage of passengers that can execute their trips without transferring) reaches $42.6\%$. Note that this percentage is very high given the non-structured demand mobility pattern considered.
• The undirected demand (percentage of passengers that need to perform at least one transfer) is $56\%$ of the total demand. Although this percentage may seem relatively high, it is worth mentioning that, according to the 2021 statistical report of the public bus operator, the number of single-trip cards (10 trips without transfer) sold during the studied year was about $23.44\%$ and the number of single-trip tickets was $10.6\%$. Taking into account the rest of the ticket types, the total number of trips with transfers in the network reaches annually $50\%$, which is quite similar to the result obtained with our unfavorable demand matrix for a peak-hour scenario.
• The total covered demand (percentage of passengers that can take a bus to complete their trip from their origin to their destination) reaches $98.67\%$, i.e., the model produces a set of lines that traverse the nodes with higher demand, and the design of lines (including the frequencies assigned to each line) is sufficient to cover the trips in less time than using the pedestrian mode.
• The rest of the passengers’ trips corresponding to the uncovered nodes are $1.33\%$ but all of them have a bus stop within 400 m. In this sense, the transit network covers $100\%$ of the passenger demand.
• The average travel time (including waiting and transfer times) is 32.7 min, approximately 2 min less than the current average travel time in the city, which indicates the quality of the solution obtained after applying the proposed methodology.

Figure 6. The line plan for the Seville experiment.

Figure 6. The line plan for the Seville experiment.

Figure 7 and Figure 8 show the final assignment of passengers on the transit network. The width of the lines used to represent the flows is proportional to their values. Although the OD demand matrix does not reflect the real behavior of passengers (the current trends of movement in the city), as previously mentioned, the global results exceed the values of several transport indicators for the city.

8. Conclusions

In this paper, the use of a multilayer network structure to model and solve the line planning problem has been examined. The main objective consists of reformulating the line planning problem using a bimodal network composed of a pedestrian layer and a bus superlayer that contains as many layers as lines in the set of candidate lines. The inclusion of the pedestrian layer, which allows users to move from their origins to their destinations, ensures feasibility even when the capacity of the transport system is not sufficient.

An extensive review of previous works in the literature has been presented, classifying the contributions into different groups according to the means of approaching the problem.

A relatively simple mechanism is proposed to compute the pool of candidate lines. For each OD pair, the method first computes a predefined number of the best shortest paths that must fulfill certain length conditions. Then, the list of paths is ordered according to the direct demand that they can move and finally pruned to remove the worst paths.

The pedestrian network and the candidate lines are used to construct the multilayer network that later will be used to formulate the line planning problem. The use of the multilayer structure drastically simplifies the process of modeling waiting and transfer times, avoiding the need to consider the strong non-linearities that appear in classical formulations, thus reducing the complexity of the model.

Since, for each OD pair, the number of trip alternatives over the network is usually quite large, in order to obtain a sparse model and with the objective of working with medium-sized scenarios, only a subset of efficient paths, the strategy subgraph, is considered before performing the transit assignment.

Two mathematical formulations, arc-based and path-based, are presented. In both formulations, a first objective function that minimizes the total travel time is considered. In the second formulation, the use of path flow variables allows us to define alternative objective functions in terms of the average travel times, as well as to introduce explicit expressions to measure the number of passengers transferring one, two or more times.

For validation purposes, a first set of experiments was performed using the Mandl network as a benchmark scenario. Although our formulations do not impose a specific number of lines but an upper bound, in order to make a fair comparison, the problem was solved for 4, 6 and 8 lines, and the results were compared to those reported by 19 works extracted from the scientific literature, most of them using metaheuristics. For this comparison, several common performance indicators were considered: the percentage of demand that is served directly, the percentage of demand that is served with one transfer, the percentage of demand that is served with two transfers and the percentage of unattended demand. In the four-line experiment, our approach generated a solution with a value of the attended demand that was $3\%$ worse than the best reported solution and outperformed 15 of the 19 approaches analyzed. In the six-line experiment, our solution was $1\%$ worse than the best reported one, outperforming 17 of the previous approaches. For the 8-line experiment, the results were similar to the 4-line case.

To test the applicability of our formulations in larger scenarios, we performed a new experiment based on the Seville transit network, a network with 141 nodes, 454 links and 19,440 OD pairs. For this experiment, we considered an OD matrix with a total number of trips slightly larger than that corresponding to the rush hour period, but using a random distribution of trips, thus considering an unfavorable situation where no clear mobility patterns were represented. Despite the size of the problem, we obtained optimal solutions in a reasonable computation time, about 800 s. Our results exceed the values of several mobility indicators of the city, even considering an unfavorable demand scenario.

We believe that the use of a multilayer network representation can help to model and solve several problems in the field of transportation where different transport modes are involved. For instance, it is possible to analyze the expansion of a specific mode supposing that the rest of the transport modes remain invariable, but considering the interactions between the different modes. Other problems, such as the location of interchanges and the design of bicycle networks, seem to also be approachable using the methodology presented in this paper.