Models and a solution algorithm for planning transfer synchronization of bus timetables

Transportation Research Part E 131 (2019) 247–266

Contents lists available at ScienceDirect

Transportation Research Part E journal homepage: www.elsevier.com/locate/tre

Models and a solution algorithm for planning transfer synchronization of bus timetables

T

James C. Chu, Kanticha Korsesthakarn, Yu-Ting Hsu , Hua-Yen Wu ⁎

Department of Civil Engineering, National Taiwan University, Taiwan

ARTICLE INFO

ABSTRACT

Keywords: Bus timetabling Passenger transfer Transfer synchronization Mixed-integer programming Heuristic algorithm

Mixed-integer linear programming models and a heuristic algorithm are proposed for the transfer synchronization planning of bus timetables. The most important novelty of the proposed methodology is that the bus timetables and passenger choices of travel paths are simultaneously optimized in the model. The possibilities of transfer synchronization are fully explored in the timetable optimization because passengers select their paths to the destination in response to the timetables. Hence, bus timetables can be planned with high precision, and bus timetabling can be improved.

1. Introduction Timetabling is one of the major tasks in the design of bus systems. In any bus system of medium size or larger, the timetable should allow for the fact that many passengers will require two or more bus routes to reach their destination. In particular, a transfer is required when a direct service for an origin and destination (OD) pair is unavailable. In a dense network with carefully designed transfer stops and demand sufficient to allow bus service to be provided at high frequency, passengers may be able to transfer from one route to another with short waiting periods. However, in other situations, as considered here, with sparser routes and longer headways, transfers may involve longer waiting periods, and can make the service unattractive to passengers. Therefore, limiting the number and duration of transfers required to complete any given trip is desirable. Another important factor that must be considered in bus timetabling is the capacity of buses. Bus boarding will fail if the bus that the passengers are attempting to board at an origin or transfer stop is full. This situation is highly undesirable because it causes a delay in service for passengers, and hence late arrival at their respective destinations. The optimization of transit timetables has been studied extensively in literature. Comprehensive reviews of transit timetabling problems were conducted by Desaulniers and Hickman (2007), Guihaire and Hao (2008), and Ibarra-Rojas et al. (2015). Selected studies on transit timetabling problems are reviewed below. The tackled problems include different types of transit systems and are not limited to buses. Many of the studies on transit timetabling problems focused on the efficient operation of buses or trains and/or serving the travel demands of the entire system; transfers were not explicitly considered (Ceder, 2007; Ceder et al., 2013; Wang et al., 2017, 2018; Zhang and Xu, 2017; Yang et al., 2017; Chen et al., 2015). Other studies aimed to improve the transfer timing in transit systems. In this respect, timetabling problems can be grouped into two major categories. One category focuses on transit operations, particularly the synchronization of bus arrivals at transfer stations. Transfer synchronization refers to the timing of services to provide convenient interchange for passengers transferring from one route to another. Examples of studies under this category include those of Ceder, Golany, and Tal (2001), Ibarra-Rojas and Rios-Solis

⁎

Corresponding author at: No. 1, Sec. 4, Roosevelt Road, Taipei 10617, Taiwan. E-mail address: [email protected] (Y.-T. Hsu).

https://doi.org/10.1016/j.tre.2019.10.001 Received 27 June 2018; Received in revised form 7 October 2019; Accepted 9 October 2019 1366-5545/ © 2019 Published by Elsevier Ltd.

Transportation Research Part E 131 (2019) 247–266

J.C. Chu, et al.

(2012), Guo et al. (2017, 2018), and Kang et al. (2015). However, studies in this category have not explicitly considered passenger demand or vehicle capacity, leading potentially to unrealistic outcomes. The other category of problems focuses on reducing the inconvenience of transfers for passengers. In these problems, the number of successful transfers is maximized, and the inconvenience of transfers, usually the passenger waiting times, is minimized. This category includes the studies of Schroder and Solchenbach (2006), Wu et al. (2015), Kwan and Chang (2008), Tilahun and Ong (2012), Kang et al. (2019), Yin et al. (2018), Wong et al. (2008), Sun et al. (2019), and Xiong et al. (2015). Most of the models in these studies formulated the minimum walking time from one bus or train stop to another in the event of a transfer, indicating that this factor is critical in timetabling. Except for a few studies (Kwan and Chang, 2008; Sun et al., 2019; Xiong et al., 2015), capacity was assumed to be sufficient for all buses/trains in other studies. Notably, passengers’ choices of bus routes, including transfer, are usually assumed and regarded as constant regardless of the timetable design. In reality, when multiple choices of bus routes and transfer stops are available, passengers decide on travel paths based on timetables. The traveling decisions of passengers have long been recognized in the planning and management of transportation systems in classical problems, such as traffic assignment, network design, and signal setting (Migdalas, 1995; Wardrop, 1952; LeBlanc, 1975; Marcotte, 1983). However, only a few studies in bus timetabling have incorporated a transit assignment model into the transfer synchronization problem to consider changes in passenger demand patterns in response to timetable design. An example is the study of Xiong et al. (2015), in which two notable assumptions were made to simplify the synchronization of transfers. First, all passengers are traveling from shuttle bus stops to a small number of metro stations, implying that passengers are traveling in the same general direction. Second, the time required to make a transfer (i.e. to alight from the first bus, and walk to and board the second) is assumed zero. To the best of the authors’ knowledge, a transit timetabling problem that allows passengers to respond to a timetable design without assuming passengers traveling in the same direction and zero transfer time has not been studied yet. This paper presents a bus timetabling model that incorporates details of passenger behavior, with a particular focus on the synchronization of services at transfer points. The model is applicable to networks with relatively low demand and high headways because synchronization to minimize waiting times is clearly of minimal consequence when headways are less than the maximum acceptable waiting times. In this study, a bus transfer is considered synchronized when the departure time of a bus is later than the arrival time of another bus, with allowance for the time required to make the transfer. The most important novelty of the proposed methodology is the simultaneous optimization of timetables and passenger paths by the model. Passengers select their paths to the destination in response to the timetables. The feasible paths that can be selected by passengers are restricted by criteria representing travel preferences, such as total number of transfers, maximum waiting time, and maximum total travel time, to produce reasonable travel behavior. The passenger paths generated by the proposed methodology are also useful for accurately estimating the transfer demand at transfer stops. The information can be used for guidance in all stages of bus system planning and operations, such as adjustment of bus routes and frequencies and determination of the capacity of passenger facilities at transfer stops. The present model is formulated as a mixed-integer linear programming (MILP) problem. A notable assumption in the model formulation may affect the practicality of planning results. The path choice behavior of passengers is determined by a transit assignment model proposed by Carraresi et al. (1996), which is selected because of its simple assumption on passenger behavior; it can also be represented by a network flow model and formulated with linear constraints. The transit assignment model is coupled with the minimization of the total passenger travel time. The implied assumption of this objective is that passengers can find the shortest travel time path among feasible ones during traveling. This assumption is valid when the complete information of the bus system is available to the passengers. Another factor that is disregarded in the transit assignment model is the crowding condition of buses, that is, the discomfort experienced by passengers when passenger loads are near full capacity. Note that crowding in this sense is distinct from failure of service in a situation where a driver cannot pick up passengers because the bus is already full (this is dealt with explicitly in the proposed model). A solution approach, which includes a set-partitioning (SP) reformulation and a heuristic algorithm based on the reformulated models, is developed to solve the problems. In addition, the technique of big-M reduction is adopted to improve the problem-solving efficiency. This paper is organized as follows. In the next section, the assumptions and the formulations of the optimization model for transfer synchronization are explained. The solution approach is proposed in Section 3. In Section 4, numerical examples used to validate the proposed model are discussed. Section 5 provides the summary of findings and the suggestions for future research. Appendix A describes the technique of big-M reduction, which improves the efficiency of problem-solving. Appendix B specifies and gives sources for the parameters used in Section 4. 2. Mathematical model This section describes the problem statement, defines the notations, and proposes the mixed-integer linear programming model. 2.1. Problem statement A graph is used for representing the physical road network, wherein the links are the roads on which the buses travel, and the nodes are the bus stops. A bus route is defined as a sequence of links and bus stops connected by those links. The bus routes and their frequencies are known and constant, whereas the headways and running times between stops can be adjusted within predetermined bounds. The dwelling time of a bus at a stop is constant. A fleet of buses with equal passenger capacity is considered. The timetable of a bus route describes the arrival and departure times at all stops for all bus dispatches of the route. An origin and destination (OD) group 248

Transportation Research Part E 131 (2019) 247–266

J.C. Chu, et al.

represents passengers with the same OD and desired departure time. The size of an OD group is limited to the capacity of a bus: if the number of passengers of an OD group initially exceeds the bus capacity, then the passengers are divided into multiple groups according to the bus capacity. The passengers can transfer between buses provided sufficient time for transfer is available. The sequence of buses that an OD group takes to complete the travel is called a passenger path or path. All the passengers of an OD group are assumed to take the same path. Passengers choose their paths to the destination in response to the timetables. Following the transit assignment model by Carraresi et al. (1996), each OD pair has a set of feasible paths. A path is feasible to a passenger if the following conditions are satisfied.

• The total travel time does not exceed the maximum acceptable total travel time of the OD group. • The number of transfers does not exceed the maximum number of total transfers. • The initial waiting time at the origin stop does not exceed the maximum acceptable initial waiting time. • The waiting time at any transfer stop does not exceed the maximum acceptable transfer waiting time. Any solution in which all OD groups complete their trips with feasible paths while satisfying the bus capacity constraint is a feasible instantiation of the passengers’ path choices (Carraresi et al., 1996). The objective in the timetable optimization model proposed here is to minimize total travel times of passengers. The model ensures that all passengers accepted for travel can reach their destinations as required by time-window and other passenger-travel requirements, with penalties applied where passengers are not accepted. The model also ensures that a passenger waiting at a bus stop can always take the first available bus: this, with the bus capacity constraint, ensures that service is never denied due to buses running at, or close to, full capacity. The outputs of the model are the bus timetables and the chosen passenger paths. 2.2. Mixed-integer programming model formulation Here, the MILP formulation of the timetable optimization model is presented. The following notations of indices, sets, parameters, and variables are adopted in the MILP model. Note that several variables (Sg0 , Eg0 , dg0 , lg0 ) are involved in the construction of dummy routes, which are assigned notionally to OD groups for which travel is infeasible. Sets S: Sl :

the set of physical bus stops, indexed by i, j, s. the set of physical bus stops on route l, indexed by i, j, s. the set of dummy stops of OD group g, indexed by i, j, s.

E: El :

the set of physical network links, indexed by a pair of stops. the set of physical links on which bus route l travels, indexed by a pair of stops. the set of dummy links for OD group g, indexed by a pair of stops.

L: Ps : G:

the set of bus routes, indexed by l, r. the set of route pairs that intersect at stop s, indexed by (l, r ) the set of OD groups, indexed by g.

Sg0 :

Eg0 :

Ps , l

r, s

Parameters c: d: dg :

the bus capacity. the dwelling time at a stop. the destination stop of OD group g.

dg0 :

the dummy destination stop of OD group g.

hllb :

the lower bound of the bus headway on route l.

fl :

hlub :

the total number of bus dispatches of route l (i.e., dispatch frequency). the upper bound of the bus headway on route l.

k:

the maximum number of total transfers. the dummy route of OD group g.

M: ng : og : pg :

a large positive number. the number of passengers in OD group g.

rijub :

the upper bound of the bus running time from stop i to stop j.

lg0 :

q: rij :

rijlb :

tg : ug : w: zg :

gs :

the origin stop of OD group g.

the penalty for an unserved passenger in OD group g.

the maximum acceptable initial waiting time. the average running time of a bus from stop i to stop j.

the lower bound of the bus running time from stop i to stop j. the desired departure time of OD group g.

the maximum acceptable total travel time of OD group g.

the maximum acceptable transfer waiting time. the sink stop of OD group g.

the minimum transfer walking time of OD group g at stop s.

249

S.

Transportation Research Part E 131 (2019) 247–266

J.C. Chu, et al.

Decision variables Bus-related B : Dsln Hsln : Rijln :

Passenger-related P Ags :

P Dgs :

the departure time of the nth bus dispatch on route l at stop s. the headway of the nth bus dispatch on route l at stop s. the running time of the nth bus dispatch on route l from stop i to stop j. the arrival time of OD group g at stop s. the departure time of OD group g at stop s.

Sgs : Wgs : Xgij :

the total time that OD group g spends at stop s (transfer plus waiting).

Zgslr :

1 if OD group g transfers from route l to route r at stop s; 0, otherwise.

Ygijln :

the transfer waiting time of OD group g at stop s.

1 if OD group g travels from stop i to stop j; 0, otherwise.

1 if OD group g takes the nth bus dispatch of route l from stop i to stop j; 0, otherwise.

Before introducing the model formulation, the network used to define the bus routes and feasible paths of passengers in the model formulation is described. This network is designed to allow formulation as a network flow model. A network contains physical bus route links and stops shared by all OD groups and OD group-specific dummy links and stops. Fig. 1 shows an example of such a network. The solid lines and circles respectively represent the physical bus route links and stops shared by all OD groups. The arrows of the lines indicate the route directions. For OD group g with origin stop og and destination stop dg , the feasible paths of OD group g have two types. First, for the paths dg z g and the where OD group g takes physical buses to complete the trip, the feasible paths are defined by og 1 2 z g is physical buses that are taken, where 1, 2 , … are the first and second intermediate stops and so on, and z g is the sink node. dg the dummy collection link for OD group g, which is defined for all g G . Second, for the paths where OD group g cannot complete the dg0 z g , with dummy destination stop dg0 . og dg0 is the dummy bus trip with the physical buses, the feasible path is defined by og 0 z g is the dummy collection link for OD group g, which is defined for all g G . In Fig. 1, the dummy links and stops for link, and dg OD 1 → 7 are represented by dashed lines and circles, respectively. Those of other ODs are not displayed for simplicity. A single dummy bus is running for each of the dummy links. The dummy bus is not restricted by any timetables; therefore, traveling via the dummy buses easily satisfies all the requirements of a feasible path defined in Section 2.1. However, dummy buses are associated with a high penalty. Thus, an OD group will use the dummy bus to complete the trip only if the OD group cannot be served by the physical buses. Moreover, dummy destination stop dg0 is introduced to distinguish a trip via a dummy bus from a feasible single-link

Fig. 1. Bus route network (small). 250

Transportation Research Part E 131 (2019) 247–266

J.C. Chu, et al.

dg z g via the physical bus. For example, Stops 1 and 4 in Fig. 1 are adjacent. For OD group g with Stop 1 as the origin and path og Stop 4 as the destination, the path of Stop 1 (og ) → dummy destination stop (dg0 ) → sink node (z g ) takes the dummy bus while that of Stop 1 (og ) → Stop 4 (dg ) → sink node (z g ) utilizes a physical bus. Additional constraints are introduced to ensure that the bus timetables meet the required headways and running times (Section 2.2.1), the movements of the OD group constitute a feasible path (Section 2.2.2), and the movements of the OD group are consistent with the movements of buses (Section 2.2.3). The complete details of the objective function and the constraints follow. The objective function of the model (Eq. (1)) is to minimize the total travel time of the passengers. Unserved passengers (i.e. OD groups using dummy bus links) experience a very long travel time.

ng × (Agzp g

min

tg )

(1)

g G

2.2.1. Bus timetables Constraints (2)–(5) govern the timetables for all bus physical routes. The constraints are modified from those used by Wong et al. (2008). Constraint (2) defines the headway between departures. Constraint (3) requires that the headway must be between the minimum and maximum headways. Constraint (4) states that the departure time of a bus at stop j is the sum of the departure time at the previous stop i, a bus running time from stop i to stop j, and the dwelling time at stop j. Constraint (5) states that the bus running time on a link must be between its lower and upper bounds. The dummy bus routes do not have timetables. Each of the dummy routes has only one bus dispatch, i.e., flg0 = 1, g G . Constraint (6) states that the dummy bus link has a long travel time ( pg ) as a penalty for the unserved passengers, forcing OD groups to take a physical bus if possible, rather than using the dummy bus. Constraints (7) and (8) state that both dummy collection links have a travel time of 0. B Dsln = DslB(n

1)

hllb

hlub,

Hsln

+ Hsln,

l

l

rijub,

Rij ln

Rog dg0 lg0 1 = pg ,

Rdg z g ln = 0,

l g

g

Rdg0 z g lg0 1 = 0,

Sl, n

L, s

Sl, n

l

L , (i , j )

B B Djln = Diln + Rij ln + d ,

rijlb

L, s

L , (i , j )

(2)

{2,..,fl }

(3)

{1. .fl } El, n

(4)

{1, ..,fl }

(5)

El

G

(6)

G, l

g

L, n

{1, ...,fl }

(7)

G

(8)

2.2.2. Feasible passenger paths Passengers choose their paths from the feasible path set in response to the timetables through Constraints (9)–(18). Constraints (9) and (10) ensure the flow conservation on the bus route network for each passenger OD group. Constraint (11) states that the departure time of each OD group at a stop is the arrival time at the stop plus the transfer and waiting times at the stop. Constraint (12) states that the arrival time of an OD group at the origin stop is the desired departure time of the group, i.e., tg . According to Constraint (11), the initial waiting time (i.e., the waiting time at the origin stop, Sgog ) is the time difference between the desired and the actual Zgslr = 1, meaning that OD group g transfers from route l to route departure times at a stop. Constraints (13) and (14) state that if (l, r ) Ps

r at stop s, then the time the passengers spend at the stop (Sgs ) is the sum of the transfer walking time gs and the transfer waiting time Wgs . These constraints imply that when passengers transfer at a stop, the time spent at the stop must be greater than or equal to the minimum transfer walking time ( gs ). The determination of Zgslr is explained in Section 2.2.3. The remaining constraints guarantee that the paths taken by the passengers are feasible in terms of service quality. Constraint (15) states that the total time of the trip cannot exceed the maximum acceptable travel time that the OD group can tolerate. Constraint (16) ensures that the number of transfers does not exceed the maximum allowable number of transfers. Constraints (17) and (18) restrict the initial waiting time at the origin stop and the transfer waiting times at the remaining stops so as not to exceed their maximum allowable values.

Xgis = (i, s ) E Eg0

(og, j) E Eg0

Xgsj ,

Xgog j = 1,

AgsP + Sgs = DgsP ,

AoPg = tg ,

s

S

(s, j) E Eg0

g

s

g

S

Sg0

{og , z g }, g

G

(9)

G Sg0, g

(10) (11)

G

G

(12) 251

Transportation Research Part E 131 (2019) 247–266

J.C. Chu, et al.

M× 1

Zgslr + Sgs

gs

+ Wgs,

g

G, s

S

(13)

(l, r ) Ps

M× 1

Zgslr + Sgs

+ Wgs,

gs

g

G, s

S

(14)

(l, r ) Ps

Agzpg

tg

ug ,

g

Zgslr

G

k,

g

(15)

G

(16)

s S,(l, r ) Ps

Sgog

q,

g

G

Wgs

w,

g

G, s

(17)

S

(18)

{og , d g }

2.2.3. Connection between passengers and bus timetables The following constraints connect the passenger movements and the bus timetables. Constraint (19) states that if an OD group moves from stop i to stop j, it must choose one bus from the ones that also travel from stop i to stop j. Constraints (20) and (21) ensure that the departure time of the passengers at a stop is identical to the departure time of the chosen bus. Constraints (22) and (23) assert that if an OD group chooses to take a bus to move from one stop to another, then its arrival time at the next stop is its departure time Ygisln = 1 and at this stop plus the running time of the chosen bus. Constraints (24)–(26) state that l L,(i, s ) El, n {1,.., fl }

r L,(s, j) Er , n {1,.., fr }

Ygsjrn = 1 if and only if Zgilr = 1. That is, a transfer is identified at a stop for an OD group if the OD group arrives at the

stop by a route and leaves the stop by a different route. Constraint (27) imposes that the bus capacity is not violated. The remaining constraints of the model require that an OD group boards the first bus after the group arrives at the stop. Note that the effect of constraints (28)–(32) is to ensure that the capacity constraint does not delay service to passengers. In particular, when a bus traversing route r arrives at a stop where an OD group g is waiting for a bus on r, if the bus has insufficient capacity to take in all members of g, any further progress by g via r (e.g. by waiting for the next bus traversing r) is ruled out. P DoBg lm + 1 for all m {1, ..,n 1} . That is, if an OD group boards a bus at the Constraint (28) states that if Ygog j ln = 1, then Ago g origin stop, then the arrival time of the group must be later than the departure time of any bus preceding the boarded bus (of the same B Zgslr = 1 and Ygsjrn = 1, then AgsP + gs Dsrm + 1 for all m {1, ..,n 1} . gsrn is an route). Constraints (29)–(32) state that if (l, r ) Ps

auxiliary variable that is only used in constraints (29)–(32), which is 1 if passengers of OD group g take nth bus dispatch of route r at stop s via transfer and 0 otherwise. That is, if an OD group boards a bus via transfer, then the arrival time plus the walking transfer time of the group must be later than the departure time of any bus preceding the boarded bus (of the same route).

Xgij =

Ygij ln,

g

G , (i , j )

Eg0

El

l L {lg0}, n {1,.., fl }

Ygij ln ) + DgiP

M × (1

Ygij ln ) + DgiP

M × (1

M × (1

B Diln ,

Ygij ln ) + AgjP

M × (1

B Diln ,

Ygij ln ) + AgjP

g

g

G, l

Rij ln + DgiP ,

Rij ln + DgiP ,

Ygisln

G, l

Zgslr ,

g

L, n

L, n g

g G, s

{1, ..,fl }, (i, j)

{1, ..,fl }, (i , j )

G, l

G, l

(19)

L

S, (l, r )

{lg0}, (i , j )

(21)

El

{lg0}, (i, j)

L

(20)

El

El

El

Eg0, n

Eg0, n

{1, ..,fl }

{1, ..,fl }

Ps

Zgslr ,

g

G, s

S , (l , r )

Ps

(25)

r L,(s, j) Er , n {1,.., fr }

Ygis ln + l L,(i, s ) El, n {1,.., fl }

(Ygij ln × ng )

Ygsjrn

Zgslr + 1,

g

G, s

r L,(s, j ) Er , n {1,.., fr }

c,

l

L , (i , j )

El, n

{1. .fl }

g Ygsjrn

G, r gsrn,

P Ygog jrn ) + Ago g

L , (og , j) g

G, s

Er , n S

S, (l, r )

Ps

(26) (27)

g G

M × (1

(23) (24)

l L,(i, s ) El, n {1,.., fl }

Ygsjrn

(22)

DoBg rm + 1,

{2, ..,fr }, m {og , dg }, r

{1, ..,n

L , (s , j )

(28)

1} Er , n

{2, ..,fr } 252

(29)

Transportation Research Part E 131 (2019) 247–266

J.C. Chu, et al.

Zgslr

gsrn,

g

G, s

S

{og , dg }, r

L, n

{2, ..,fr }

G, s

{og , dg }, r

(30)

(l, r ) Ps

Ygsjrn +

Zgslr

gsrn

+ 1,

g

S

L , (s , j )

Er , n

(l, r ) Ps

M × (1 g

G, s

S

gsrn )

+ AgsP +

{og , d g }, r

(31)

B Dsrm + 1,

gs

L, n

{2, ..,fr }

{2, ..,fr }, m

{1, ..,n

(32)

1}

3. Solution approach This section introduces the solution approach to the problem described in Section 2. The solution approach includes a setpartitioning (SP) reformulation and a heuristic algorithm based on the reformulated models. In addition, the technique of big-M reduction is adopted to improve the problem-solving efficiency. The application details of the technique on the original formulation in Section 2.2 and the reformulated one in Section 3.1 are explained in Appendix A. 3.1. Problem reformulation Inspired by the set-partitioning reformulation of vehicle routing problems (Balinski and Quandt, 1964), the solution approach first enumerates the feasible paths (sequences of buses) for all OD groups and then reformulates the model given these paths. The two major advantages of the set-partitioning reformulation are that it generally provides tight lower bounds and avoids the modeling of route constraints when they are complex (Toth and Vigo, 2015). The procedure of reformulating the model into a set-partitioning one is described in Section 3.1.1. The algorithm for enumerating the feasible paths is presented in Section 3.1.2. 3.1.1. Set-partitioning model The set-partitioning formulation is logically equivalent to the original formulation given in Eqs. (1)–(32). It involves the dropping of some constraints from the original formulation, and the addition of some new ones. The additional constraints are specified in Eqs. (33)–(43). The reformulated model then comprises Eqs. (1)–(8), (11), (12), (15), (17), (18) and (33)–(43). Y Given the passenger paths, the buses on the kth path of OD group g are collected in set Cgk and indexed by (i,j,l,n). The stops where Z the transfers are made and the connecting route on the kth path of OD group g are stored in set Cgk and indexed by (s,r). A newly introduced binary variable, gk indicates whether the kth path is chosen by OD group g or not; 1 if chosen and 0 if not. In the reformulated model, the objective function (Eq. (1)) and the constraints that do not involve Xgij , Ygijln , and Zgslr remain the same Y (Constraints (2)–(8), (11), (12), (15), and (17), (18)). The constraints that involve Ygijln and Zgslr are reformulated using gk , Cgk , and Z Cgk . More specifically, Constraints (13), (14), (20)–(23), and (27)–(32) are replaced by Constraints (33)–(42), where Ng is the number of paths of OD group g. Constraint (43) is a newly added constraint to ensure that only one path is chosen for each OD group. The rest Y Z of the constraints in the original model (Constraints (9), (10), (16), (19), (24)–(26)) are implicitly considered in Cgk and Cgk , so they are not included in the reformulated model. The reformulated set-partitioning problem can be directly solved as a MILP problem and can also be solved with the heuristic algorithm proposed in Section 3.2.

M × (1

M × (1

gk )

gk )

M × (1 M × (1

M × (1 M × (1

+

gk ) gk )

gs

+ Sgs

gk ) gk )

+ Sgs gs

+ DgiP DgiP

+

+ Wgs,

+ Wgs, B Diln ,

B Diln ,

AgjP

+ AgjP

g

g

G, k

g g

G, k G, k

Rij ln +

DgiP ,

Rij ln + DgiP , (n g ×

G, k

gk )

Z Cgk

{1, ..Ng }, s

{1, ..Ng}, s

Z Cgk

{1, ..Ng }, (i , j, l, n) {1, ..Ng }, (i, j, l, n)

g

G, k

g

G, k

c,

l

(34) Y Cgk

{1, ..Ng }, (i, j, l, n) L , (i , j )

El, n

(35)

Y Cgk

{1, ..Ng }, (i , j, l, n)

Y g G, k {1,.. Ng},(i , j , l , n ) Cgk

(36) Y Cgk Y Cgk

(37) (38)

{1. .fl } (39)

s . t . (i , j , l , n ) = (i, j,l, n )

M× 1

(33)

Y k {1,.. Ng },(s , j , r , n ) Cgk

gk

P + Ago g

DoBg rm + 1,

s . t . s = og , r = r , n = n

g

G, r

L, n

{2, ..,fr }, m

{1, ..,n

(40)

1} 253

Transportation Research Part E 131 (2019) 247–266

J.C. Chu, et al.

gk

=

gsrn,

g

G, s

S

{og , d g }, r

L, n

Y k {1,.. Ng },(s , j , r , n ) Cgk

{2, ..,fr } (41)

Z, s . t . s = s, r = r , n = n,(s , r ) Cgk

M × (1 g

G, s gk

S = 1,

gsrn )

+ AgsP +

{og , d g }, r g

L, n

B Dsrm + 1,

gs

{2, ..,fr }, m

{1, ..,n

(42)

1}

G

(43)

k {1.. Ng }

3.1.2. Generation of passenger paths Algorithm 1 is proposed to enumerate the feasible paths efficiently. An illustrative example of the algorithm can be found in Section 4.1.1. The algorithm utilizes a tree structure and is developed from the breadth-first search and path enumeration algorithms in graph theory. A node of the tree stores a tuple representing the previous stop, the current stop, the arriving route, the arriving bus, the lower bound of passenger arrival time, the upper bound of passenger arrival time, and number of transfers so far. Line 1 of the algorithm prepares the two data inputs for enumeration. First, the network links (physical and dummy) are prepared for all OD groups. Second, given the upper and lower bounds of the headway and bus running time, the time windows of the departure of every B B bus at every stop are calculated. Dsln and Dsln , the lower and upper bounds of the departure time window of the nth bus dispatch on route l at stop s, respectively, are calculated in Appendix A. After the data are prepared, the following steps are conducted for each of the OD groups. Line 3 states that T0 stores the tree nodes that have not been considered yet and is initialized with the tree node that represents the origin stop. Line 4 states that Tg stores the tree nodes that have been added to the path enumeration tree of OD group g and is initialized as an empty set. In the loop in Lines 5–20, tree nodes are generated to form passenger paths. Before a tree is added to Tg, the tree is checked by the three following conditions to ensure that the generated paths are feasible for the OD group. (1) The node represents the origin stop of the OD group, and the time window of the departing bus overlaps the time window of the OD group (Lines 11–13). (2) The node does not represent the origin stop, and no transfer occurs (arriving and departing on the same bus) (Lines 15–17). (3) The node does not represent the origin stop of the OD group, a transfer occurs (arriving and departing on different buses), the number of transfers is smaller than the limit, and the time window of the departing bus overlaps the time window of the OD group (Line 18–20). At the end of the algorithm, the path enumeration trees are constructed for all OD groups. The trees are used in the SP model introduced in Section 3.1.2 and the heuristic solution procedure in Section 3.2. Algorithm 1 (Passenger path enumeration for OD pairs). 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21

INPUT: Eg = E + Eg0,

g

B B , Diln ], time windows of departure times at all nodes for all buses and routes G . [Diln

for g G do T0:={(-,o,-,-,−∞,+∞,0)}. Tg:={} while T0 not empty do:

From T0, select tree node α=(previous stop, i, larr, narr, [AiP , AiP ], xfer). Remove it from T0 and add to Tg. for (i, j) in Eg: for l ∈ L:: if network link (i, j) is passed by route l: for n {1, ...,Fl} : B B , Diln ] if i = o and [Diln

then

[tg, tg + q]

θ:= (i , j, l, n, AiP + d + rijlb, AiP + d + rijub, xfer )

T0:= T0 { θ } elif i ≠ o then if (l,n) = (larr, narr) then

θ:=(i , j, l, n, AiP + d + rijlb, AiP + d + rijub, xfer )

T0:= T0 { θ }

B B , Diln ] elif xfer < k and [Diln

θ:=(i , j, l, n,

AiP

+

gi

+

rijlb,

, narr

AiP )

d+

+w+

rijub,

xfer + 1)

[max( D Barr il

AiP

+

gi

T0:= T0 { θ} OUTPUT: Path enumeration tree Tg for all OD groups

gi,

min (D Barr

254

il

, narr

AiP )

d+

gi

+ w]

then

Transportation Research Part E 131 (2019) 247–266

J.C. Chu, et al.

3.2. Heuristic procedure to solve the reformulated models Although the reformulated problem can be directly solved as a MILP one, this problem can still be computationally difficult due to a large number of alternative paths of each OD group. A heuristic solution procedure that combines incremental assignment and local search (Algorithm 2) was developed as an alternative to solve the reformulated model via MILP. The procedure first solves the problem incrementally; that is, the problem in each iteration is only solved for a portion of OD groups, and the paths of those OD groups are fixed. Assigning OD groups incrementally is particularly effective for large problems, in which the number of OD groups is large and solving the problems is time-consuming. After the paths of all OD groups are determined, a local search procedure is conducted to improve the solution by repeatedly resolving the path choice variables of OD groups with common destinations or origins. The motivation of the procedure is that passengers with common origins or destinations are likely to compete for the same bus service. Resolving the timetable and paths of these conflicting OD groups has the potential for improving the solution. The case of resolving OD groups of two destination stops in Fig. 1 is considered to illustrate the local search procedure. Let Nodes 1, 2, 6, and 7 be the origin and destination stops of passengers. The combinations of two destination stops include (1,2), (1,6), (1,7), (2,6), (2,7), and (6,7). In each iteration, the paths of OD groups that correspond to one of the combinations (as well as the timetables) are resolved. For example, given the combination of (1,2), all the paths of OD groups that are destined for Stop 1 or 2, including those of ODs 2 → 1, 6 → 1, 7 → 1, 1 → 2, 6 → 2, and 7 → 2, are resolved. The greater improvement can be expected when more paths are relaxed; however, relaxing more paths leads to a longer solution time. The maximum number of destination or origin stops of which OD groups are resolved should be determined depending on the problem scale. The details of the heuristic method are described in Algorithm 2. In Step 1, the paths are enumerated. In Step 2, the OD groups are sorted in decreasing order according to the number of passengers per group. Steps 3–5 state that if any OD group is unassigned, then some OD groups are selected in decreasing order, and the problem is solved for these OD groups. The corresponding path choice variables ( gk ) are fixed after solving the problem. These steps are repeated until all the OD groups are assigned to paths and all the path choice variables are fixed. In Step 6, a feasible solution, that is, the current best solution, is identified. In Step 7, the number of destination or origin stops (b) of which OD groups are resolved is set, increasing from 1 to B. In Step 8, all the combinations of origin and destination stops are considered unselected. In Steps 9–11, the OD groups are selected based on the current combination of destination stops. The timetable and path choice variables of the selected OD are relaxed, and the problem is solved again. If the objective value is unimproved, then the ODs of the next combination of destination stops are resolved. If the objective value is improved, then the new solution is adopted as the best one in Step 10 and the algorithm goes back to Step 8. Steps 12–14 are similar to Steps 9–11 but based on origin stops. In Step 15, if no improvement is found after all combinations of destination and origin stops are selected, then the algorithm terminates. Finally, considering that only an improved solution will be adopted, the solution after the improvement must be at least as good as the original one. Algorithm 2 (Heuristic algorithm). Input:

Bus routes and OD groups

Step 1 Step 2 While

Enumerate paths for all OD groups (Algorithm 1). Sort OD groups in decreasing order according to their size At least one of the OD groups is not added to the problem. Step 3 OD groups are selected in decreasing order. If the number of remaining OD groups is less than , select all of them. Step 4 Add the selected OD groups to the problem and solve the problem. Step 5 Fix the values of the path choice variables ( gk ).

Step 6 Adopt current solution as the best solution. Step 7 For b = 1 to B (the maximum value of stops of which OD groups are resolved in each iteration): Step 8 All the combinations of b destination stops and those of b origin stops are marked unselected. While At least one of the combinations of b destination stops is unselected. Step 9 Select an unselected combination of b destination stops. Select OD groups that are destined for the selected destination stops. Step 10 Relax the timetable and path choice variables of the selected OD groups. Solve the problem. Step 11 If the objective value is improved, then adopt the current solution as the best solution and proceed to Step 8. While At least one of the combinations of b origin stops is unselected. Step 12 Select an unselected combination of b origin stops. Select OD groups that depart from the selected origin stops. Step 13 Relax the timetable and path choice variables of the selected OD groups. Solve the problem. Step 14 If the objective value is improved, then adopt the current solution as the best solution and proceed to Step 8. Step 15 Terminate. Output: Bus timetables and path of each OD group

4. Computational tests The proposed methodology is validated using two examples. A small network is studied in Section 4.1, and a large example is studied in Section 4.2. The MILP models and the algorithms are respectively developed in Python Programming Language version 2.7 and Gurobi Optimizer version 8.1. The hardware environment for the numerical example is Intel Core i7 3.4 GHz CPU and 32 GB RAM on Microsoft Windows 8. 255

Transportation Research Part E 131 (2019) 247–266

J.C. Chu, et al.

Table 1 Number of bus dispatches (small network).

No. of dispatches

Route

S

S*

G

G*

K

K*

Low frequency High frequency Mixed frequency

5 10 5

5 10 5

5 10 5

5 10 5

5 10 10

5 10 10

Three solution methods are tested in this section: the original model introduced in Section 2 solved via MILP solver, denoted by Original (MILP), the reformulated SP model introduced in Section 3.1 solved via MILP solver, denoted by SP (MILP), and the reformulated SP model solved by the heuristic algorithm, denoted by Heuristic. 4.1. Small network The proposed methodology is validated using the bus route network shown in Fig. 1. The network includes six routes (S, S*, K, K*, G, and G*) and seven stops, two of them are transfer stops (4 and 5). For the investigation of the influence of the important parameters on the model results and computational times, three levels of bus frequency, two types of passenger flow, and two values of transfer walking time, which total eight scenarios, are considered. First, three levels of bus frequency, namely, low frequency, high frequency, and mixed frequency, which are defined in Table 1, are considered. The main purpose of exploring these three levels of frequency is to understand the importance of the frequency values of bus routes. Second, two types of passenger flow are considered, namely, unidirectional and bidirectional. In the unidirectional flow, only one direction between two stops has demand. In the bidirectional passenger flow, both directions have demand between two stops. The purpose of this analysis is to understand the effect of direction on the ease of timetabling. Third, the necessity of modeling the transfer walking time of passengers is considered by comparing models without explicit transfer walking time constraints and the ones with constraints. The passenger demand is assumed as follows. In the bidirectional demand scenarios, the following 12 OD pairs exist: 1 → 6, 6 → 1, 1 → 7, 7 → 1, 2 → 6, 6 → 2, 2 → 7, 7 → 2, 3 → 1, 1 → 3, 3 → 2, and 2 → 3. In the unidirectional demand scenarios, the following six OD pairs exist: 1 → 6, 1 → 7, 2 → 6, 2 → 7, 3 → 1, and 3 → 2. Notably, all the above OD pairs require a transfer to complete the trips. The OD pairs without transfer requirement to complete the trips are hardly affected by the optimization of transfer synchronization; for example, OD 1 → 4 can be served by Route S without transfer. These pairs are excluded from the example for the ease of analysis and presentation. All the passengers of an OD pair are grouped based on the largest number of bus dispatches of all the bus routes in all the scenarios, that is, 10 groups. The departure times of OD groups are evenly distributed within the day. The number of passengers of each OD group is uniformly distributed from 8 to 15, and the total number of OD groups and passengers is listed in Table 2. Other parameters of the example are defined in Appendix B 4.1.1. Computational study (small network) The parameters adopted in the solution approaches are as follows. The solution time is limited to 1 h for all the solution approaches. Based on the network in Fig. 1, Fig. 2 shows the path enumeration tree for the OD group from Stops 1 to 7 departing at 12:00 and the tree containing four paths as an example. In the heuristic algorithm, the number of OD groups to be added in each iteration (α) is 12, and the maximum value of stops of which OD groups are resolved in each iteration (B) is 2. The computational results for the eight scenarios solved by the three solution approaches are summarized in Table 2. The big-M reduction introduced in Appendix A is adopted in the aforementioned approaches. Appendix A also shows the computational results without big-M reduction. The results demonstrate that the optimality of the solution can be proven for only some of the scenarios. The difficulty of proving the optimality of the solutions is due to the weak lower bounds caused by several logical constraints, such as Eqs. (13), (14), (20)–(23), (28), and (32), and their equivalents in the reformulated model. This issue has been observed in similar classes of problems, such as the capacitated network design problems considered in Magnanti and Wong (1984), Yang and Bell (1998), and Hewitt et al. (2010, 2013). Moreover, optimality is not always proven in these problems. The three approaches are compared in this study. According to the comparison result of Original (MILP) and SP (MILP), given 3600 s of solution time, the Original (MILP) has an average of 27.79% optimality gap while that of SP (MILP) has 28.91%. However, SP (MILP) finds the better feasible solutions in most of the scenarios, thereby indicating the benefit of the reformulation. The heuristic algorithm has better performance in efficiently finding good feasible solutions than the MILP models. The CPU times are between 47.7 s and 309.0 s. The differences between the heuristic and the best feasible solutions identified by the two MILP models are minimal. The largest difference is 6.51%, and nine of the differences are within 5.00%. The computational times of the heuristic algorithm in finding the final feasible solution are much shorter than those of the MILP models. The results are analyzed from the following three aspects: frequency, direction of demand, and transfer time. First, the frequency of bus service has a remarkable effect on the difficulty of solving problems in terms of computational time and solution quality. Although the scenarios with high frequency have large problem size, the results show that the scenarios of low and mixed frequencies are more difficult to solve than those of high frequency. Solving the scenarios of low frequency is difficult due to the high headway and insufficient capacity, which can be verified from the high numbers of unserved OD groups in these scenarios. The bus frequency and capacity are higher in the scenarios of mixed frequency than those in the low-frequency scenarios. However, solving the scenarios of mixed bus frequency is still difficult because the bus frequencies of the connecting routes are different, thus complicating the 256

136453.39 89130.26

10 89183.79 60179.58 32.52% 9 685 3600.0 90464.79 55152.66 39.03% 11 293 3600.0 89595.2 0.45% 9 199.0

60227.99 55569.40

7.73%

8

184

3600.0

59708.95 49068.83

257

17.81%

7

1025

3600.0

62213.0 3.19%

6

47.7

96.5

18

135160.8 0.94%

3600.0

2879

21

41.21%

137233.77 80671.86

3600.0

1865

17

34.67%

12x10 1379 0

Number of OD groups Number of passengers Methodology Min. walking time (min) Original (MILP) Objective (Upper/Lower bounds) Optimality gap Unserved OD groups CPU time to find UB (s) Total CPU time (s) SP (MILP) Objective (Upper/Lower bounds) Optimality gap Unserved OD groups CPU time to find UB (s) Total CPU time (s) Heuristic Objective Gap to UB of MILP Unserved OD groups Total CPU time (s) 54.4

41

194900.6 5.49%

3600.0

2707

37

52.28%

196532.59 93768.34

3600.0

515

43

52.05%

206240.99 98890.09

10

Bi-direction (12 stop pairs)

6x10 711 0

Single direction (6 stop pairs)

Direction

Low frequency

Table 2 Computational results (small network, with big-M reduction).

93.4

0

37499.0 0.00%

3600.0

1496

0

3.43%

37498.96 36212.60

3600.0

131

0

3.08%

37498.98 36342.20

6x10 711 0

107.2

0

51500.2 0.00%

3600.0

130

0

17.83%

51500.19 42314.26

3600.0

456

1

15.65%

51500.19 43439.63

10

Single direction (6 stop pairs)

High frequency

309.0

0

77042.0 0.83%

3600.0

473

0

14.01%

77041.97 66248.07

3600.0

2427

0

13.80%

77690.00 66966.89

12x10 1379 0

281.2

1

106782.4 3.68%

3600.0

2120

2

27.72%

108196.20 78193.40

3600.0

687

3

31.37%

110870.20 76087.10

10

Bi-direction (12 stop pairs)

271.7

0

49290.9 1.23%

3600.0

1635

0

22.74%

49290.97 38077.92

3600.0

423

0

20.79%

49908.99 39529.42

6x10 711 0

157.5

0

59989.6 6.51%

3600.0

2597

0

22.62%

59989.59 46414.39

3600.0

1927

2

34.63%

64172.19 41948.91

10

Single direction (6 stop pairs)

Mixed frequency

158.0

2

116892.8 2.61%

3600.0

1215

2

42.57%

120984.97 69476.63

3600.0

1381

4

38.25%

113837.3770287.84

12x10 1379 0

113.0

10

163457.8 6.19%

3600.0

2354

5

45.64%

151733.38 82473.00

3600.0

3436

8

50.14%

153336.40 76440.51

10

Bi-direction (12 stop pairs)

J.C. Chu, et al.

Transportation Research Part E 131 (2019) 247–266

Transportation Research Part E 131 (2019) 247–266

J.C. Chu, et al.

Fig. 2. An example of path enumeration tree of the small network (OD group from stop 1 to stop 7 departing at 12:00).

synchronization of transfers. The second important factor to the computational complexity is the direction of passenger demand. The results clearly show that the problems in which the demand flow is unidirectional, that is, problems in which the synchronization of transfers is simple, are significantly easier to solve than the ones with bidirectional demand flows. The third contributing factor to problem complexity is the requirement for transfer walking time at stops. For the relatively realistic problems, in which transfer walking times are required to transfer successfully, the optimality gap and the solution time increase, meaning the problems are difficult to solve. The problem that does not require a transfer walking time has a much lower objective than the one with this requirement. Thus, if the requirement is ignored, the objective will be seriously underestimated. Notably, the scenarios with mixed frequency and bidirectional flows are clearly more difficult to solve than are the others (using the MILP models), according to the optimality gaps and CPU times required in finding the final feasible solutions. However, the performance of the heuristic algorithm is still better in these scenarios than that of the MILP models. The mixed frequency and requirement for transfer walking time impose reduced difficulty on the heuristic algorithm. 4.1.2. Result analysis (small network) In this section, the results for the scenario with mixed bus frequency, bidirectional flow, and transfer walking time requirement are discussed. This scenario is selected because satisfying the passenger transfer demand is difficult; thus, understanding the realization of synchronization is valuable. The complete bus timetables and passenger itineraries are not reported due to space limitation. As an example, the trajectories of the buses and passengers from Stops 7 to 1 are illustrated in Fig. 3. Notably, the solutions and timetable are solved by SP (MILP), which generated the best results among all approaches. The solid lines are the passenger paths, and each line represents an OD group. The value in the square is the number of passengers, and then the succeeding value is the departure time. The dashed lines represent bus movements. First, the effect of maximum acceptable transfer waiting time is analyzed. Two values of maximum acceptable transfer waiting time, 60 min and 30 min, are tested. The former and the latter represent a situation in which the passengers have high and low tolerance for transfer waiting, respectively. OD 7 → 1 is selected to demonstrate the effect of the maximum acceptable transfer waiting time. By definition of transfer synchronization in this study, if the time difference between the arrival bus of Route K and the departure bus of Route S* is between the transfer walking time (10 min) and the transfer walking time plus the maximum transfer waiting time (70 min), then the two buses are synchronized. Fig. 3 indicates that buses can be easily synchronized under this condition. For example, most buses of Route K are synchronized with buses of Route S* and G* at Stops 4 and 5, respectively. Given these opportunities of transfers, all passengers of OD 7 → 1 complete the trips via Stops 7 → 4 → 1, transferring at Stop 4. Next, the maximum acceptable transfer waiting is diminished from 60 min to 30 min. The objective value of the best feasible solution increased from 151,733.38 to 160,453.76 when the maximum acceptable transfer waiting is reduced. The numbers of unserved OD groups and passengers increased from 5 and 16 to 47 and 163, respectively. The results coincided with the expectation 258

Transportation Research Part E 131 (2019) 247–266

J.C. Chu, et al.

Fig. 3. Trajectories of buses and passengers of OD pairs 7–1 (maximum acceptable transfer waiting time = 60 min).

because a 30 min maximum acceptable transfer waiting time is more restricted than 60 min. Fig. 4 shows the trajectories of buses and passengers of OD pairs 7 → 1 considering the 30 min maximum transfer waiting time. The transfer of five of the OD groups is no longer synchronized due to the reduction of the maximum transfer waiting time, and their previously selected paths become infeasible. Rather than considering them unserved, these OD groups change paths to complete their trips. The changes can be categorized into three types. First, because Buses K-3 and K-9 cannot be synchronized with Buses S*-2 and S*-5 at Stop 4, the departure times of Buses K-3 and K-9 are adjusted to early times. Consequently, OD groups departing at 12:00 and 18:00 will respectively board buses K-4 and K-10 upon arrival at the origin stop, which will be respectively synchronized with Buses S*-2 and S*-5 at Stop 4. Essentially, the total travel times of the two OD groups do not increase. The waiting of the OD groups is shifted from the transfer to the origin stop because the maximum acceptable time of initial waiting is larger than that of transfer waiting in this example. Second, Buses K-1 and K-5 cannot be synchronized with Buses S*-1 and S*-3, respectively, at Stop 4. Consequently, the OD groups at 10:00 and 14:00 must take a detour and complete the trips via Stops 7 → 4 → 3 → 5 → 1, transferring at Stop 5. That is, the two groups board the same buses but alight the buses at Stop 5 to make the transfer. The in-vehicle time increases for the two OD groups, thus raising the total travel time. Third, no feasible paths via physical buses can be identified for the group departing at 16:00; therefore,

Fig. 4. Trajectories of buses and passengers of OD pairs 7–1 (maximum acceptable transfer waiting time = 30 min). 259

Transportation Research Part E 131 (2019) 247–266

J.C. Chu, et al.

Table 3 Sensitivity analysis of capacity constraints. Half demand

Original (MILP)

SP (MILP)

Heuristic

Objective (Upper/ Lower bounds) Optimality Gap Unserved OD groups/passengers CPU time to find UB (s) Total CPU time (s) Objective (Upper/ Lower bounds) Optimality Gap Unserved OD groups/passengers CPU time to find UB (s) Total CPU time (s) Objective Gap to UB of MILP Unserved OD groups/passengers Total CPU time (s)

Original demand

Double demand

With capacity constraint

Without capacity constraint

With capacity constraint

Without capacity constraint

With capacity constraint

Without capacity constraint

43951.59 27018.72 38.52% 4/15

43040.39 27075.58 31.71% 2/5

153336.40 76440.51 50.14% 8/93

125609.39 77110.36 38.60% 0/0

313820.40 237008.45 24.47% 43/804

228164.39 127921.34 43.93% 2/34

881

1607

3436

3548

2896

1801

3600.0 43040.39 29710.45 30.97% 2/5

3600.0 43828.58 29300.97 33.14% 4/18

3600.0 151733.38 82473.03 45.64% 5/47

3600.0 134950.99 82777.22 38.66% 0/0

3600.0 312172.13 240982.66 22.80% 43/797

3600.0 210604.54 136764.35 35.06% 4/78

973

3198

2354

2030

2834

2462

3600.0 49761.8 11.67% 9/31

3600.0 44024.6 2.23% 2/13

3600.0 163457.8 6.19% 10/120

3600.0 137393.8 8.57% 3/38

3600.0 298825.4 4.77% 43/781

3600.0 211781.2 7.18% 4/78

79.0

36.1

113.0

172.1

136.5

47.2

the OD group completes the trip via the dummy bus, indicating that the OD group is unserved. These results show that the proposed methodology is capable of exploring the additional possibilities of passenger choices and transfer connections. The optimized timetables provide systematic synchronized transfers to the OD groups. This scenario is an excellent example showing that passengers must change their paths in response to the timetables for increasing transfer synchronization flexibility and thus timetable effectiveness. Second, a sensitivity analysis is conducted to evaluate the effectiveness of the capacity constraints in the model formulation. Two additional levels of passenger demand are tested. In the case of “half demand,” the number of passengers of each OD group is uniformly distributed from 1 to 7, which is approximately half of the original demand. In the case of “high demand,” the number of passengers is uniformly distributed from 16 to 22, which is approximately double of the original demand. Each level of demand is solved with and without the capacity constraints (Eqs. (27)–(32) in Original (MILP) and Eqs. (39)–(42) in SP (MILP) and the heuristic algorithm). The results are summarized in Table 3. The effectiveness of the capacity constraints is evaluated by comparing the results with and without the capacity constraints. The difference between with and without constraints in the “half demand” and “original” cases is small in terms of the objective and numbers of unserved OD groups and passengers. The difference between with and without constraints becomes clear as demand increases. In the “double demand” case with capacity constraints, the number of unserved OD groups and passengers are 43 and 781–804, respectively. When capacity constraints are not imposed, the corresponding values are only 2–4 and 34–78, respectively. These findings clearly indicate that considering bus capacity is critical for achieving realistic results when passenger demand is high. Disregarding capacity constraints underestimates the unserved OD groups and passengers and leads to seriously overoptimized results by an order of approximately 10 in this particular example. 4.2. Mandl’s network This section presents the results obtained by applying our proposed heuristic algorithm to a larger network. This network was first presented by Mandl (1980) and later used in various studies on transit network design and timetabling, such as those of Baaj and Mahmassani (1991) and Chakroborty (2004). In this example, we apply the proposed heuristic algorithm to the aforementioned network and the four round-trip bus routes (eight one-way routes) proposed in Mandl (1980) and shown in Fig. 5. The network has 15 stops; among which, 6 are transfer stops. Although smaller than many urban bus networks, Mandl’s network is similar in size to the networks studied in other research on timetabling problems formulated in mathematical programming (e.g., Ceder et al., 2001, Ibarra–Rojas et al., 2014, Wong et al., 2008). The four round-trip routes reported in Mandl (1980) are adopted. The original headway of 10 min is too short for the proposed methodology, which aims for long headway bus systems. Thus, the headway is increased to 30 min. Accordingly, the total daily demand is reduced to one-third to balance bus supply and travel demand. The number of nonzero OD pairs in the network is 172, as indicated in Table 4. The total number of passengers is 5134, and 1534 of these passengers require transfer to complete their trip. The bus frequency of each route is set to 20 given the planning horizon duration of 600 min and the headway of 30 min. The passengers of each OD pair are evenly distributed among up to 20 OD groups assuming that demand is uniformly distributed within the day. Consequently, the total number of OD groups in the network is 2084. The maximum acceptable initial and transfer waiting time is set 260

Transportation Research Part E 131 (2019) 247–266

J.C. Chu, et al.

1 8

9 2

2 5

3

8

3

4

4

4

15

3

6

10

12

2

2

2

7

8

8

10 5

8

Routes (two-way)

11

Route A Route B Route C Route D

14 5

2

13

Fig. 5. Swiss road network in Mandl (1980). Table 4 Travel demand matrix (derived from Mandl (1980)).

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15

1

2

3

4

5

6

7

8

9

10

11

12

13

14

15

0 133 66 20 26 50 25 25 10 53 10 8 11 0 0

133 0 16 40 6 60 30 30 5 43 6 3 3 1 0

66 16 0 13 20 60 30 30 5 15 6 3 3 1 0

20 40 13 0 16 33 16 16 5 80 13 8 3 1 0

26 6 20 16 0 16 8 8 3 40 6 5 1 0 0

50 60 60 33 16 0 33 33 10 293 20 5 5 3 0

25 30 30 16 8 33 0 16 5 146 11 3 3 1 0

25 30 30 16 8 33 16 0 5 146 11 3 3 1 0

10 5 5 5 3 10 5 5 0 46 6 1 0 0 0

53 43 15 80 40 293 146 146 46 0 200 83 166 66 0

10 6 6 13 6 20 11 11 6 200 0 25 31 5 0

8 3 3 8 5 5 3 3 1 83 25 0 23 0 0

11 3 3 3 1 5 3 3 0 166 31 23 0 15 0

0 1 1 1 0 3 1 1 0 66 5 0 15 0 0

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

Note: The demands that require a transfer are shown in bold font.

to 30 min considering the headway of 30 min. The remaining parameters are set to be the same as those described in Appendix B. A large size of the Mandl’s network will impose excessive demands on the MILP solver, and tests on this network are conducted using only a heuristic procedure. The number of OD groups that will be added to each iteration (α) is set to 60, and the maximum number of stops of which OD groups are resolved in each iteration (B) is set to 2. The problem is solved within 54173 s with an objective value of 175856.2. In particular, the incremental assignment takes 1820 s to obtain an objective of 196408.6. The local search with b = 1 takes 691 s to reduce the objective to 189369.7 (3.6% improvement). The local search with b = 2 takes 51662 s and further reduces the objective to 175856.2 (7.1% improvement). The detailed timetables and passenger paths cannot be reported in this paper due to page limitation. The transfer synchronization of selected OD groups is analyzed as follows. The network has multiple transfer stops as shown in Fig. 5. Three of the four routes (Routes A, B, and C) meet at Stop 6. Two routes meet at Stops 4, 8, 10, 13, and 15. Most of the OD groups have only one choice of transfer stop, and thus, their transfer synchronization exhibits minimal flexibility. For example, Routes A and C meet only at Stop 6; 261

Transportation Research Part E 131 (2019) 247–266

J.C. Chu, et al.

hence, passengers transferring between the two routes must make their transfer at Stop 6, such as OD pairs 1–12 and 12–1. For several OD groups, however, multiple choices of route combinations and transfer stops are available. For example, OD pairs 10–4 and 4–10 have three means of transfer. Transfer can be made between Routes A and B at either Stops 6 or 8. Passengers can also transfer between Routes A and C at Stop 6. Evidently, assuming that these OD groups transfer at a predetermined stop, similar to that in previous studies, considerably restricts the possibility of transfer synchronization. In the optimized bus timetables and path choices, all three possible means of transfer are utilized by OD pairs 10–4 and 4–10, providing another proof that flexibility in changing paths is beneficial for the optimization of transfer synchronization and the efficiency of a timetable. The test results have shown that optimized timetables are reasonable and provide insights into transfer synchronization planning. The results also show that the proposed heuristic algorithm can still solve the problem within a reasonable time even when the magnitude of the problem is significantly increased, e.g., the number of one-way routes increased from 6 to 8, the number of transfer stops increased from 2 to 6, and the number of OD groups increased from 120 to 2084. 5. Conclusions This paper has presented a model for planning bus timetables with detailed attention to the synchronization of transfers between routes and satisfaction of capacity constraints. The model has been formulated initially as a MILP and then reformulated in setpartitioning terms, and has been solved using big-M reductions for improved problem-solving efficiency. Alongside these exact techniques, a heuristic procedure has been developed which can achieve near-optimal results with much-reduced solution times. The computational test results show that bus frequency, demand direction, and transfer walking time affect bus system performance. The comparison of the test results using MILP and the proposed heuristic procedure on a smaller network shows that heuristic can be relied upon to produce near-optimal solutions, at least at this scale. A detailed examination of the test results also confirms the plausibility of the fine-grained modeling of passengers’ use of bus timetables. Two future research directions are identified as follows. First, the transit assignment model might be extended so as to improve realism with respect to the assumptions discussed in Section 1 above (e.g. to allow for the effects of crowding). Second, testing on actual bus networks could be carried out to investigate the scalability of the proposed methodology and to help establish the credibility of the model as a potential decision-making tool in practice. Declaration of Competing Interest The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper. Acknowledgments The authors are supported by the Minister of Science and Technology of Taiwan through research grants 105-2628-E-002-004MY3, 106-2119-M-002-018, 107-2119-M-002-044 and 108-2628-E-002-003-MY3 and by National Taiwan University through research grants NTUCDP-108L7710 and NTUCDP-108L7751. Appendix A. Big-M reduction The constraints, namely, Eqs. (13), (14), (20)–(23), (28), and (32), and the reformulated Eqs. (33)–(38), (40), and (42), are frequently called big-M constraints. Mathematically, the big-M parameters in these constraints can be set to arbitrarily large values. However, such condition leads to weak linear programming (LP) relaxation and dual bounds and make the problem computationally inefficient and difficult to solve using MILP solvers. The reduction of big-M parameters during problem formulation can considerably strengthen the LP relaxation and tighten the dual bounds of the solution. Consequently, the solution time of MILP with big-M constraints can be remarkably shortened. The big-M reduction technique has been widely used in MILP problem solving, such as in the studies of Dietrich and Escudero (1994), Fügenschuh (2009), and Fischetti and Monaci (2017). For the numerical examples in the current study, a naïve choice for the M-parameters is the length of the planning horizon (denoted by T). As shown in the succeeding analysis, the solution of the model remains inefficient given this parameter. The arbitrarily large or naïve values in big-M constraints can be replaced with tighter values without changing the original problem formulation by carefully estimating the bounds of the variables in such constraints. In the following description, the bounds are determined first, and the calculation of M-parameters is subsequently introduced. Let rllb and rlub be the lower and upper bounds, respectively, of the running times of route l from the first stop to the final one. The two values are the sums of the minimum and maximum travel times (rijlb and rijub ), respectively, that route l traverses and the dwell times at stops except for the first stop. Let rillb and rilub be the lower and upper bounds of the running times of route l from the first stop to stop i. These values are the sums of rijlb and rijub , respectively, that route l traverses from the first stop to stop s and the dwell times at all but

262

Transportation Research Part E 131 (2019) 247–266

J.C. Chu, et al.

Table A1 Summary of big-M reduction. Variable

Lower bound

Upper bound

B Dsln

B Dsln =hllb (n

P Dgs

P Dgs =tg

+ tog s

d (general)

P Dgs =tg

P Ags

P Ags =tg + tog s

d (general)

Wgs

Wgs = 0

P Ags =tg + pg (general)

Sgs

Sgs = 0

Equations

Big-M parameter

rllb

hllb (fl

n) + rslub

+ pg (general)

Wgs =w

Sgs =q

(13)(33)

Sgs

(14)(34)

(

gs

+ Wgs )

(20)(35)

P Dgs

(21)(36)

DsBln

DsBln

P Dgs

(22)(37)

P P Ags ( Rij ln + Dgs )

(23)(38)

P (Rij ln + Dgs )

(

B Dsln =T

1) + rsllb

gs

+ Wgs )

(28)(40)

(DoBg l m

(32)(42)

(DsBl m + 1)

+ 1)

Sgs

P A gs

P Ago g P ( Ags +

gs )

B B B the first stop. Let Dsln and Dsln be the lower and upper bounds of Dsln (the departure time of the nth bus dispatch on route l at stop s), B B B n) + rslub and Dsln respectively. Their bounds are Dsln = T rllb hllb (fl = hllb (n 1) + rsllb . In Dsln , the first two terms indicate the latest possible bus departure time of route l at the first stop. Subtracting the third term from the first two gives the latest possible departure time of the nth bus dispatch of route l at the first stop. Adding the fourth term to the first three gives the latest possible B departure time of the nth bus dispatch of route l at stop s. In Dsln , the first term is the earliest possible departure time of the nth bus dispatch of route l at the first stop. Adding the second term to the first one gives the earliest possible departure time of the nth bus dispatch of route l at stop s. First, we calculate the bounds of passenger travel times between stops to calculate the bounds of the departure and arrival times at the stops for passengers. A network is developed wherein arcs are the union of the physical links of all bus lines (i.e., solid lines in Fig. 1). The travel time of the arc is the lower bound of the bus running time of the arc and dwell time. The passenger travel time from stop i to j is the travel time of the shortest travel time path from stop i to j in the network, which is denoted by tij . Let DgsP and DgsP be the

lower and upper bounds, respectively, of the departure time of OD group g at stop s (DgsP ). Let AgsP and AgsP be the lower and upper

bounds, respectively, of the arrival time of OD group g at stop s ( AgsP ). The general bounds of the departure and arrival times at stop s

are AgsP = tg + tog s

d , DgsP = tg + tog s

d , AgsP = tg + pg , and DgsP = tg + pg . AgsP and DgsP are obtained by the assumptions that the

initial waiting time is zero, the bus running times take the lower bounds, and the time spent at each of the stops before s is d. AgsP and

DgsP are tg + pg because pg is the largest possible travel time for the OD group. The bounds should be modified on the basis of the following three cases. First, if s is the sink node, then AgsP and DgsP are replaced by tg + t og dg d and tg + t og dg d , respectively. The approach is conducted because the sink node is excluded from the aforementioned network for calculating tij . Thus, tog s , where s is the sink node, does not exist and can be replaced by tog dg given that the bus running time from node dg to the sink node is zero. Second, if s is an origin stop, then an OD group always arrives at time tg and can wait for q time units at most. Thus, in the case that s is an origin stop, the upper bounds can be simplified as AgsP = tg and DgsP = tg + q . Third, in the case that s is the dummy destination stop, the OD group is served by the dummy bus, and the arrival time is determined. Thus, if s is a dummy destination stop, two additional lower bounds can be added: AgsP = tg + pg and DgsP = tg + pg , which make AgsP and DgsP constant. Finally, the upper bounds of Sgs and Wgs (i.e., the transfer waiting and total times of OD group g at stop s) are explicitly set in Constraints (17) and (18). The lower bounds are both zero. The M-parameters in the constraints can be reduced on the basis of the above bounds, as summarized in Table A1. The effectiveness of the big-M reduction introduced above is evaluated numerically. Table A2 shows the computational results obtained by the three solution approaches without reducing big-M, and Table 2 illustrates the outcome acquired with big-M reduction. According to the comparison results of the two tables, the big-M reduction is highly beneficial for tightening the lower bounds of the solutions and optimality gaps. In many of the scenarios (for example, the mixed frequency, bi-direction, and 10 min transfer time), the lower bounds remain negative and their optimality gaps exceed 100% after 1 h when the big-M reduction is not adopted in the MILP models. By contrast, the lower bounds and optimality gaps are greatly reduced when the big-M reduction is used.

263

264

89009.52 −48787.35

154.80%

10

3387

3600.0

93879.2

59838.71 −51692.01

186.37%

8

2202

3600.0

61811.0

515.0

3600.0

3600.0

188.6

725

442

16

13

8

6

54.00%

41.45%

6.37%

87895.99 40431.20

59838.99 35033.69

3.19%

10

6 × 10 711 0

Number of OD groups Number of passengers MethodMin. olowalking gy time (min) Original Objective (MI(Upper/ LP) Lower bounds) Optimality gap Unserved OD groups CPU time to find UB (s) Total CPU time (s) Set Objective Part(Upper/ itioLower bounds) ning Optimality (MIgap LP) Unserved OD groups CPU time to find UB (s) Total CPU time (s) Heuristic Objective

Gap to UB of MILP Unserved OD groups Total CPU time (s)

Single direction (6 stop pairs)

Direction

Low frequency

411.5

18

5.08%

134997.8

3600.0

2499

27

213.1

41

1.98%

194900.3

3600.0

2237

47

196.0%

831.0

0

0.00%

37035.0

3600.0

3341

0

492.70%

37034.75 −145450.16

207412.14 −199155.95

154585.29 −201812.38 230.55%

3600.0

1719

0

30.29%

37034.99 25816.80

3600.0

723

45

85.83%

198848.14 28164.39

10

6 × 10 711 0

2211.3

0

11.22%

51500.2

3600.0

286

0

1097.82

0

23.00%

77428.6

3600.0

2280

2

448.37%

1574.7

1

2.86%

105977.5

3600.0

1616

2

355.23%

3205.6

0

10.82%

48827.0

3600.0

3288

1

261.27%

900.8

0

0.00%

59989.6

3600.0

2609

4

197.21%

67913.54 −66024.13

55886.76 −90132.28 123804.79 −315998.78

93431.44 −325494.07

52159.80 −130842.78 350.84%

3600.0

1919

3

63.75%

59989.59 21744.85

3600.0

3287

2

62.35%

54752.98 20614.15

10

3600.0

3407

6

81.42%

119404.60 22176.94

10

6 × 10 711 0

Single direction (6 stop pairs)

3600.0

3069

9

79.72%

100561.40 20390.65

12 × 10 1379 0

Bi-direction (12 stop pairs)

Mixed frequency

3600.0

1638

2

56.52%

58011.99 25222.21

10

Single direction (6 stop pairs)

High frequency

3600.0

1482

23

81.75%

142227.55 25838.11

12 × 10 1379 0

Bi-direction (12 stop pairs)

Table A2 Computational results (small network, without big-M reduction).

670.3

1

7.98%

119155.5

3600.0

1968

10

519.1

10

162326.8 1.39%

3600.0

1665

10

157701.91 −271627.15 272.24% 135287.83 −287394.28 312.43%

3600.0

1908

14

160065.59 −5382.74 103.36%

10

3600.0

3529

7

106.29%

129499.00 −8157.75

12 × 10 1379 0

Bi-direction (12 stop pairs)

J.C. Chu, et al.

Transportation Research Part E 131 (2019) 247–266

Transportation Research Part E 131 (2019) 247–266

J.C. Chu, et al.

Appendix B. Parameters of numerical tests The parameters adopted in the computational tests in Section 4 are defined and they are summarized in Table B1. Most of the parameters are decided according to sources. Few of the parameters are given subjectively by the authors. In future studies, additional data can be collected for accurately determining these parameters. In the example, the minimum transfer walking time ( gs ) at every stop is 10 min for all OD groups. The value is calculated based on the walking time of 5 km/h and the distance of 0.8 km (Browning et al. 2006; Yu et al., 2012). The other parameters are as follows. The dwelling time at every stop (d ) is 1 min (Kittelson & Associates, 2003). The penalty value ( pg ) for OD group g is assumed to be the maximum acceptable travel time of the OD group plus 1 min. The maximum number of total transfers k is 1. The value is chosen because the trips can be completed with at most one transfer with particularly few exceptions in this example. The maximum acceptable transfer waiting time w of 60 min is adopted as a design parameter of the timetables. This value is consistent with the maximum transfer waiting time observed in the real world (Ma et al., 2013) and the long headway in the example. The value of this parameter can be adjusted depending on the problem. A substantially smaller value (e.g., 30 min) is also used in an experiment to demonstrate the capability of the proposed methodology in Section 4.1.2. The maximum acceptable initial waiting time of passengers q is the largest average bus headway at the origin stop of the passengers. The initial waiting time cannot be greater than the bus headway at the origin stop of the passengers because passengers can board the first bus when they arrive at the stop. Furthermore, when a timetable is available to the passengers, they can wait at home and only go to the stop when departure time is near. That is, waiting occurs at home. Therefore, initial waiting time has relatively minimal negative effect on passengers, and thus, its maximum acceptable value is not as restricted as that of transfer waiting time. The bus capacity c is 45 persons (Nookala and Khan, 1987), and the bus operation period is 600 min (10 h). This condition implies that the average headway is 60 min for a frequency of 10, and 120 min for a frequency of 5. The bounds of bus headway [hllb , hlub] are the average headway (total time/frequency) minus and plus 5%. The average bus running times are indicated in Figs. 1 and 5. The bounds of the bus running times [rijlb, rijub] are the average running time minus and plus 20%. The total travel time limit of OD group g is set according to ug = max{3 × bog dg , bog dg + q}, g G , where bog dg is the shortest possible travel time between og and dg (i.e., the travel time by a private vehicle). The main idea of the formula is that the acceptable maximum travel time should not be constant for all passengers. Instead, it should be proportional to the length of the trip. Thus, the first criterion of the maximum acceptable travel time is a multiple of the shortest possible travel time between the same OD pair. In this study, three times of the shortest possible travel time is adopted. However, in cases where the headway of the bus service is large but the trip length is short, a multiple of the shortest possible travel time is small compared with the possible waiting time of the trip and might be insufficient for the passengers to complete the trip. Therefore, the second criterion of the maximum acceptable travel time is set to the maximum acceptable initial waiting and shortest possible travel times of the OD pair. The larger of the two criteria is selected.

Table B1 Summary of parameters in the numerical tests. Parameter

Value

Justification

Transfer walking time

10 min

Bus capacity Dwelling time

45 passengers 1 min

Maximum acceptable transfer waiting time

30, 60 min

Maximum acceptable initial waiting time

30, 60, 120 min

Maximum number of transfer

1

According to Chinese standards and good service level, the station spacing is set to 0.5–0.8 km for the skeleton routes and to 0.3–0.5 km for the main and branch routes. We assume that the speed of pedestrians is 5 km/h with a maximum station spacing distance of 800 m. Thus, the transfer walking time should be 9–10 min (Yu et al., 2012). Nookala and Khan (1987) The door opening and closing time is assumed to be 4 s. The boarding and alighting times per passenger are assumed to be 3 and 2 s, respectively. The 1 min dwelling time is assumed on the basis of these values (Kittelson and Associates, 2003). Ma et al. (2013) analyzed the ridership data of Beijing in 2010. They found that the average transit transfer time is 25.4 min. They also determined that over 94% of transfer activities take less than 60 min. From the preceding statistics and the headway of the bus routes in the computational tests, 30 min and 60 min are adopted in the small network, while 30 min is used in Mandl’s network. The value is set to the largest average headway of all the routes in the network. In particular, 60 min is used for high-frequency cases and 120 min is used for low- and mixed-frequency cases in the small network. The value is 30 min in Mandl’s network. In the examples, the trips can be completed with at most one transfer with particularly few exceptions.

Appendix C. Supplementary material Supplementary data to this article can be found online at https://doi.org/10.1016/j.tre.2019.10.001.

265

Transportation Research Part E 131 (2019) 247–266

J.C. Chu, et al.

References Baaj, M., Mahmassani, H., 1991. An AI-based approach for transit route system planning and design. J. Adv. Transportation 25 (2), 187–209. Balinski, M., Quandt, R., 1964. On an integer program for a delivery problem. Oper. Res. 12 (2), 300–304. Browning, R., Baker, E., Herron, J., Kram, R., 2006. Effects of obesity and sex on the energetic cost and preferred speed of walking. J. Appl. Physiol. 100 (2), 390–398. Carraresi, P., Malucelli, F., Pallottino, S., 1996. Regional mass transit assignment with resource constraints. Transportation Res. Part B: Methodological 30 (2), 81–98. Ceder, A., 2007. Public Transit Planning and Operation: Theory, Modeling and Practice. Taylor & Francis, Oxford. Ceder, A., Golany, B., Tal, O., 2001. Creating bus timetables with maximal synchronization. Transportation Res. Part A: Policy Practice 35 (10), 913–928. Ceder, A., Hassold, S., Dano, B., 2013. Approaching even-load and even-headway transit timetables using different bus sizes. Public Transport 5 (3), 193–217. Chakroborty, P., 2003. Genetic algorithms for optimal urban transit network design. Comput.-Aided Civ. Infrastruct. Eng. 18 (3), 184–200. Chen, J., Liu, Z., Zhu, S., Wang, W., 2015. Design of limited-stop bus service with capacity constraint and stochastic travel time. Transportation Res. Part E: Logistics Transportation Rev. 83, 1–15. Desaulniers, G., Hickman, M., 2007. Handbooks in Operations Research & Management Science - Public Transit. Elsevier, Burlington, pp. 69–127. Dietrich, B.L., Escudero, L.F., 1994. Obtaining clique, cover and coefficient reduction inequalities as Chvatal-Gomory inequalities and Gomory fractional cuts. Eur. J. Oper. Res. 73 (3), 539–546. Fischetti, M., Monaci, M., 2017. Using a general-purpose mixed-integer linear programming solver for the practical solution of real-time train rescheduling. Eur. J. Oper. Res. 263 (1), 258–264. Fügenschuh, A., 2009. Solving a school bus scheduling problem with integer programming. Eur. J. Oper. Res. 193 (3), 867–884. Guihaire, V., Hao, J., 2008. Transit network design and scheduling: A global review. Transportation Res. Part A: Policy Pract. 42 (10), 1251–1273. Guo, X., Sun, H., Wu, J., Jin, J., Zhou, J., Gao, Z., 2017. Multiperiod-based timetable optimization for metro transit networks. Transportation Res. Part B: Methodological 96, 46–67. Guo, X., Wu, J., Zhou, J., Yang, X., Wu, D., Gao, Z., 2018. First-train timing synchronisation using multi-objective optimisation in urban transit networks. Int. J. Prod. Res. 1–16. Hewitt, M., Nemhauser, G., Savelsbergh, M., 2010. Combining exact and heuristic approaches for the capacitated fixed-charge network flow problem. INFORMS J. Computing 22 (2), 314–325. Hewitt, M., Nemhauser, G., Savelsbergh, M., 2013. Branch-and-price guided search for integer programs with an application to the multicommodity fixed-charge network flow problem. INFORMS J. Computing 25 (2), 302–316. Ibarra-Rojas, O., Rios-Solis, Y., 2012. Synchronization of bus timetabling. Transportation Res. Part B: Methodological 46 (5), 599–614. Ibarra-Rojas, O.J., Delgado, F., Giesen, R., Muñoz, J.C., 2015. Planning, operation, and control of bus transport systems: A literature review. Transportation Res. Part B: Methodological 77, 38–75. https://doi.org/10.1016/j.trb.2015.03.002. Ibarra-Rojas, O., Giesen, R., Rios-Solis, Y., 2014. An integrated approach for timetabling and vehicle scheduling problems to analyze the trade-off between level of service and operating costs of transit networks. Transportation Res. Part B: Methodological 70, 35–46. Kang, L., Wu, J., Sun, H., Zhu, X., Wang, B., 2015. A practical model for last train rescheduling with train delay in urban railway transit networks. Omega 50, 29–42. Kang, L., Zhu, X., Sun, H., Wu, J., Gao, Z., Hu, B., 2019. Last train timetabling optimization and bus bridging service management in urban railway transit networks. Omega 84, 31–44. Kittelson and Associates, 2003. Transit Capacity and Quality of Service Manual. Transportation Research Board. Kwan, C., Chang, C., 2008. Timetable synchronization of mass rapid transit system using multiobjective evolutionary approach. IEEE Transactions Syst., Man, Cybernetics, Part C (Appl. Rev.) 38 (5), 636–648. LeBlanc, L.J., 1975. An algorithm for the discrete network design problem. Transportation Sci. 9 (3), 183–199. Ma, X., Wu, Y.J., Wang, Y., Chen, F., Liu, J., 2013. Mining smart card data for transit riders’ travel patterns. Transportation Res. Part C: Emerging Technol. 36, 1–12. Magnanti, T., Wong, R., 1984. Network design and transportation planning: models and algorithms. Transportation Sci. 18 (1), 1–55. Mandl, C., 1980. Applied Network Optimization. Academic Press, London. Marcotte, P., 1983. Network optimization with continuous control parameters. Transportation Sci. 17 (2), 181–197. Migdalas, A., 1995. Bilevel programming in traffic planning: Models, methods and challenge. J. Global Optim. 7 (4), 381–405. Nookala, M., Khan, A., 1987. Cost-efficiency of intercity bus technology innovations. Transp. Res. Rec. 1125, 57–63. Schroder, M., Solchenbach, I., 2006. Optimization of Transfer Quality in Regional Public Transit. Technical report, Fraunhofer ITWM, Kaiserslautern. Sun, H., Wu, J., Ma, H., Yang, X., Gao, Z., 2019. A bi-objective timetable optimization model for urban rail transit based on the time-dependent passenger volume. IEEE Trans. Intell. Transp. Syst. 20 (2), 604–615. Tilahun, S., Ong, H., 2012. Bus timetabling as a fuzzy multiobjective optimization problem using preference-based genetic algorithm. PROMET - Traffic & Transportation 24 (3), 183–191. Toth, P., Vigo, D., 2015. The Vehicle Routing Problem. Society for Industrial and Applied Mathematics, Philadelphia. Wang, H., Yang, X., Wu, J., Sun, H., Gao, Z., 2018. Metro timetable optimisation for minimising carbon emission and passenger time: a bi-objective integer programming approach. IET Intel. Transport Syst. 12 (7), 673–681. Wang, Y., Tang, T., Ning, B., Meng, L., 2017. Integrated optimization of regular train schedule and train circulation plan for urban rail transit lines. Transportation Res. Part E: Logistics Transportation Rev. 105, 83–104. Wardrop, J.G., 1952. Road paper. Some theoretical aspects of road traffic research. Proc. Inst. Civil Engineers 1 (3), 325–362. Wong, R., Yuen, T., Fung, K., Leung, J., 2008. Optimizing timetable synchronization for rail mass transit. Transportation Sci. 42 (1), 57–69. Wu, J., Liu, M., Sun, H., Li, T., Gao, Z., Wang, D., 2015. Equity-based timetable synchronization optimization in urban subway network. Transportation Res. Part C: Emerging Technol. 51, 1–18. Xiong, J., He, Z., Guan, W., Ran, B., 2015. Optimal timetable development for community shuttle network with metro stations. Transportation Res. Part C: Emerging Technol. 60, 540–565. Yang, H., Bell, M., 1998. Models and algorithms for road network design: a review and some new developments. Transport Rev. 18 (3), 257–278. Yang, X., Chen, A., Ning, B., Tang, T., 2017. Bi-objective programming approach for solving the metro timetable optimization problem with dwell time uncertainty. Transportation Res. Part E: Logistics Transportation Rev. 97, 22–37. Yin, H., Wu, J., Sun, H., Kang, L., Liu, R., 2018. Optimizing last trains timetable in the urban rail network: social welfare and synchronization. Transportmetrica B: Transport Dyn. 1–25. Yu, B., Yang, Z., Jin, P., Wu, S., Yao, B., 2012. Transit route network design-maximizing direct and transfer demand density. Transportation Res. Part C: Emerging Technol. 22, 58–75. Zhang, W., Xu, W., 2017. Simulation-based robust optimization for the schedule of single-direction bus transit route: The design of experiment. Transportation Res. Part E: Logistics Transportation Rev. 106, 203–230.

266