Synchronization of bus timetabling

Transportation Research Part B 46 (2012) 599–614

Contents lists available at SciVerse ScienceDirect

Transportation Research Part B journal homepage: www.elsevier.com/locate/trb

Synchronization of bus timetabling Omar J. Ibarra-Rojas ⇑, Yasmin A. Rios-Solis Graduate Program in Systems Engineering, Universidad Autónoma de Nuevo León (UANL), Mexico

a r t i c l e

i n f o

Article history: Received 4 November 2010 Received in revised form 13 January 2012 Accepted 13 January 2012

Keywords: Bus timetabling Synchronization Passenger transfer Bus bunching NP-hard Feasible solution space

a b s t r a c t Timetable generation is a subproblem of bus network strategic planning, in which the departure time of each trip is determined. We study the bus network of Monterrey, Mexico, which is similar to those of other cities in Latin America. It is a large bus network where passenger transfers must be favored, almost evenly spaced departures are sought, and bus bunching of different lines must be avoided. We formulate the timetabling problem of this network with the objective of maximizing the number of synchronizations to facilitate passenger transfers and avoid bus bunching along the network. We define these synchronizations as the arrivals of two trips with a separation time within a time window to make a flexible formulation. This flexibility is a critical aspect for the bus network, since travel times vary because of reasons such as driver speed, traffic congestion, and accidents. By proving that our problem is NP-hard we answer a 10-year-old open question about the NP-hardness of similar problems present in literature. Next, we analyze the structural properties of the feasible solution space of our model. This analysis leads to a preprocessing stage that eliminates numerous decision variables and constraints. Moreover, this preprocessing defines feasible synchronization and arrival time windows that are used in a new metaheuristic algorithm. Empirical experimentation shows that our proposed algorithm obtains high-quality solutions for real-size instances in less than one minute.  2012 Elsevier Ltd. All rights reserved.

1. Introduction The planning process of a bus network is complicated. Therefore, it is often divided into several subproblems such as line planning, timetable generation, vehicle scheduling, and crew scheduling (Ceder, 2007). In the timetabling problem we determine departure times for trips from different lines to achieve a specific goal, which in our case is to maximize the number of synchronizations, to favor passenger transfers and avoid bus bunching. Preparing a timetable is a critical step in the planning process of a bus network since service quality and subsequent planning steps such as vehicle and crew scheduling depend on its solution. We focus on Monterrey, Mexico, where the bus transit system is private. Moreover, there are different bus companies sharing the demand. Competition among them creates a particularly large bus network. The left side of Fig. 1 shows the entire bus network of Monterrey’s metropolitan area. Our case has the following characteristics:  The bus network is private and has more than 300 bus lines (Monterrey has a population close to the 4 million).  There is almost full coverage of the territory. The area marked with a black box on the left part of Fig. 1 shows the line concentration downtown since there is a tendency to design lines crossing this area. Therefore, many transfers are located there. ⇑ Corresponding author. Address: Av. Universidad s/n, 66450 Monterrey, Nuevo León, Mexico. Tel.: +52 81 1492 0363. E-mail addresses: omar@yalma.fime.uanl.mx (O.J. Ibarra-Rojas), yasmin@yalma.fime.uanl.mx (Y.A. Rios-Solis). 0191-2615/$ - see front matter  2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.trb.2012.01.006

600

O.J. Ibarra-Rojas, Y.A. Rios-Solis / Transportation Research Part B 46 (2012) 599–614

Fig. 1. The left panel shows the bus network of Monterrey’s metropolitan area. The marked area shows the line concentration downtown (Rabut, 2010). The right panel shows a bus line i consisting of three sub-lines a, b, and c.

 Many lines have ‘‘sub-lines’’ which share a common route segment but have variations at the beginning and at the end of the routes (see right part of Fig. 1).  The timetable is only for the company’s administration. Passengers do not know the time at which bus arrives at each stop. They only have an estimate of their waiting time to take the next bus. Considering the bus network with the above-mentioned characteristics, what can we optimize? To answer this question, we noted the need to provide quality service and meet other company priorities. We found the following key components to optimize through timetable generation.  Passenger transfers: Travel from one point to another might imply transfer between lines. We seek short waiting times at some key nodes like downtown.  Bus bunching: Due to the presence of sub-lines, different lines commonly converge on a specific node (called a bus bunching node). We aim to avoid bus bunching between sub-lines or between different lines using trip separation.  Almost evenly spaced departures: A large variation in the time between consecutive trips affects the behavior of passenger demand, even in small planning periods. Because of the uncertainties inherent in a public transportation system such as Monterrey’s, passengers prefer a flexible transfer rather than an instant transfer. For example, different travel times due to variations in driver speeds are drawbacks for achieving punctual synchronization. In our case of study, a flexible model that results in a more representative formulation is required. Therefore, we consider two trips to be synchronized if the difference between arrival times at a transfer or a bus bunching node is within a given time window. Passenger activity during the day also affects the bus network behavior. Peak hours imply a large amount of bus bunching of different lines on multiple nodes of the network and a large number of passenger transfers. Therefore, we categorize the days (weekend, Monday, holidays, etc.) and then, we split each day into planning periods (peak hours, morning, afternoon, etc.) to achieve more accurate deterministic travel times. The characteristics of the Monterrey case are also present in other bus networks in Latin America, and thus any flexible model that optimizes the synchronization between lines is more relevant and useful. Maximizing line synchronization to reduce bus bunching between different lines and optimize passenger transfers is an important objective for both; companies and passengers. Efficient timetabling increases the quality of service and might bring several concomitant benefits such as passengers switching from cars to public buses, and hence a reduction in the number of vehicles in the area. On the basis of the above, we propose a formulation for the synchronization bus timetabling problem to reduce bus bunching and to optimize passenger transfer. Our formulation is based on Ceder et al. (2001), which was later improved by Eranki (2004). The objective is to maximize the number of synchronizations along the network during a planning period. The enhancements we introduce to the integer linear programming are the following: departure times distribution along the planning period (in Eranki’s case, one can see that departures might be concentrated at the beginning of the planning period), oriented synchronization (transfer form line i to line j, but not vice versa), and almost evenly spaced departure times to avoid variations in passenger demands. These improvements constitute a new mathematical model that is more suitable to analyze the Monterrey case. The remainder of this paper is organized as follows. First, we present a literature review in Section 2. Our improved formulation is presented in Section 3. In Section 4, we present the idea for the complexity proof of our timetabling problem (the details are presented in Appendix A). This proof can be used for similar formulations such as Ceder et al. (2001) and Eranki (2004), therefore, it answers a 10-year-old open question about the complexity of these problems. In Section 5, we present some characteristics of the feasible solution space of our problem which generates a preprocessing stage that leads to a

O.J. Ibarra-Rojas, Y.A. Rios-Solis / Transportation Research Part B 46 (2012) 599–614

601

smaller formulation. The convergence of the branch and bound (B&B) algorithm is improved since we eliminate a large number of constraints and decision variables. In Section 6, we introduce a new metaheuristic based on the preprocessing stage that help us solve real-size instances as can be seen in Section 7. Finally, Section 8 summarizes our results and suggests areas of further research. 2. Literature review The strategic planning of transit networks (network design, setting frequencies, and generating timetables) has been addressed in several studies, as it can be seen in survey of Guihaire and Hao (2008a). There are different formulations for every subproblem that corresponds to particular characteristics of the transit network. The timetabling problem itself has different formulations that factor in: frequencies, evenly spaced departure times, passenger activity-based planning, and passenger transfer at one or multiple stations, among others. Passenger transfer is crucial in public transit systems. In several studies (Bookbinder and Desilets, 1992; Cevallos and Zhao, 2006; Chakroborty et al., 1997; Chakroborty et al., 1995; Daduna and Voß, 1995; Jansen and Nielsen, 2002; Wong et al., 2008) the authors seek to minimize the overall transfer waiting times and consider evenly spaced departure times. The minimization of overall waiting time could not be a representative objective in our study since we also address bus bunching. Schröder and Solchenbach (2006) remark on the necessity of better measurement of the transfer function, since minimizing waiting times may lead to risky passenger transfers due common delays in bus arrivals at stops. They make a formulation to minimize the cost of different kinds of transfers using time windows. Ceder and Tal (2001) address the problem of synchronizing arrivals of different lines at some nodes to achieve more efficient passenger transfer. They define synchronization as the simultaneous arrival of two buses and consider limits for consecutive trips’ separation times (referred to as headway times) and frequency (number of trips for each line) as given. The objective of the problem is to maximize the number of simultaneous arrivals. The authors note that their problem seems intractable since solving small instances lead to large computational times. They design constructive heuristic procedures to generate timetables. In Ceder et al. (2001), the mathematical formulation of the previous problem is presented. In addition, they show the heuristic implementation on a real case study, where they found good solutions within seconds. One extension of Ceder et al. (2001) is presented in Eranki (2004), where synchronization is redefined as the arrival of two trips at a synchronization node with a separation time within a small time window instead of simultaneous arrivals. Eranki designs a constructive heuristic algorithm on the basis of Ceder’s algorithm to solve her formulation. We based our model on this formulation, and we enhanced several of its components like redefining the variables to have oriented synchronization (passengers switch from line i to line j, but not necessarily vice versa), taking into account almost evenly spaced departures, and avoiding bus bunching between sub-lines. These enhancements make our model more realistic. Here, we prove that indeed the problems presented in Eranki (2004) and Ceder et al. (2001) are NP-hard. Liu et al. (2007) address the synchronization timetabling problem based on the formulation presented in Ceder et al. (2001). They define an objective function for simultaneous arrivals. The function consists of a ratio of the number of routes where there are vehicles arriving simultaneously at a connection stop to the number of all routes passed at the same stop. Hall et al. (2001) address a model for holding times to ensure passenger transfers. Note that we do not factor in this last point of view in our problem analysis, since bus holding at transfer nodes could imply bus bunching of sub-lines. Nevertheless, considering bus holding times would be an interesting research area. There are also formulations that consider synchronization between different transport systems. Clever (1997) considers transfers in intermodal transportation systems and implements minimum headway times to avoid loss of passenger transfer. Chowdhury and Chien (2001) also focus on the intermodal case, and consider fixed demand. They propose a formulation to minimize the operational cost. Bruno et al. (2009) address the transfer optimization problem at only one key node for the case of an intermodal transit system. They present two formulations to minimize the operating cost of activating a trip and waiting. They also consider an initial timetable as given. The bus bunching problem has been addressed by Daganzo (2009), where trips of some lines are not separated, and thus they often arrive as a bunch. He develops a dynamic method based on real-time data to generate holding times for bus runs to avoid bunching. Adamski (1993) addresses the synchronization problem of lines that share common route segments. He presents a mathematical formulation to harmonize headway times of different lines to separate trips in a planning period. He considers possible departure times as given and solves the problem with an integrated tabu search and genetic approach using small computational resources. This is different from our case study since it might be difficult to harmonize headway times for bus lines of different companies. In addition, he considers evenly spaced departures, i.e., fixed headway times for each line. We are interested in optimizing passenger transfer and reducing congestion or bus bunching with almost evenly spaced departure times. Strategic planning is also based on passenger behavior (not constant demand) and real-time (travel time) data. Li et al. (2010) study the transit scheduling problem to optimize the interaction of different services at an intermodal transport network considering demand variability and therefore variable headways. They use a heuristic solution algorithm that combines the Hooke–Jeeves method and an iterative supply-demand equilibrium approach to solve their model. Dessouky et al. (1999) address the minimization of transfer delays: they make a probabilistic model based on real-time data for delays, and then determine the departure times to achieve their objective. Chung (1990) addresses a travel time forecasting process in the bus network for peak and non-peak hours. On the one hand, the author solves the simultaneous synchronization

602

O.J. Ibarra-Rojas, Y.A. Rios-Solis / Transportation Research Part B 46 (2012) 599–614

problem with the objective of minimizing overall waiting time using genetic algorithms. On the other hand, the author presents a real-time control strategy to maintain transfers against unexpected delays. Integration of subproblems is an important research area in transit network planning. Van den Heuvel et al. (2008) address the integration of the clock-face timetabling problem (trips depart at regular intervals, and thus at the same number of minutes past each hour) and multiple-depot multiple-vehicle-type scheduling problem. They design a tabu search to solve their proposed integration. Guihaire and Hao (2008b) also implement a similar idea. Their timetabling problem consists of minimizing the overall waiting time considering a given frequency and fixed headway times. They design an iterated local search for solving their formulation. Fleurent and Lessard (2009) propose a measure function for transfers based on ideal waiting times. They design an optimization approach to minimize vehicle costs such as the number of vehicles and unproductive time. Fleurent et al. (2005) propose an integral formulation for timetabling and vehicle scheduling problems that considers weights on the objective function. However, these weights reflect the planner’s necessity, which is a complex characteristic if two or more objectives are in conflict. The literature review shows different characteristics of the timetabling problem. However, to the best of our knowledge, there is no formulation that represents all the characteristic of our case study, which is similar to the case of some cities in Latin America, such as Mexico City. 3. Integer programming formulation In this section, we present our integer programming formulation for the Bus Synchronization Timetabling Problem (BTP). The bus network is represented by a set of nodes and a set of lines. There are only a few specific stops for the bus lines because a bus can stop at almost every corner along its route. Therefore, in our case, a node does not represent a single stop but the beginning of a shared route section or a route section where transfers occur. We consider the existence of two types of key nodes: transfer and bus bunching nodes. Fig. 2 shows these two cases. Case (a) represents two lines, i and j, which converge at node b and have a common line segment. Avoiding bus bunching of these lines implies separation of their trips on node b. This represents the case of sub-lines, where for example, the separation time between the arrivals of two buses is 5 min to avoid bunching. Case (b) represents a transfer node where many passengers would like to go from trips on line i to that on line j. A small waiting time, for example, 3 min is needed for passenger transfers. In this study, the events of avoiding bus bunching of different lines and favoring passenger transfers are called synchronization events. This study is based on real data obtained from a management company of an urban transit system that has a large data-base of passenger demand for bus lines. Therefore, it is possible to divide any day into several planning periods on the basis of demand and travel time behavior to obtain accurate deterministic data. For the formulation, we use the following notations. The set of lines of the bus network is denoted by I. For each line i 2 I, we define J(i) as the set of lines that may have synchronization nodes with line i. Set Bij represents all the synchronization nodes for lines i and j 2 J(i). T is the planning period in minutes. The number of trips needed to fulfill the demand (frequency) of each line i 2 I is represented by fi. Travel times from depot of each line i to each synchronization node b are represented by tib . The interval [wb;Wb] is the waiting time window for each transfer or bus bunching node b. Usually, within a bus company, consecutive trips have a separation time (called headway times) between a minimum hi and maximum Hi time. To achieve almost evenly spaced departure times, we define an ‘‘ideal headway’’ time (evenly separates consecutive trips) given by i Hideal ¼ fTi . On the basis of this parameter, we define the minimum and maximum headway times as hi = Hideali  di and i H = Hideali + di, respectively, where di is a given headway amplitude factor. The decisions are represented by two types of variables. The first type is to determine the departure time for every trip of each line, and the second one identifies the existence of a synchronization between two trips of two different lines. These variables are represented as follows.  X ip : (integer) departure time of the p-th trip of line i.  1 if trip p of line i arrives first at node b and synchronizes with trip q of line j;  Y ijpqb ¼ 0 otherwise:

i

i 14

14

j

b 19 Wb

17

b Wb

5

3

j

(a)

(b)

Fig. 2. Two types of synchronization nodes on the bus network. Case (a) represents a bunching node and Case (b) represents a transfer node.

603

O.J. Ibarra-Rojas, Y.A. Rios-Solis / Transportation Research Part B 46 (2012) 599–614

Then, BTP is given by

F BTP ðYÞ ¼ max

fi X fj XX XX Y ijpqb i2I j2JðiÞb2Bij p¼1 q¼1

subject to X i1 6 Hi i

T H 6

8i2I X if i

6T

ð1Þ

8i2I

ð2Þ

i

h 6 X ipþ1  X ip 6 Hi 8 i 2 I; p ¼ 1; . . . ; f i  1       8 i 2 I; j 2 JðiÞ; p ¼ 1; . . . ; f i ; X jq þ tjb  X ip þ t ib P wb  M 1  Y ijpqb       8 i 2 I; j 2 JðiÞ; p ¼ 1; . . . ; f i ; X jq þ tjb  X ip þ t ib 6 W b þ M 1  Y ijpqb X ip 2 f0; 1; . . . ; Tg i Y ijpqb 2 f0; 1g

8 i 2 I; p ¼ 1; . . . ; f

i

8 2 I; j 2 JðiÞ; p ¼ 1; . . . ; f i ;

ð3Þ q ¼ 1; . . . ; f j ; b 2 Bij

ð4Þ

q ¼ 1; . . . ; f j ; b 2 Bij

ð5Þ ð6Þ

q ¼ 1; . . . ; f j ; b 2 Bij

ð7Þ

The objective function maximizes the number of synchronizations. Constraints (1) force the first trip of each line to depart at the beginning of the planning period T. Similarly, constraints (2) force the last trip of each line to depart at the end of the planning period. Constraints (3) separate consecutive trips of each line i by a number of minutes between minimum hi and maximum Hi headway times. Constraints (4) and (5) allow the synchronization variables Y ijpqb be activated if the difference between arrivals of qth trip of line j and p-th trip of line i at node b is between wb and Wb. Although M is a big number, we can bind it by the maximum difference of arrival times for every pair of lines (i, j) that synchronize at every node b, that is, n n  o n oo  M ¼ max maxi;j2JðiÞ;b2Bij T þ t jb  tib ; maxi;j2JðiÞ;b2Bij T þ t ib  t jb . Finally, constraints (6) and (7) represent the domain of the decision variables. The previous model is an enhancement of the model proposed by Eranki (2004). The first modification to adapt Eranki’s model to Monterrey’s case is to redefine the synchronization variables. We consider oriented synchronizations, i.e., pth trip     of line i synchronizes with qth trip of line j at node b if the difference X jq þ t jb  X ip þ tib is between wb and Wb. Therefore, Y ijpqb models transfer from line i to line j and not vice versa. Since we seek trip separation in these type of nodes, we can model synchronization at a bunching node as a double oriented transfer. We also add constraints (2) to achieve departure dispersion along the planning period T. Finally, we consider the characteristics of almost evenly spaced headway times using the headway amplitude factor. 4. Computational complexity of BTP Computational complexity theory states that problems can vary in the effort required to solve them (e.g. Garey and Johnson, 1979). An optimization problem that belongs to the NP-hard class means, in simple terms, that an efficient algorithm for the exact solution of this problem does not exist. Proving that a problem is NP-hard is not a negative result. Rather it is an important result because researchers are unlikely spend time searching for an efficient algorithm but instead will seek approximate solutions (like the one we present in Section 6) or resort to other ways to tackle the problem. Studies of Ceder et al. (2001) and Eranki (2004) remark on the intractability of closely related problems of BTP with the help of empirical results. Nevertheless, the question whether these problems were part of the NP-hard class was a 10-year-old open question. In this paper, we prove that BTP belongs to the NP-hard class by proving that its decision version belongs to NPcomplete class. We propose a polynomial reduction from the monotone NAE-3SAT, whose NP-completeness is assured by Schaefer (1978), to the decision version of BTP. With this proof we also guarantee that indeed problems of Ceder et al. (2001) and Eranki (2004) are also NP-hard. The complexity proof of BTP is not trivial and its details can be found in Appendix A. 5. Preprocessing stage Since many synchronization variables are forced to be zero, the solution space of BTP has an interesting structure. Indeed in our BTP formulation, we define synchronization variables Y ijpqb for every node b 2 Bij, trips p and q of lines i and j 2 I(i), respectively. Nevertheless, it is improbable that the first trip of line i could synchronize with the last trip of line j due their departure time constraints (1) and (2). Fig. 3 shows a typical structure of the matrix formed by synchronization variables for the pair of lines (i,j) at some node b where the element in pth row and qth column represents to variable Y ijpqb . Zeros in the matrix represent the variables forced to be zero. The variables between the curved lines are the only ones whose values could be 1, i.e., they might satisfy the constraints of BTP. To identify these synchronization variables forced to be zero, we deduce an implicit synchronization window for all trips p of each line i. To achieve this goal, we first determine a feasible departure time window for trip p using arrival times at nodes, headway constraints, and waiting time parameters. We denote this departure time window as Dip . Procedure 1 shows the steps to obtain the Dip (this is also illustrated by Fig. 4).

604

O.J. Ibarra-Rojas, Y.A. Rios-Solis / Transportation Research Part B 46 (2012) 599–614

Procedure 1. Generate Dip . 1. Suppose first trip departs at X i1 ¼ 0, calculate the earliest departure time for trip p of line i given by earliest1 = (p  1)hi. 2. Suppose first trip departs at X i1 ¼ Hi , calculate the latest departure time for trip p of line i given by latest1 = min{T, pHi}. i 3. Suppose last trip departs at X fp ¼ T  Hi , calculate the earliest departure time for trip p of line i given by i earliestf i þ1 ¼ maxf0; T  ðf  ðp  1ÞÞHi g. i f i þ1 i 4. Suppose last trip departs at X fp ¼ T, calculate the latest departure time for trip p of line i given by latest ¼ T  ðf i  pÞh . i i 5. Define Dp as ½earliest1 ; latest 1  \ ½earliestf i þ1 ; latest f i þ1 . Therefore, the departure time window Dp for every trip p of each line i can be defined as follows: i

i

Dip ¼ ½maxfðp  1Þh ; T  ðf i  ðp  1ÞÞHi g; minfpHi ; T  ðf i  pÞh g

ð8Þ

Fig. 4 shows an example of the departure time window construction for trip p = 8 of line i. This example has a planning period of T = 30 minutes, a frequency of fi = 10, and headway amplitude factor of di = 33%, i.e., hi = 2 and Hi = 4. Fig. 4 has four time lines. First line shows the earliest departure time 7hi = 14 for the eighth trip, assuming that first trip departs at X i1 ¼ 0 (step 1). Second time line shows the latest departure time min{T,8Hi} = 30 for the eighth trip, assuming that first trip departs at X i1 ¼ Hi ¼ 4 (step 2). Third line shows the earliest departure time max{0,30  3Hi} = 18 for the eighth trip, assuming that the last trip departs at X if i ¼ T  Hi ¼ 26 (step 3). Finally, fourth line shows the latest departure time T  2hi = 26 of the eighth trip, assuming that the last trip departs at X if i ¼ T (step 4). Therefore, the intersection of the earliest and latest departure times [max{14,18};min{30,26}] results in the feasible departure time window Di8 (see marked area of Fig. 4). By definition of window Dip we deduce an interesting properties. Property 1. Time window Dip determines all departure time values of X ip that might be feasible. If we consider the departure time window Dip of each trip p of line i, we can define the arrival time window Aipb for this trip at node b. We obtain Aipb by shifting Dip by t ib time units. It can be expressed as

h     i Aipb ¼ left Dip þ t ib ; right Dip þ tib ;

    where left Dip and right Dip represent lower and upper limit of time window Dip , respectively. Similarly, we define a synchronization window Sipb for each trip p of line i and node b as

h     i Sipb ¼ left Aipb þ wb ; right Aipb þ W b :

Fig. 3. Typical structure of the matrix formed by the synchronization variables Y ijpqb related to the trips of the lines i and j 2 J(i), and synchronization node b.

Fig. 4. Feasible departure time window Di8 for the eighth trip of line i corresponding to an instance with a planning period T = 30 min, fi = 10, and di = 0.33%.

O.J. Ibarra-Rojas, Y.A. Rios-Solis / Transportation Research Part B 46 (2012) 599–614

605

We can determine if synchronization between pth trip of line i and qth trip of line j is possible by looking at synchronization time window of trip p and arrival time window of trip q. We clarify this idea with the following theorem. T Theorem 1. For any trips p and q of lines i and j 2 J(i), respectively, and any synchronization node b 2 Bij of BTP, Sipk Ajqk ¼ ;, if ij and only if, Y pqk is forced to be zero due feasibility, and constraints (4) and (5) related to this indexes are redundant. T Proof. ()) Let p and q be two trips of lines i and j, respectively, such that Sipb Ajqb ¼ ; for a synchronization node b. Now, ij suppose that Y pqb ¼ 1 is feasible.     By Property 1, there exist feasible values X ip 2 Dip and X jq 2 Djq , such that Y ijpqb ¼ 1. Then, X ip þ tib 2 Aipb , X jq þ tjb 2 Ajqb and

h    i X ip þ tib þ wb ; X ip þ t ib þ W b # Sipb :     We have Y ijpqb is feasible, therefore wb 6 X jq þ t jb  X ip þ tib 6 W b . Then,

  h    i X jq þ tjb 2 X ip þ t ib þ wb ; X ip þ tib þ W b # Sipb :

  T Therefore, it follows that X jq þ tjb 2 Sipb Ajqb ¼ ;, which is a contradiction. () Let p and q be two trips of lines i and j, respectively, such that Y ijpqb must be zero for a synchronization node b. Now, T suppose Sipb Ajqb –;. Then, there exists an element a such that, a 2 Ajqb and a 2 Sipb .   By definition of window Ajqb , element a can be expressed as a ¼ X jq þ tjb , for some feasible X jq 2 Djq . By definition of   window Sipb , element a can be expressed as a ¼ X ip þ tib þ c, where wb 6 c 6 Wb, for some feasible X ip 2 Dip . Then,

            X jq þ tjb ¼ X ip þ t ib þ c ) X jq þ tjb  X ip þ tib ¼ c ) wb 6 X jq þ t jb  X ip þ t ib 6 W b :         Therefore, Y ijpqb ¼ 1 is feasible and wb  M 1  Y ijpqb 6 X jq þ t jb  X ip þ tib 6 W b þ M 1  Y ijpqb , i.e., satisfies constraints (4) and (5) which is a contradiction. h Using Theorem 1, we can eliminate useless variables and constraints. The preprocessing stage shown in Algorithm 1, calculates all the feasible departure, arrival, and synchronization time windows for each trip. Then, it implements Theorem 1 to identify when the synchronization of two trips is impossible, in this case, its corresponding synchronization variable and constraints of type (4) and (5) are eliminated. Algorithm 1. Preprocessing(BTP) Input: BTP instance Output: smaller formulation BTP0 1: for (i 2 I, j 2 J(i), p 2 {1, . . . , fi}, q 2 {1, . . . , fj}, and b 2 Bij) do 2: 3: 4:

Calculate Aipb and Sjqb   T if Aipb Sjqb ¼ ; then Eliminate variables Y ijpqb and constraints (4) and (5) for indexes i, j, p, q, and b.

5: end if 6: end for

After the preprocessing stage, the new BTP0 formulation has considerably fewer variables and constraints than the original BTP formulation. Therefore, the B&B performance can be improved. We present some preliminary results in Table 1 which shows the implementation of the B&B on BTP and BTP0 formulations using GAMS/CPLEX configured with default options, except for execution time limited to 3 hours. We used a Sun Ray terminal connected to a Sun Fire V440 server with 4 Ultra SPARC III processors running at 1602 GHz, and fitted with 8 GB RAM. Each test instance has a planning period of T = 200 min, three different lines, a frequency of fi = 20 for each line i, and one synchronization node for pair of lines (1, 2) and (1, 3), and the same waiting time windows for both nodes. Column one shows the minimum and maximum headway times. Column two shows the waiting time window for all synchronization nodes b. Columns three and five show the relative gap (difference between the best feasible solution and the best upper bound obtained by CPLEX) for BTP and BTP0 , respectively. Columns four and six show the computational time needed by the B&B to solve BTP and BTP0 , respectively. Finally, column seven shows the percentage of synchronization variables and constraints eliminated from BTP to obtain BTP0 . We do not specify the computational time of the preprocessing stage since it is negligible (less than one second for all experiments in this study).

606

O.J. Ibarra-Rojas, Y.A. Rios-Solis / Transportation Research Part B 46 (2012) 599–614

Table 1 Results of solving small instances of BTP and BTP0 using a B&B algorithm. [hi; Hi]

[9; 11] [8; 12] [5; 15] [9; 11]

[wb; Wb]

[3; 5] [3; 5] [3; 5] [1; 8]

BTP0

BTP

Reduction (%)

Gap (%)

Time (s)

Gap (%)

Time (s)

0 108.4 410.5 35

694 10,800 10,800 10,800

0 11 66 6.5

72 10,800 10,800 10,800

79 69 44 77

Note that even for small instances presented in Table 1, the formulations are sensitive to headway time variation (rows 2 and 3). Although, the variation in waiting time windows also generates difficult instances (row 4), the formulations are not as sensitive as they are in headway time variations. This behavior is also present in larger instances of the problem. However, note that in some instances that are difficult to solve (rows 2–4), there are high levels of eliminated variables using the preprocessing stage. The new BTP0 formulation has at least 44% fewer synchronization variables and constraints than the BTP formulation. This behavior is also present in larger instances. In particular, we see high percentage levels of eliminated variables for instances with small headway amplitude factor (di 6 30%). Therefore, the gap of the solutions obtained using BTP0 formulation are improved. However, there are still large gaps inclusive in small instances. We present an alternative solution methodology in Section 6 to obtain good quality solutions using few computational resources. 6. Multi-start iterated local search Obtain high quality solutions for timetabling problems in a short time is critical for transit networks. Usually, re-optimization is needed due to unacceptable solutions in subsequent problems of the entire planning process, or disturbances in the system such as failure of buses, and alterations of routes. Moreover, in our case of study, we must compute a timetable for each scenario based on different planning periods like rush hours, rainy day, holiday, and week-end. In this section, we propose a multi-start iterated local search to solve real-size instances of BTP using small computational resources. An iterated local search algorithm proceeds as follows (Lourenço et al., 2010). Given an initial solution s (which can be a local optimum) a perturbation is applied to it to obtain an intermediate feasible solution s0 . Then, a local search algorithm finds a new local optimum s0 ⁄. If this new solution meets an acceptance criterion, the incumbent solution is now s0 . We iterate until the algorithm reaches a stop criterion. The iterated local search is a simple and fast tool that can be improved by increasing the quality of each one of its modules, where interaction between intensification and diversification is critical, but challenging. Our proposed solution algorithm is exhibited by Algorithm 2. We define a number of main iterations (step 1). In each iteration, we do the following. First, we apply the preprocessing stage that we presented in Section 5. Then, we build an initial solution with Constructive(BTP0 ) (see Algorithm 3). Afterward, we try to improve the current solution by a local search algorithm at step 4 (see Algorithm 4). Next, we define the iterations of the iterated local search. In each one, we implement the perturbation phase (step 7) that searches for trips that can be forced to synchronize, followed by the local search to improve the perturbed solution (step 8), and we accept the new solution only if it represents an improvement in the total number of synchronizations (step 9). Algorithm 2. MILS(BTP) Input: BTP instance Output: a solution s for BTP 1: for (it = 1 to it = iterMULTI) do 2: BTP0 Preprocessing(BTP) 3: s Constructive(BTP0 ) 4: s Improve(s) 5: for (iter = 1 to iterILS) do 6: for (each pair of lines i and j that can synchronize) do 7: s0 Resynchronize(s,i,j) 8: s0 ⁄ Improve(s0 ) 9: s argmax{FBTP(s0 ⁄),FBTP(s)} 10: end for 11: end for 12: end for

In the following subsections, we present the details of the components Constructive(BTP0 ), Improve(s), and Resynchronize(s, i, j) of Algorithm 2. MILS is an algorithm based in the preprocessing idea that allows to determine the feasibility and search space along the algorithm steps. In Section 7, we show the performance of MILS.

O.J. Ibarra-Rojas, Y.A. Rios-Solis / Transportation Research Part B 46 (2012) 599–614

607

6.1. Construction of an initial solution Algorithm 3 shows our constructive phase of MILS. For each line in the bus network, we randomly generate the departure time of each one of its trips within its feasible departure time window Dip (step 3). Then, we update the departure time window of next trip (step 4). This manner we always obtain feasible solutions. Algorithm 3. Constructive(BTP0 ) Input: BTP0 formulation Output: solution s of BTP0 1: for (all lines i 2 I) do 2: for (p = 1 to p = fi) do   3: X ip Unif Dip h i i Dipþ1 \ X ip þ h ; X ip þ Hi 4: Dipþ1 5: end for 6: end for

6.2. Hill climbing local search Given an initial solution s, our hill climbing local search Improve(s) shown in Algorithm 4 identifies if it is possible to improve the objective function by forcing synchronization between pair of lines i and j at some node b using slight shifts. The steps of Improve(s) are the following. For all pair of lines i and j to be synchronized at node b (step 4), we identify each possible synchronization Y ijpqb (step 5). Then, we compute the new feasible departure time window Djq of q-th trip of line j (step 6). Then, we compute the departure time window D0jq that ensures the synchronization Y ijpqb (step 7). If D0jq is nonempty (step 8), we randomly generate the departure time of trip q of line j within D0jq (step 9). The new solution obtained by forcing the synchronization of lines i and j is accepted if it increases the objective function value (steps 12). Note that the implementation of this local search could lead to different solutions since there is a random factor in the departure time generation for all trips. Algorithm 4. Improve(s) Input: Feasible solution s Output: local optimum s⁄ 1: s⁄ s 2: while (flag_improve = 1) do 3: flag_improve = 0 4: for (each pair of lines (i,j) to be synchronized at each node b 2 Bij) do 5: 6: 7: 8: 9:

for (each feasible synchronization Y ijpqb ¼ 0) do h i h i j j Djq X jq1 þ h ; X jq1 þ Hj \ X jqþ1  Hj ; X jqþ1  h    i Th i  i Djq D0jq X p þ tb  tjb þ wb ; X ip þ t ib  t jb þ W b 0

if (D j – ;) then   X jq Unif D0jq ; Y ijpqb

1

10: end if 11: end for 12: Update incumbent solution s⁄ and flag_improve 13: end for 14: end while

6.3. Perturbation and acceptance criterion The perturbation is represented in Algorithm 5. The main idea is to force the synchronization of a pair of lines (i0 , j0 ) at 0

0 0

some synchronization node b 2 Bi j . This algorithm synchronizes a specific pair of lines, but in general, the objective value could decrease as is the case with any perturbation movement. First, there is an initialization phase that clears the departure

608

O.J. Ibarra-Rojas, Y.A. Rios-Solis / Transportation Research Part B 46 (2012) 599–614

times for all trips of lines i0 and j0 (step 1). Then, for each trip p of line i0 (step 2), we compute its new feasible synchronization h 0 0 0 i0 time window Sipb0 considering the last assigned trip (steps 3), i.e., Sipb0 ¼ X ilast p þ ðp  last pÞh þ 0

0

tib0 þ wb0 ; X ilast

0

p

0

þ ðp  last pÞHi þ t ib0 þ W b0 . Next, for trip q of line j0 possible to synchronize with trip p of line i0 (step 5), 0

we compute its new feasible arrival time window Ajqb0 considering the last assigned trip (steps 6), i.e., h 0 i 0 0 0 0 0 j0 Ajqb0 ¼ X jlast q þ ðq  last qÞh þ t jb0 ; X jlast q þ ðq  last qÞHj þ tjb0 . If the arrival time window of trip q intersects with the synchronization time window of trip p (step 7), we randomly generate the departure time of trips p and q of lines i0 and j0 , respectively, to ensure their synchronization (steps 8 and 10). Finally, we calculate the departure times of the unassigned trips considering its nearest (lower and upper) assigned trips (steps 15–19). Algorithm 5. Resynchronize(s, i0 , j0 ) 0 0

0

Input: Feasible solution s, lines i0 , j0 2 J(i0 ), and node b 2 Bi j Output: New solution s0 1: Clear departure times for all trips of lines i0 and j0 , last_p

0, last_q

1

0

2: for ðp ¼ 1 to p ¼ f i Þ do 0

0

3:

Compute the new feasible synchronization time window Sipb0 considering X ilast

4:

q

5:

while ðq 6 f j Þ do

6: 7: 8: 9: 10:

p

last_q 0

0

0

Compute the new feasible arrival time window Ajqb0 considering X jlast q  00  T j0 j i0 if Y ipqb Aqb0 –; then 0 0 ¼ 1 is feasible; i:e:; S pb  0 T 0  0 0 0 T 0 0 Unif Sipb0 Ajqb0  t jqb0 Generate X jq to arrive within Sipb0 Ajqb0 at node b, i.e., X jq h 0 i 0 0 0 0 j i0 Ensure synchronization Y ipqb Unif X jq  W b0 ; X jq  wb0  tib0 0 , i.e., X p 0 0

j Y ipqb 0

1, last_p = p, and last_q = q

11: end if 12: q q+1 13: end while 14: end for 15: for (each unassigned trip p of line l 2 {i0 ,j0 }) do 16: Search nearest assigned trips p1 < p < p2 17: 18:

Calculate new feasible departure time window Dlp considering X lp1 and X lp2   X lp Unif Dlp

19: end for

7. Experimental results The Bus Synchronization Timetabling problem is NP-hard as we discussed in Section 4. Nevertheless, it is crucial to solve real-size instances for bus companies. Therefore, in this section we show the limits of a B&B algorithm and experimentally analyze the performance of MILS metaheuristic. The algorithms were implemented as follows. B&B is the branch and bound algorithm to solve mixed integer linear programs implemented by GAMS/CPLEX 11.2 using default options, except for optimal the criterion fixed at 0.0001 percent and the execution time limited to two hours. It was executed on a Sun Ray terminal connected to a Sun Fire V440 server, equipped with 4 Ultra SPARC III processors running at 1602 GHz, and fitted with 8 GB of RAM. MILS After some preliminary tests, the iteration limits of MILS were set to iterILS = 200 and iterMULTI = 10, and we take the best solution with its computational time. MILS was coded in C++ and executed on a iMac OS X with an Intel Core 2 Duo 3.06 GHz processor, and 4 GB of RAM. To perform the experimental analysis, we design an instance generator (coded in C++) based on information of the parameters provided by an expert from a Monterrey’s bus network. We build instances of different sizes. The size is represented for the number of lines jIj and the number of different pairs of lines to synchronize along the bus network. As we have already pointed out in Table 1, the headway amplitude factor di is directly responsible for the difficulty of the BTP instances as well as

609

O.J. Ibarra-Rojas, Y.A. Rios-Solis / Transportation Research Part B 46 (2012) 599–614

the number of bus lines, and synchronization nodes. Therefore, the larger the headway amplitude factor is, the harder is the instance. We use two intervals to generate the headway amplitude factor, details for the eight instance types are shown in Table 2. All instance types have the following common characteristics. A planning period of T = 240 min to represent, for example, the period corresponding to an afternoon. The frequency fi for each line i is between [7; 11]. The travel time t ib from the depot to each synchronization node b for each line i is within [20; 60]. The minimum and maximum waiting times for each synchronization node b are within [3; 5] and [9; 12], respectively. Finally, the number of different pairs of lines to synchronize at each node b is within [1; 7]. We randomly generate ten instances for each one of the eight instance types to analyze the algorithms performance. We suspect that the upper bounds based on linear relaxations may be of poor quality. To palliate this behavior, we can deduce a simple but effective upper bound as follows:

XX X

n o min Pijb ; Q ijb ;

i2I j2JðiÞb2Bij

where, for each pair of lines (i, j) to synchronize at some node b; Pijb represents the quantity of trips of line i that are possible to synchronize with some trip of line j at node b. Similarly, Q ijb represents the quantity of trips of line j that are possible to synchronize with some trip of line i at node b. Notice, that for all instance types we use small waiting time windows compared with the minimum headway times of lines. Since departure times of consecutive trips must be separated at least by the minimum headway time which is larger than the length of waiting time windows, it is possible to synchronize a trip p of some line i with at most one trip q of each line j 2 J(i). To compare MILS with B&B, we consider ten instances of each type and report their gaps (calculated as the relative difference between the best feasible solution and the best upper bound between our proposed upper bound and the one obtained by CPLEX). Table 3 shows the relative mean gap and the mean computational time in seconds for the implementation of algorithms B&B (column ‘‘B&B’’) and MILS (column ‘‘MILS’’), respectively. Although the B&B performance is significantly improved by the preprocessing stage, we can see that except for the smallest instance type A1, the B&B did not solve optimally any of the instances in two hours. The gaps of MILS are small for large instances (A6–A8) and are obtained in less than one minute. Notice that B&B algorithm gives slightly better gaps for small instances (A1–A5), but the computational times of MILS are less than 15 s which is very important for the companies in our case of study. Our proposed algorithm MILS is the first attempt to obtain solutions for the synchronization timetabling problem that fits the transit network of Monterrey City. This algorithm is beneficial to the planner of real bus networks because real instances comprise around 200 lines, and it is often necessary to re-calculate timetables to consider subsequent planning subproblems like vehicle scheduling. We can also rapidly generate a timetable for each scenario considered by different companies. The bus network expert of Monterrey who participated in this study confirmed that the solutions generated by MILS are of good quality (based in the number of synchronizations), and more important, these solutions are obtained fast enough for their requirements.

8. Conclusions In this study, we addressed Monterrey’s bus network which is similar to transit networks present in Latin American cities. The bus network is composed of several private companies and the government does not have a strong influence in the operation of the routes. We focus on timetables generation meant for companies, while passengers only have knowledge of an estimated waiting time with respect to the bus frequency. We solve two main problems at some specific nodes of the network: avoid bus bunching of different lines and allow passenger transfer. Bus network planners in Monterrey do not have access to real-time optimization tools. Therefore, one sought to optimize the performance of the whole bus network by controlling the departure time of the buses in the planning stage. We propose a new formulation for the Synchronization Bus Timetabling Problem (BTP) that is sufficiently flexible to model Monterrey bus network. We prove that BTP is NP-hard, and this complexity proof also ensures the NP-hardness of problems presented in Ceder et al. (2001) and Eranki (2004). We find a special structure of the solution space of BTP, which gives us the possibility of eliminating many decision variables and constraints using a preprocessing stage. Although this preprocessing stage improves the B&B performance, the intractability of small instances remains.

Table 2 Size of the eight instance types for the BTP. Instance

A1

A2

A3

A4

A5

A6

A7

A8

jIj Nodes % di

15 3 [5; 15]

15 3 [10; 25]

40 8 [5; 15]

40 8 [10; 25]

100 20 [5; 15]

100 20 [10; 25]

200 40 [5; 15]

200 40 [10; 25]

610

O.J. Ibarra-Rojas, Y.A. Rios-Solis / Transportation Research Part B 46 (2012) 599–614 Table 3 Results of the instances solved by CPLEX’s B&B and MILS. Instance

A1 A2 A3 A4 A5 A6 A7 A8

B&B

MILS

Gap (%)

Time

Gap (%)

Time

0 8.54 12.41 11.84 14.76 22.49 22.04 22.64

546.82 7200 7200 7200 7200 7200 7200 7200

1.51 12.35 14.23 11.85 16.37 14.48 17.57 15.55

0.42 0.67 2.46 1.58 12.87 11.41 45.09 52.09

Obtaining a fast solution algorithm is beneficial to the planner in real-life transit networks since re-optimization is needed, for example, to improve the quality of a timetable considering other subproblems of the entire planning process. In this study, we use our preprocessing stage to define a multi-start Iterated Local Search algorithm (MILS) that obtains solutions with less than 18% of relative gap in less than one minute of computational time. These gaps have to be carefully analyzed since BTP is a very complex problem and it is quite difficult to obtain an accurate measure of the solution quality. We recall that formulation BTP was intractable even with the implementation of the preprocessing stage since small instances (three lines, and two synchronization nodes) take hours to be solved by a CPLEX’s B&B. Now, we are able to solve instances of 200 lines and 40 synchronization nodes using MILS. Other aspects that should be considered for further studies are the following:  We present a simple calculation to obtain an upper bound. In half of the experiments, our proposed upper bound is tighter than the bound obtained by CPLEX in two hours of computational time. Thus, it is important to implement other methodologies to obtain tighter bounds.  As we mentioned in Section 6, enhancements of MILS can be done by improving the quality of its modules. It could be useful to explore alternative strategies to calculate departure times replacing the random factor present in our local search algorithm.  In real cases, it is necessary to integrate passenger demand behavior to determine more accurate planning periods. Also, forecasting travel times might bring important improvements to the formulation representativeness.  Obviously, it would be necessary to consider subsequent steps, such as vehicle scheduling, in the timetabling problem. However, this latter point of view for the synchronization timetabling lacks treatment in the literature, since resource considerations are rarely taken into account (Guihaire and Hao, 2008a).

Acknowledgments This study was partially supported by Grants C003-2009-01/113032 and 101857 from the Mexican National Council for Science and Technology (CONACYT). We thank Das-Sistemas, Fernando López and Ada Álvarez for the fruitful discussions. Omar J. Ibarra-Rojas wishes to acknowledge graduate scholarship from CONACYT. Finally, we express our gratitude to M.A. Salazar-Aguilar and to the anonymous referees for their valuable comments that helped to improve the quality of this paper. Appendix A. Computational complexity proof Studies of Ceder et al. (2001) and Eranki (2004) remark on the intractability of similar problems as we noted with BTP on empirical results. In this appendix, we prove that BTP belongs to the NP-hard class. With this, we also prove that formulations of Ceder et al. (2001) and Eranki (2004) are NP-hard. To prove that an optimization problem X is NP-hard, we have to prove that its decision version is NP-complete. As it can be seen in Garey and Johnson (1979) and Papadimitriou and Steiglitz (1998), there are basic steps to prove that a decision problem X is NP-complete. These steps are the following. 1. Prove that problem X 2 NP. 2. Exhibit that a NP-complete problem Y can be polynomially reduced to X (denoted as Y  pX). 3. Prove the solution equivalence, i.e., the answer of an instance of problem X is ‘‘YES’’ if and only if, the answer of the related instance of problem Y is ‘‘YES’’. We propose a polynomial reduction from monotone NAE-3SAT, whose NP-completeness is assured by Schaefer (1978), to the decision version of BTP. First of all, we present the decision versions of both problems.

O.J. Ibarra-Rojas, Y.A. Rios-Solis / Transportation Research Part B 46 (2012) 599–614

611

Monotone NAE-3SAT: INPUT: A set of r literals and a set {C1, C2, . . . , Cc} of c clauses. Each clause has three unnegated variables. A clause is true if the values of its literals are not all equal. QUESTION: Is there an assignment of the literals values (false as 0 and true as 1) such that the sentence / = C1 ^ . . .^ Cc is true?

An example of this problem is the sentence: / = (x1 _ x2 _ x3) ^ (x4 _ x2 _ x3) ^ (x1 _ x4 _ x3) that has a true value when x1 = x3 = 1 and x2 = x4 = 0. Bus synchronization timetabling problem (decision version called dBTP): INPUT: Parameters I; Bij ; T; f i ; t ib , [hi;Hi], [wb; Wb] " i, j 2 I(i), b 2 Bij, and a scalar K > 0. QUESTION: Are there any values for variables X ip and Y ijpqb of BTP "i, j 2 I, p = 1, . . . , fi, j 2 J(i), q = 1, . . . , fj, b 2 Bij that satisfy constraints (1)–(7) and such that fi X fj XXXX

Y ijpqb P K?

i2I j2JðiÞ b2Bij p¼1 q¼1

Next, we show the theorem of the complexity for dBTP. Theorem 2. dBTP is NP-complete. Proof. The description of the steps for the complexity proof are the following. 1. dBTP 2 NP Assume there exist an algorithm that generates a solution for dBTP. Determining if the solution is feasible for dBTP implies several steps. Constraints (1)–(3) represent a quantity of inequalities in the order of the number of variables. Constraints (4) and (5) verify two inequalities for each synchronization variable Y ijpqb . Finally, the last step is to verify if the sum of the synchronization variables is greater than or equal to K (one constraint). Considering all the constraints, we need a polynomial number of steps to verify feasibility of a dBTP solution, i.e., dBTP 2 NP. 2. Monotone NAE-3SAT p dBTP Consider an arbitrary instance of monotone NAE-3SAT with c clauses. We build a particular instance of dBTP, called dBTP⁄, as follows. Set the number of lines with jIj = r (one line corresponds to each literal of monotone NAE-3SAT). The planning period is T = 1 min. The frequency of each line is fi = 1. The headway-time windows are [hi;Hi] = [0; 1]. The waiting times are set to wb = Wb = 0, "b, i.e, the synchronization of two trips exists when they arrive simultaneously at some node. To define the synchronization nodes, we make Bfg = ; for all pair of lines. Then, for each clause Cl = (xf _ xg _ xh), we create six synchronization nodes labeled as {l1, l2, . . . , l6}. If clause Cl = (xf _ xg _ xh) with three different literals, we add two synchronization nodes for each pair of lines as follows: S S S Bfg = Bfg {l1, l2}, Bgh = Bgh {l3, l4}, and Bfh = Bfh {l5, l6}. We define a disjoint synchronization time window of 9-min for these clauses as follows. The beginning of the 9-min period of clause Cl is ol = 9(l  1). Next, we define the arrival times for each node shifted by 1 minute to define a synchronization if two trips are separated by 1 min. Explicitly, for the pair of lines (f, g): tfl1 ¼ ol þ 1; tgl1 ¼ ol þ 0, t fl2 ¼ ol þ 2, and t gl2 ¼ ol þ 3. For the pair of lines (g,h): tgl3 ¼ ol þ 4, thl3 ¼ ol þ 3; tgl4 ¼ ol þ 5, and t hl4 ¼ ol þ 6. Finally, for the pair of lines (f,h): t fl5 ¼ ol þ 6; thl5 ¼ ol þ 7, t fl6 ¼ ol þ 9, and thl6 ¼ ol þ 8. Fig. A.5 shows the 9-min synchronization time window of the first clause (l = 1 and ol = 0). The quantities above nodes represent arrival times when a trip departs at time 0. If clause Cl has two of its three literals equal, we consider only two different lines. Explicitly, if Cl = (xf _ xf _ xg), we add S only synchronization nodes for the pair of lines (f, g) as Bfg = Bfg {l1, l2} using the same route design. Fig. A.6 shows the 9-min synchronization time interval for first clause (l = 1 and ol = 0). Finally, we define K = 2m + n for instance dBTP⁄, where m and n are the number of clauses of monotone NAE-3SAT instance with three and two different literals, respectively (c = m + n). The case of a clause with three equal literals is not considered, since it is a trivial case and there is no solution that provides a ‘‘YES’’ answer to the monotone NAE-3SAT. Note that this

612

O.J. Ibarra-Rojas, Y.A. Rios-Solis / Transportation Research Part B 46 (2012) 599–614

Fig. A.5. Line design for the 9-min interval related to the first clause C1 with three different literals.

Fig. A.6. Route design for the 9-min interval related to the first clause C1 with two different literals.

reduction implies a number of steps bounded by the number of literals and the number of clauses of monotone NAE-3SAT. Also, the size of instance dBTP⁄ depends on this numbers. 2.1 Properties of dBTP⁄ We assign the value of xi 2 {0,1} in monotone NAE-3SAT to the variable X ip in dBTP⁄. Therefore, if xi – xj, lines i and j synchronize at each common node. Fig. A.7 shows the cases of a true clause C1 = (xf _ xg _ xh) where the literal xf has a different value than xg and xh. A dashed line represents the synchronization of two lines. Notice that the synchronized pairs of lines are (f, g) and (f, h) since xf – xg and xf – xh. By dBTP⁄ definition, instance dBTP⁄ has the following properties which are useful to prove the solution equivalence. h Property 2. Any clause Cl of monotone NAE-3SAT with three different literals implies two synchronizations in instance dBTP⁄ if it is true, and 0 synchronizations if it is false. Property 3. Any clause Cl of monotone NAE-3SAT with two of its three different literals equal, implies one synchronization in instance dBTP⁄ if it is true, and 0 synchronizations if it is false. Fig. A.8 shows the cases of a true clause C1 = (xf _ xg _ xh). For the dBTP⁄ instance, the nodes represent the lines and the edges represent a synchronization between the connected lines.

Fig. A.7. Possible synchronizations corresponding to the pairs of lines (f, g) and (f, h). Case (a) represents lines f, g, and h with departure times of 0, 1, and 1, respectively. Case (b) represents lines f, g, and h with departure times of 1, 0, and 0, respectively.

Fig. A.8. Possible values for a true clause of monotone NAE-3SAT and the related pairs of synchronized lines.

O.J. Ibarra-Rojas, Y.A. Rios-Solis / Transportation Research Part B 46 (2012) 599–614

613

3. Solution equivalence ()) Let X = (x1, . . . , xc) be a solution of the monotone NAE-3SAT problem with a ‘‘YES’’ answer. Then, each clause Cl is true and we make X i1 ¼ xi ; 8i 2 I in instance dBTP⁄. We also assign the value of 1 to variables Y ij11b ; 8 i–j 2 JðiÞ; b 2 Bij , such that xi, xj 2 Cl, and xi – xj (constraints of dBTP⁄ are trivially satisfied). By Property 2 and 3, there are 2m + n synchronizations in instance dBTP⁄ where c = m + n is the number of clauses; m clauses with three different literals, and n clauses with two different literals, i.e.,

XXX

Y ij11b ¼ ð2m þ nÞ P K:

i2I j2JðiÞ b2Bij

() Let X and Y be a solution of instance dBTP⁄, with a ‘‘YES’’ answer. By Properties 2 and 3,

XXX

Y ij11b ¼ ð2m þ nÞ ¼ K:

i2I j2JðiÞ k2Bij

Now, suppose there is a false clause Cl. By Properties 2 and 3, the number of synchronizations related to this false clause is 0. P P P This implies i2I j2JðiÞ b2Bij Y ij11b < ð2m þ nÞ ¼ K, which is a contradiction. Considering an optimization problem, where its decision version belongs to the NP-complete class, then, the optimization problem belongs to the NP-hard class. Thus, the following corollary is deduced. Corollary 1. BTP is NP-hard. The concept of the previous polynomial reduction from monotone NAE-3SAT to dBTP can be used for the BTP formulation with real variables X ip and also for the problems presented in Ceder et al. (2001) and Eranki (2004). On the one hand, the particular instance dBTP⁄ is precisely an instance of Ceder’s synchronization timetabling problem. On the other hand, the instance construction of Eranki’s problem can be done in a manner similar to dBTP⁄, but considering a constant waiting time window [w;W] for each node b, and shifting arrival times at synchronization nodes by w + 1 minutes instead of 1 minutes. Thus, the polynomial reduction from NAE-3SAT to Eranki’s instance is straight forward. References Adamski, A., 1993. Public transport lines synchronization problem. In: Automatics, Scientific Bulletins of AGH, vol. 64, pp. 221–234. Bookbinder, J.H., Desilets, A., 1992. Transfer optimization in a transit network. Transportation Science 26 (2), 106–118. Bruno, G., Improta, G., Sgalambro, A., 2009. Models for the schedule optimization problem at a public transit terminal. OR Spectrum 31 (3), 465–481. Ceder, A., 2007. Public Transit Planning and Operation: Theory, Modeling and Practice. Elsevier, Butterworth-Heinemann. Ceder, A., Golany, B., Tal, O., 2001. Creating bus timetables with maximal synchronization. Transportation Research Part A 35 (10), 913–928. Ceder, A., Tal, O., 2001. Designing synchronization into bus timetables. Transportation Research Record: Journal of the Transportation Research Board 1760, 28–33. Cevallos, F., Zhao, F., 2006. Minimizing transfer times in public transit network with genetic algorithm. Transportation Research Record 1971, 74–79. Chakroborty, P., Deb, K., Porwal, H., 1997. A genetic algorithm based procedure for optimal transit system scheduling. In: Proceedings of 5th International Conference on Computers in Urban Planning and Urban Management, Mumbai, India. pp. 330–341. Chakroborty, P., Deb, K., Subrahmanyam, P.S., 1995. Optimal scheduling of urban transit systems using genetic algorithms. Journal of Transportation Engineering 121 (6), 544–553. Chowdhury, S., Chien, S., 2001. Optimization of transfer coordination for intermodal transit networks. In: Proceedings of the 80th Annual Meeting, Transportation Research Board. Paper No. 01-0205, Washington, DC. Chung, E., 1990. Transfer coordination model and real-time strategy for inter-modal transit services. Ph.D. thesis. Graduate Department of Civil Engineering, University of Toronto. Clever, R., 1997. Integrated timed transfer: a European perspective. Transportation Research Record 1571, 107–115. Daduna, J., Voß, S., 1995. Practical experiences in schedule synchronization. In: Daduna, J., Branco, I., Paixao, J. (Eds.), Computer-Aided Scheduling of Public Transport, Lecture Notes in Economics and Mathematical Systems, vol. 430. Springer, Berlin, pp. 39–55. Daganzo, C.F., 2009. A headway-based approach to eliminate bus bunching: Systematic analysis and comparisons. Transportation Research Part B 43 (10), 913–921. Dessouky, M., Hall, R., Nowroozi, A., Mourikas, K., 1999. Bus dispatching at timed transfer transit stations using bus tracking technology. Transportation Research Part C 7 (4), 187–208. Eranki, A., 2004. A model to create bus timetables to attain maximum synchronization considering waiting times at transfer stops. Master’s thesis. Department of Industrial and Management Systems Engineering, University of South Florida. Fleurent, C., Lessard, R., 2009. Integrated Timetabling and Vehicle Scheduling in Practice. Technical Report. GIRO Inc. Montreal, Canada. Fleurent, C., Lessard, R., Séguin, L., 2005. Transit timetable synchronization: evaluation and optimization. In: Proceedings of the 9th International Conference on Computer-Aided Scheduling in Public Transport (CASPT), San Diego, California. pp. 9–11. Garey, M.R., Johnson, D.S., 1979. Computers and Intractability: A Guide to the Theory of NP-Completeness, first ed. W.H. Freeman. Guihaire, V., Hao, J.K., 2008a. Transit network design and scheduling: a global review. Transportation Research Part A 42 (10), 1251–1273. Guihaire, V., Hao, J.K., 2008b. Transit network re-timetabling and vehicle scheduling. In: Proceedings of the 2nd International Conference on Modelling, Computation and Optimization in Information Systems and Management Sciences (MCO’08). Springer, Heidelberg, pp. 135–144. Hall, R., Dessouky, M., Lu, Q., 2001. Optimal holding times at transfer stations. Computers &Industrial Engineering 40 (4), 379–397. Van den Heuvel, A., Van den Akker, J., Van Kooten, M., 2008. Integrating timetabling and vehicle scheduling in public bus transportation. Technical Report UU-CS-2008-003. Department of Information and Computing Sciences, Utrecht University, Utrecht, The Netherlands. Jansen, L.N., Nielsen, O.A., 2002. Minimizing passenger transfer times in public transport timetables. In: Proceedings of the 7th Conference of the Hong Kong Society for Transportation Studies, Transportation in the Information Age, Hong Kong. pp. 229–239. Li, Z.C., Lam, W., Wong, S., Sumalee, A., 2010. An activity-based approach for scheduling multimodal transit services. Transportation 37 (5), 751–774. Liu, Z., Shen, J., Wang, H., Yang, W., 2007. Regional bus timetabling model with synchronization. Journal of Transportation Systems Engineering and Information Technology 7 (2), 109–112. Lourenço, H., Martin, O., Stützle, T., 2010. Iterated Local Search: Framework and Applications. Handbook of Metaheuristics, vol. 146, second ed. pp. 363–397.

614

O.J. Ibarra-Rojas, Y.A. Rios-Solis / Transportation Research Part B 46 (2012) 599–614

Papadimitriou, C.H., Steiglitz, K., 1998. Combinatorial Optimization: Algorithms and Complexity. Dover Publications. Rabut, P., 2010. Système de transport urbain de l’aire métropolitaine de Monterrey (Mexique) Analyse de l’offre et de la demande. Master’s thesis. Université de Franche-Comté. France. Schaefer, T., 1978. The complexity of satisfiability problems. In: Proceedings of the 10th Annual ACM Symposium on Theory of Computing, San Diego, California, US. pp. 216–226. Schröder, M., Solchenbach, I., 2006. Optimization of Transfer Quality in Regional Public Transit. Technical Report 84. Berichte des Fraunhofer-Instituts für Techno- und Wirtschaftsmathematik. Wong, R.C.W., Yuen, T.W.Y., Fung, K.W., Leung, J.M.Y., 2008. Optimizing timetable synchronization for rail mass transit. Transportation Science 42 (1), 57– 69.