Towards transfer synchronization of regularity-based bus operations with sequential hill-climbing

In this work we model and discuss how we can achieve coordination between different bus service lines. Key problem challenges are (a) the multiple conflicting priorities (on one hand the improvement of bus service regularity and on the other hand the reduction of passenger transfer waiting times) and (b) the computational complexity for re-scheduling the dispatching times of bus trips for meeting the conflicting priorities. Initially, a model for reducing the waiting times at bus transfer stations while also improving the operations of regularity-based bus services subject to operational constraints is introduced. Conflicting priorities are handled with the introduction of weight factors that allow bus operators to decide the trade-off between improvement of regularity-based operations and reduction of passenger waiting times at transfer stations. After that, an exterior point penalty function is introduced for handling operational constraints and a sequential hill-climbing search strategy is applied for converging to an approximate optimal solution. For our case study, we utilize general transit feed specification data from two regularity-based bus services in central Stockholm that intersect in five transfer stations. Experimental tests showcase a 13% potential waiting time improvement at transfer stations while sacrificing only 2.8% of service regularity and satisfying all operational constraints.


Introduction
Attempts to improve the service quality of high and low frequency bus services have led to the introduction of several key performance indicators (KPIs) for bus operations such as the On Time Adherence (OTA) to the planned schedule or the excess waiting time (EWT) of passengers at bus stations.OTA is a widely used KPI for punctuality-based bus services where the performance of the bus operations is measured by comparing the real arrival times of buses at bus stops against the planned ones.On the other hand, high-frequency services in metropolitan areas have difficulties adhering to the scheduled arrival times at stops and for this reason the bus operations are run on a regularity basis.Regularity-based bus operations, which are the main focus of this work, do not try to adhere to the scheduled arrival times at the bus stops but aim to achieve a certain level of regularity by avoiding bus bunching and significant headway variations.For instance, each bus trip does not have a specific planned arrival time at the bus stops and the only objective is the minimization of the headway variations among buses (keeping the actual headways as close to the scheduled ones as possible).For this reason, regularity-based services utilize different KPIs such as the EWT which represents the variation of the actual waiting times of passengers at bus stops from the planned ones.
Offline and real-time control measures such as departure time changes, speed control and bus holding have been proposed from several studies that try to meet the Public Transport Authorities (PTA) KPIs related to passenger waiting times at stops and in-vehicle travel times (Dessouky et al. 2003;Cats et al. 2012;Gkiotsalitis and Maslekar 2015;Fu and Yang 2002;Xuan et al. 2011).Other works, such as Yu et al. (2011), Chowdhury and Chien (2011), were more general by not focusing on PTA KPIs but trying to find an acceptable balance between passenger and operational costs to maximize service quality while reducing operational costs.Apart from meeting the PTA KPIs for each individual bus service, bus operators should also reduce the total travel time of passengers to ensure that their services remain competitive.This was also discussed in Botea et al. (2013), Botea and Braghin (2015) that focused on the significant travel time increase of passengers who change several modes due to the operational uncertainty.
Reducing the overall passenger travel time for passengers that use multiple lines for completing their journeys requires the minimization of the waiting times at bus stations that serve as interchange points among different bus services.Meeting the bus service KPIs while ensuring acceptable waiting times at bus interchange points are often contradictory problems because attempts to improve the coordination of bus services might deteriorate the individual performance metrics of each bus service.For instance, a bus service A can improve its EWT by re-scheduling its planned operations until the EWT of passengers is minimized.The same holds true for another bus service B which intersects with service A at a common bus stop T (transfer station).However, the optimal schedules that minimize the EWTs of services A and B might lead to significant waiting times at the transfer stop T for the passengers who are willing to transfer from bus service A to service B and vice versa.If one wants to reduce the transfer 1 3 Transfer synchronization of regularity-based bus operations waiting times at the stop T for improving the coordination between the two services, the services A and B should be re-scheduled and, consequently, they will no longer be optimal regarding the passengers' EWT.Therefore, it is evident that the coordination of various bus services is a trade-off between the reduction of waiting times at the transfer stops and the minimization of the passengers' EWT of each specific service.

Related studies
In the work of Ceder et al. (2001) bus timetables were re-scheduled for achieving maximal synchronization and the problem was formulated as a mixed integer programming one.However, the synchronization was achieved by simultaneous arrivals at transfer stations without considering the implications to the performance of each bus line.Daduna and Voß (1995) and Domschke (1989) worked also on scheduling timetables with a specific emphasis on synchronizing the transfers among different lines.The perplexities of such task included the headway differences among different lines that complicate the harmonization of passenger transfers.For this reason, Daduna and Voß (1995) adopted a multi-objective approach where the minimization of the sum of all transfer waiting times and the operating costs were considered by a heuristic optimization approach while an upper threshold for the maximum allowed waiting time was adopted.
In their work, Bookbinder and Desilets (1992) tried to reduce the transfer times of passengers by modifying the scheduled departure times of bus trips.For selecting the optimal departure time changes, a simulation-based approach was employed for considering the stochastic nature of the actual travel times of bus trips.Voss (1992) focused also on minimizing the transfer times of passengers in transit systems (buses and trains).The decision variables of the transfer time minimization problem were the departure time changes of scheduled trips (that affect the arrival times of trips at the transfer stops) and the problem was considered at its simpler version where all transit lines are met only at certain points and its more complicated version where they partially use the same tracks.Jansen et al. (2002), Wong et al. (2008), Vansteenwegen and Van Oudheusden (2006) worked also on the problem of reducing the waiting time of passengers at the transfer stops while keeping the daily trips evenly spaced.Similarly to the approach of Bookbinder and Desilets (1992), the more recent studies of Cevallos and Zhao (2006a) and Cevallos and Zhao (2006b) tried also to reduce the complexity of the problem by merely shifting the departure times of the scheduled trips for finding the optimal solution for the passenger transfers with the use of a genetic algorithm.
Apart from the synchronization of the bus timetables, Zhao et al. (2003) focused on the more computationally complex problem of real-time coordination via a decentralized strategy where stops and buses communicate in real-time to achieve dynamic bus holding coordination at a local level.
There is a series of works also on bus and rail operations coordination.Coffey et al. (2012) optimized the time schedules of public transport modes for reducing missed connections under a multi-modal setting by trying to match the passenger demand expressed via journey planners with the public transport supply.Chien and Schonfeld (1998), Sun and Hickman (2004), Shrivastava and O'Mahony (2009), Sivakumaran et al. (2012), Verma and Dhingra (2006) proposed multi-modal coordination approaches by considering the 'feeder model' where the bus system feeds a rail transit line and bus services have to adjust their schedules to the rail schedules.Nevertheless, the 'feeder line' concept is not applicable if the system is not multi-modal, but contains only bus transfers which are expected to have almost similar importance/weight.This is also the direction of our work which models the bus coordination problem considering the passenger waiting times at transfer stations and the adherence to the KPIs of each bus service.The motivation behind our work is the non-evenly distributed passenger waiting times at transfer stations that can penalize the total passenger travel times even if the regularity of single bus services is optimal.In Sect. 3 the bus coordination problem is formulated considering the waiting times of passengers at transfer stations, the regularity-based KPIs of bus services (passenger EWT) and the operational constraints.The modeled bus coordination problem has exponential computational complexity and is computationally intractable even for small scenarios with two bus lines.For this reason, in Sect. 4 we present a heuristic solution space search strategy for converging to a nearly optimal solution that improves coordination by rescheduling the departure times of the scheduled bus trips.Experimental results from two high-frequency, bi-directional services in Stockholm are presented in Sect. 5 showcasing a ∼ 13% transfer waiting time improvement for only 2.8% deterioration of ser- vice regularity.

Modeling the bus coordination problem
Let us assume that are the bus services on a city net- work, where |Q| is the total number of available bus services.Each bus service is the total number of daily trips for bus service R Q j .R Q j is already defined from the tactical planning phase of the broader Transit Network Planning problem (TNP) including the network design, timetable development and crew/vehicle scheduling for each bus line (Kepaptsoglou and Karlaftis 2009;Farahani et al. 2013;Ceder 2007;Gkiotsalitis and Cats 2017).The tactical planning phase defines also a number of operational constraints such as the expected service frequency during different daily periods (i.e., morning peak, afternoon peak) for adjusting to the passenger demand variations during the day.Let us assume that for each bus service Q j ∈ Q the solution of the TNP has defined a frequency range set:

3
Transfer synchronization of regularity-based bus operations where F 1,Q j = (a 1 , b 1 ) is the minimum and maximum allowed departure time differ- ence between consecutive bus trips of the same service during the first time period of the day for bus service F Q j and | F Q j | is the number of time period splits within a day that enjoy different frequency settings (i.e., F 1,Q j : morning peak, F 2,Q j : midday etc.).
Each bus service Q j serves also a number of bus stations is the total number of bus stations served by service Q j .In dynamic control operations, the bus service operator can re-schedule the bus service departure times (timetable re-scheduling).Given the bus travel times from one bus station to another during different day periods, one can calculate the arrival time of a bus trip at every bus station: is the expected travel time of bus trip R Q j between two consecutive bus stations i to i + 1 plus the associated dwell time, t R Q j ,S |S Q j |,Q j the arrival time of bus trip R Q j at the terminal and t R Q j ,S 1,Q j the departure time of bus trip R Q j from the origin that can be re-scheduled.Let x R Q j be the decision variable of the problem which represents the time variation of departure time t R Q j ,S 1,Q j after the re-scheduling and receives integer values (those integer values represent the time changes in minutes), then the arrival time of the bus trips at every station after re-scheduling is: For coordinating two bus services Q i , Q j , first there should be a search for the com- mon interchange (transfer) stations.The transfer bus stations of two bus services are denoted as Similarly, for coordinating three or more bus services, the transfer stations are: Considering the operations of each service line, EWT is the performance index for regularity-based high frequency bus services operating in dense metropolitan areas.EWT is calculated at each station {S 1,Q j , S 2,Q j , … , S |S Q j |,Q j } of a bus service Q j over the day for identifying the level of service irregularity by calculating the excess waiting times (EWTs) of passengers.The EWT on a bus station S k,Q j for bus service Q j is: ) is the summary of all headways between consecutive trips of service Q j at station S k,Q j during the entire day.The first term of Eq. 4 represents the observed average waiting time of passengers at the bus station over a time period and the second the scheduled waiting time.In most cases, PTAs give different importance to the EWT of bus stations when computing the total EWT of the day for one bus service Q j by adding weight factors at each bus sta- tion according to its importance.As a result, the EWT of one service Q j is: where w i is the associated EWT weight to each station If two or more bus services Q i , Q j intersect at one or more bus stations then the coordination of services requires the optimization of waiting times at the transfer stations while also minimizing the EWT for each independent service Q i , Q j subject to tactical planning constraints such as the required bus frequencies.For passenger transfers from bus trips of service Q j to bus trips of service Q i , the waiting time at each transfer station for each trip R 1 3 Transfer synchronization of regularity-based bus operations The trip R y,Q i that minimizes Eq. 6 is the trip from service Q i that arrives at the trans- fer bus station TS j after trip R x,Q j and its arrival time is closer to trip R x,Q j than any Computing the passenger waiting time of every trip R x,Q j at each transfer station TS j requires to find first its closest trip from the connected bus service Q i and this needs |R Q i | computations.Equation 6 assumes that both services intersect at the same bus stop and does not consider service intersections where passengers have to walk from one bus stop to another as presented in Fig. 1.For computing the waiting time of all bus trips of service Q j a number of tions is required and for expanding the computation to all bus transfer stations TS, the total number of computations is The transfer waiting time at all transfer bus stations TS j ∈ (S Q i ∩ S Q j ) during the day between the bus services Q j , Q i is: Concurrently, the waiting time of passenger transfers from service Q i to service Q j is: From Eqs. 7, 8 one can derive the total transfer waiting time from all transfer stops TS 1 , … , TS |TS| of two bus services Q i , Q j by adding the expected transfer waiting times at each transfer stop.However, a simple aggregation of the transfer waiting times at the transfer stops is not an adequate metric for measuring the total transfer waiting times of passengers since some transfer stops might accommodate significantly more transfers than others.For considering the demand for transfers, the ) is assigned at every stop TS e that accommodates transfers between services Q i and Q j where is the number of daily transfers from service Q i to Q j at transfer stop TS e .This way, the sum of the importance weights of all transfers between two services are equal to one: = 1 and some transfer stops might have higher importance value (i.e., IW Q i ,Q j TS e = 0.7 ) while others a lower one (i.e., IW Q i ,Q j TS f = 0.1 ) indicating a lower number of transfers.The same approach of defining weight factors to bus stops for including the importance of the passenger demand to the waiting times is followed also by the transport authorities that use operational KPIs such as the EWT (refer to Eq. 5).Those importance weights are used to calculate the total transfer waiting times with the consideration of transfer demand after updating Eq. 7: The coordination of two (or more) bus services Q i , Q j is a multi-objective problem with the following conflicting priorities: (1) minimize the EWT for each independent service Q i , Q j ; (2) minimize the passenger waiting times at the transfer stations.Therefore, a trade-off between the improvement of each service and the improvement of passenger transfer times should be established.For this reason, one can utilize weight factors that indicate the importance of each problem objective resulting in the formation of a single scalar objective function.The decision variables of the optimization problem are the rescheduling departure times x of each bus trip.For instance, when coordinating two bus services If we consider also that from the frequency planning phase all consecutive bus trips of service Q i must have a departure time difference in the range ( a 1 , b 1 ) minutes and of service Q j in the range of ( a 2 , b 2 ) minutes for one time period of the day, then the coordination problem can be defined as: 1 3 Transfer synchronization of regularity-based bus operations where W 1 , W 2 , W 3 are weight factors denoting the importance of minimizing the EWT times at services Q i , Q j and the importance of reducing the waiting times of transfer passengers.This problem is a combinatorial problem because the minute changes of the departure times for each trip, x i , can take any integer value from a set q that contains all integer numbers since further denomination of timetable departure times is not acceptable from transport authorities for practical reasons.For instance, a typical timetable (general transit feed specification (GTFS) schedule) of a bus service is following the minutes' denomination.The denomination cannot be reduced to seconds since it is difficult to communicate this information to the bus driver who has to run the service.As might be expected, a further denomination to half-minutes can be applied (i.e., x i = 4.5 min) during the operational control.In such case, the solution set will remain discrete and the options for departure time changes will be increased.Theobjective function of the problem is nonlinear and non-convex while all inequality constraints resulting from the frequency range limits are linear.Because of the non-convexity, it should be noted that the convergence to the global optimum cannot be guaranteed even if the problem was continuous by nature.
The above problem is a nonlinear integer programming problem because the decision variables can take discrete values.A solution method that is widely used for discrete optimization is the Branch and Bound (B&B) method (Land and Doig 1960) that relaxes the discrete optimization problem of Eq. 10 to a series of continuous sub-problems.Given the nonlinearity of the objective function, any continuous sub-problem of Eq. 10 can be solved with sequential (10) quadratic programming (SQP) which is one of the most effective numerical optimization methods for nonlinearly constrained optimization problems (Powell 1978).It cannot be guaranteed, though, that the solution of the SQP is a globally optimal solution of the problem at hand, notwithstanding the fact that SQP converges always to a locally optimal point after starting from an initial solution guess, because the objective function of Eq. 10 is non-convex; thus, it cannot be proved that one of the locally optimal SQP solutions is also a globally optimal one.Due to that, it is required to use exact optimization algorithms that are more computationally expensive (such as the brute-force method) in order to guarantee the convergence to a globally optimal solution.The brute-force search evaluates the performance of all possible solutions.Finding a globally optimal solution with brute-force requires an exponential number of computations Proof Let us consider the coordination of two bus services, Q i , Q j .If we modify the departure time of one trip the required computations for finding the consecutive bus trips that can be used for a potential transfer at any transfer stop are (according to Eq. 6).If the number of transfer stops between services Q i and Q j is |TS|, this cost grows to |TS||R Q j ||R Q i | for computing the transfer waiting times of passengers at all transfer stops after the departure time modification of the examined bus trip.Therefore, for each departure time modification a number of computations is required for re-evaluating the waiting times of passengers at transfer stops.Each examined trip can modify its departure time by receiving any potential value from the set q.This results in a number of |TS||R Q j ||R Q i ||q| computations for evaluating the impact of different departure times of an examined trip at the coordination problem of Eq. 10.Now, if one evaluates the full list of potential combinations of departure time modifications for all trips; then, the number of combinations that should be explored is |q| computations which are required for evaluating the performance of each departure time modification, the number of required computations for evaluating all potential departure time modification options and selecting the globally optimal one is Transfer synchronization of regularity-based bus operations This means that a globally optimal solution cannot be computed even for very small scenarios by evaluating all possible departure time modification combinations.As an example, finding the global optimum re-scheduling solution with simple enumeration for two bus services with

Coordinating bus services with sequential hill climbing after approximating all constraints with exterior point penalties
Due to the computational intractability of the bus coordination problem, we introduce a greedy method based on the hill climbing heuristic search for finding a solution for the bus coordination problem.It cannot be guaranteed though that the solution provided by the proposed heuristic is a globally optimal solution of the bus coordination problem.In addition, its convergence rate to a globally optimal solution cannot be established because a globally optimal solution cannot be calculated with exact numerical optimization methods for realistic problem instances.As a result, the performance of the solution provided by the hill climbing heuristic search is compared against the performance of other heuristic methods.
In the beginning, we introduce an exterior point penalty function, p(x), to approximate the constrained bus coordination problem by an unconstrained one which is structured in a way that its minimization favors the satisfaction of the constraints.This penalty function adds to the objective function f(x) a term that produces a high cost for the violation of any frequency range constraint; thus, providing a direction towards constraints' satisfaction.In this way, the constrained bus coordination problem of bus services Q i , Q j in Eq. 10 is approximated by the following unconstrained one: The penalty function p(x) is structured in such a way that it is equal to the objective function value f(x) if at some point we reach a solution x for which all frequency range constraints are satisfied: ]) 2 is the added term to the objective function of the constrained optimization problem and dictates that if a constraint c i (x) has negative score, then min[−c i (x), 0] = −c i (x) since the constraint is violated for the current set of variables x.In this case, the objective function f(x) is penalized by the term (11) minimize ) 2 which returns a progressively higher penalty the further from zero c i (x) is.If a constraint c i (x) ≥ 0 for the current set of variables x, then this constraint is satisfied and is not penalizing the objective function since min[− c i (x), 0] = 0 .In this way, we allow the progressive penalization of violating constraints c i (x) < 0 (exte- rior point penalization) and we do not penalize any satisfied constraint c i (x) ≥ 0 even if it is an active one c i (x) = 0 (no interior point penalization).We should note here that any constraint violation should penalize significantly the penalty function.For this reason, we included also a weight term W c at the penalty function which ensures that if the penalty term ]) 2 > 0 , the penalty func- tion is significantly affected.
For minimizing the penalty function score, an iterative solution space search is performed based on hill-climbing with the aim of finding repeatedly a new solution x that provides a lower penalty function score than the previous one until no better solution can be found (convergence).Although there is no guarantee that this is the global optimum solution due to the heuristic nature of the search, the algorithmic search tries to return a near optimal solution.Initially, at iteration k = 0 we start with a random selection of decision variable values x k=0 .We can start for example from the "do-nothing" scenario where no re-scheduling is performed yet and Starting from x k=0 which is the initial guess solution with pen- alty function score p(x k=0 ) we try to explore new solutions and replace the existing one in a greedy manner.If we coordinate bus services Q i , Q j we have a set of trips The search starts by selecting the first trip ∈ R Q i ∪Q j and tries to update the current solution x k=0 with a new solution that reduces the value of the penalty func- tion.In this greedy local search we test all possible departure time variations x k=0 from a set q = {− 5, … , 0, … , +5} and select the one that returns the lowest penalty function score.
The set q = {− 5, … , 0, … , +5} min.does not include all integer values, as it should be according to the problem definition of Eq. 10, but a small sub-set of them.After the end of each iteration this departure time change is adopted and we search for new departure time changes in future iterations.Therefore, for a bus trip we might have a calculated departure time change of + 3 min.from the first iteration, + 2 min.from the second and + 4 min.from the third resulting in a total departure time change of + 9 min.if we had three iterations until convergence which would have been higher than the upper bound of set q (= + 5 min.).Therefore, the departure time changes for each trip can be outside those bounds (in fact, they can be any integer number).However, we use those bounds q = {−5, … , 0, … , +5} min. of departure time changes at each iteration for the reduction of the severe computational costs.

3
Transfer synchronization of regularity-based bus operations If we had allowed an unbounded departure time change for each trip at each iteration of the search, then the set q would have contained all integer numbers and the number of computations for finding the optimal departure time change would have prohibited the implementation of the algorithm.Therefore, we allow only the search of a discrete number of limited options from the set q = {− 5, … , 0, … , +5} min.at each iteration and the algorithm is considerably faster without loss of generality because if a bus trip needs to postpone its departure time significantly, this will be calculated sequentially via postponing it little by little at each iteration until we converge to the final solution which cannot be further improved.
After testing all possible values from the set q for trip , we continue sequentially by performing the same procedure for trip + 1 and this continues until we reach the last daily trip of R Q i ∪Q j .Therefore, we have an updated solution x k=1 where p(x k=1 ) ≤ p(x k=0 ) .After that, we update the departure times of all trips considering the departure time changes: After the first iteration a number of con- straints is not satisfied if the penalty function value p(x k=1 ) is still greater than the objective function score f (x k=1 ) .In any case we start a new iteration k = 2 and we perform exactly the same greedy local search for replacing the departure time values of each trip with newer ones that reduce the penalty function score.Then, we update again the departure times of all trips Iterations continue in a similar manner and after some iteration, k = k � , we would reach a point where the penalty function score is equal to the objective function score.This ensures that all constraints are satisfied.From this point onward, any improvement to the penalty function will be a direct improvement to the objective function because before iteration k = k � there was no guarantee that solution updates that reduced the penalty function score reduced also the objective function score.Finally, after a number of iterations, k = k �� > k � , it will not be possible to improve the penalty function value any further.At this point, the hill-climbing-based solution space search has converged to a locally optimal solution and the search terminates.It should be noted here that it cannot be guaranteed that the resulting solution is a globally optimal one.1: function Hill-Climbing heuristic for exterior point penalty function minimization(p(x), f(x)) 2: Set iteration k = 0 and initial random guess x k=0 ; 3: Set q = {−5, ..., 0, ..., +5} minutes; 4: while new iterations reduce the penalty function score: p(x k ) = p(x k−1 ) do 5: Set iteration k←k+1; 6: Set Store the penalty function value p ← p(x k ) 8: for ρ in range R Q i ∪Q j do 9: for m in range {1, ..., |q|} do 10: Set z ← x k ρ and x k ρ ← m; 11: if p(x k ) < p then 12: Set p ← p(x); 13: else 14: Set back x k ρ to its previous value: end if 16: end for 17: end for 18: Update the departure times of each bus trip: 21: end while 22: return solution x k=k where p(x k=k ) = f (x k=k ) if all constraints are satisfied 23: end function The computational complexity of this heuristic search can be expressed with the use of the big O notation which is used for classifying algorithms according to how their running time grows as the input size grows.This heuristic search is expected to converge into a solution in polynomial time after performing |TS||R Q j ||R Q i ||q||x|k ′′ computations which depend on the number of iterations until convergence k ′′ .This computational complexity allows us to coordinate two or more bus services by proposing new departure times for the scheduled trips.
Proof Let us consider the coordination of two bus services, Q i , Q j .As it is already demonstrated, a number of |TS||R Q j ||R Q i | computations is required for evaluating the performance of a departure time modification.In the proposed hill-climbing search, the search is performed in a sequential manner by starting from the first trip and evaluating all potential time modification options from the set q resulting in |TS||R Q j ||R Q i ||q| computations.After the first trip, the same approach is performed for all other trips ∈ R Q i ∪Q j in a sequential order requiring |TS||R Q j ||R Q i ||q||x| com- putations for establishing the departure time modifications of trips that reduce further the cost of the penalty function.After the end of this search, the established departure time modifications are considered as the incumbent problem solution and the same iteration is performed k ′′ times until it is observed that the penalty function 1 3 Transfer synchronization of regularity-based bus operations score cannot be improved any further resulting in a total number of |TS||R Q j ||R Q i ||q||x|k ′′ computations until termination.□ From the above proof, it is evident that the computational cost of the algorithm is polynomial and does not depend only on the size of the problem, |R Q j | + |R Q i | , and the total number of transfer stops |TS|, but also on the required number of iterations until the penalty function cannot be improved any further: k ′′ .The number of required iterations until convergence can vary significantly when solving different problem instances (i.e., different coordination problems).For this reason, a typical approach in numerical optimization is the use of an upper bound that determines the number of maximum allowed iterations until convergence for ensuring the termination of the optimization process even in the case that a solution cannot be found for a specific problem instance.If one uses a threshold value, K max , for the maximum number of iterations k ′′ , then the computational complexity at the worst-case sce- nario can be described with the use of the big O notation as
For deriving the planned schedules, a library was developed in Python 2.7.The library processes .txtfiles and converts/stores them to an SQL database.This facilitates data queries and enables web-based visualization of the public transport operations with the use of OpenStreetMap (via OpenLayers, an open-source JavaScript library: http://www.openlayers .org/api/OpenLayers .js).
In particular, we focused on coordinating two bi-directional central bus services ( Q 1 is bus line 1 and Q 4 is bus line 4) in Stockholm.Each direction of one bi-direc- tional service has another EWT score and set of constraints, but the service-wide EWT can be computed as the average of each direction's EWT.Bus service Q 1 com- prises of direction 1 (Essingetorget to Stockholm Frihamnen) and direction 2 (Stockholm Frihamnen to Essingetorget), bus service Q 4 comprises of direction 1 (Gull- marsplan to Radiohuset) and direction 2 (Radiohuset to Gullmarsplan).
Figure 2 depicts the examined bus lines 1 and 4. In addition, the EWT at each bus station calculated from the planned arrival times of the GTFS data according to Eq. 4 is presented in Fig. 5.The EWT of passengers at each bus station is calculated from 2 p.m.-7:30 p.m. because during the planning frequency stage, the daily operations of services Q 1 , Q 4 were split in four time periods and 2:00 p.m.-7:30 p.m. was one of them to which a frequency range constraint was assigned from the transport authority.
The allowed upper bounds for the headways between successive bus trips at the departure station are defined by the transport operator during the frequency planning phase and are presented in Table 1 for each direction of services Q 1 , Q 4 .Those upper bounds are directly derived from the GTFS trip schedules.During the coordination of services, the departure times of bus trips can vary but they have to satisfy the upper bounds.In addition, to ensure that two successive bus trips do not start at the same time, we introduce a lower headway bound of 1 min at the departure stop.
There are also five ( 5) transfer bus stations of the two bi-directional services  Table 1 Headway upper bounds at the departure stop for the successive bus trips of services Q 1 , Q 4 on both directions for the time period 2:00 p.m.-7:30 p.m.

Service
Headway upper bounds (min) 1 3 Transfer synchronization of regularity-based bus operations stations according to the GTFS data are presented in Fig. 4a.In more detail, in Fig. 4a are presented: (1) the aggregated transfer waiting times of passengers from service Q 1 to Q 4 for direction 1 trips which is equal to WT dir.1 Q 1 ,Q 4 = 717 min for all daily trips according to Eq. 7; (2) the aggregated transfer waiting times from service Q 4 to Q 1 for direction 1 trips which is equal to WT dir.1 = 1808 min for all daily trips; (3) the transfer waiting times from service Q 1 to Q 4 for direction-2-trips? which is equal to = 711 min for all daily trips; (4) the transfer waiting times from service Q 4 to Q 1 for direction-2-trips? which is equal to WT dir.2 Q 4 ,Q 1 = 1525 min for all daily trips.Those total transfer waiting times include all waiting times at each of the five transfer stops.The values of those transfer waiting times are so high because they are the sum of all waiting times at transfer stops during the entire day (not the average waiting time values).From the above, the total daily transfer waiting times from service = 1428 min and from ser- = 3333 min.The big difference to the total transfer waiting time from service Q 4 to Q 1 is due to the lower frequency of service Q 1 .
For optimizing the regularity-based bus service operations of the circular service we utilized the exterior point penalty function, p(x), and the hill-climbing heuristic search strategy described in Sect. 4 for converging to a near optimal solution.The convergence starting point is the planned GTFS schedule ( x i = 0, ∀x i ∈ x ).The planned GTFS schedule EWT scores and transfer waiting times of passengers for services Q 1 and Q 4 are presented in Table 2.In Table 2 the service-wide EWT Q 1 of service Q 1 is the average of the EWT values from direction 1 and direction 2 presented in Fig. 5; therefore, 1∕2(1.053 + 0.84) = 0.947 min.The EWT score is significantly lower than the trans- fer waiting times value because the transfer waiting times are the sum of all waiting times during the day while the EWT is the average value of excess waiting times at stations.Thus, a direct comparison between them requires the use of weight factors.Consequently, the bus service coordination objective function 000433732 for placing the same importance to the EWT scores of each service and to the total daily transfer waiting times.Finally, the exterior point penalty function ]) 2 had initial score p(x) = 30.13because 26 bus trips violated the frequency range limits of Table 1 and started their trip without any departure time difference from their preceding bus trip.Fig. 4 Passenger waiting times at every transfer station for each direction of services Q 1 , Q 4 in minutes (for the time period 2:10 p.m.-7:30 p.m.).The consecutive transfer bus trips at the transfer stations are presented at the horizontal axis, the transfer bus stations at the vertical axis and the transfer waiting time at each station is presented with a color scheme ranging from deep blue (synchronized transfers with 0 min waiting time) to deep red (12 min transfer waiting time) 1 3 Transfer synchronization of regularity-based bus operations The convergence of the hill-climbing search is implemented on a 2556 MHz processor machine with 1024 MB RAM and is presented in Fig. 3 showcasing the coordination objective function, f(x), and penalty function, p(x), improvements after re-scheduling the departure times of bus trips.In Fig. 3 it is demonstrated how the hill-climbing search finds new bus trip departure times that reduce the value of the penalty function.In the horizontal axis is the number of iterative departure time changes, c, that kept decreasing the penalty function score.For each departure time change that decreased the penalty function score, one can note the changes to the EWT, the WT and the coordination objective function f.At the 404th penalty function decrease the penalty function score is for the first time equal to the coordination objective function p(x c=404 ) = f (x c=404 ) = 2.463 .After this point all frequency range constraints are satisfied and any decrease to the penalty function results in an equivalent decrease of the coordination objective function.Finally, the penalty function value stops to decrease further after the 542nd change and we terminate the search assuming that we have reached the approximate global optimum point.At convergence, p(x c=542 ) = f (x c=542 ) = 2.216 , EWT(x c=542 ) = 0.706 min and W 3 × WT(x c=542 ) = 1.509 or WT(x c=542 ) = 3480 min.
The regularity improvement after optimization for services Q 1 , Q 4 in terms of EWT reduction of passengers at stations is presented in Fig. 5.In addition, passenger waiting times at each transfer station after optimization are presented in Fig. 4b showcasing the more evenly distributed waiting times of passengers at transfer stations compared to the previous ones from the planned schedule presented in Fig. 4a.
The original schedules of bus services Q 1 , Q 4 were manually optimized from the bus operator based on ad-hoc rules; therefore, we were able to improve further the service regularity by 50-70% as shown in Fig. 5. Nevertheless, if the original schedules of bus services Q 1 , Q 4 were optimized based on regularity, then any re-sched- uling attempt of reducing the waiting times at transfer stations would have deteriorated the regularity level of those services.For this reason, a sensitivity analysis is important to understand how the transfer waiting time reduction affects the regularity of services and find a critical point where any gains from transfer waiting times improvements are less than the losses due to deterioration of service regularity.For this reason, we changed the value of the transfer waiting times weight factor, W 3 , and optimized again the bus coordination problem starting from W 3 = 0 which returns the optimal service regularity without considering the transfer waiting times.
Results are summarized in Fig. 6 where for W 3 = 0 the EWT of passen- gers for both services took its minimum possible value of EWT = 0.591 min.At this point, the total transfer waiting times of passengers over the entire day was WT = 4204 min.After re-scheduling the bus operations by putting more emphasis on coordination, we reach a point where for W 3 = 0.0002 we have EWT = 0.609 min.and WT = 3641 .In other words, for 2.8% deterioration of the regularity of services Q 1 , Q 4 we reduced the transfer waiting times by ≃ 13% .After this critical point, reducing further the transfer waiting times deteriorates significantly the EWT of individual services leading to an EWT of 0.706 min.for W 3 = 0.0004337 .This 2.8% reduction of the regularity for services Q 1 , Q 4 is well-justified because the transfer waiting times are improved by ∼ 13% and, according to the literature Rietveld et al.   2001), Iseki and Taylor (2010), the transfer waiting times are perceived as even higher than they really are (up to 2-3 times) from the commuters.

Comparison against SoA solution methods
The scope of this analysis is to test the proposed exterior point sequential hill climbing heuristic against other heuristics or analytical optimization methods for determining the trade-off between solution accuracy (in terms of global optimum approximation) and computational costs.This is very important because the regularity-based bus service coordination problem is computationally intractable and high accuracy heuristics are required.At first, exhaustive exact optimization methods -such as simple enumeration-are rejected due to the exponential computational costs.A solution method that can theoretically approximate the global optimum is the Branch and Bound (B&B) that relaxes the discrete optimization problem of Eq. 10 to a series of continuous ones that can be solved with sequential quadratic programming (SQP) given the non-linearity of the objective function.However, Eq. 10 is also non-convex and in practical implementations the SQP method converged to local optima when trying to solve each instance of the continuous regularity-based bus coordination problems.Thus, a multi-start strategy where the SQP is implemented multiple times by utilizing different initial search points for exploring more local optima is adopted.With this strategy, multiple local optima were computed and Fig. 6 Optimization of the regularity-based bus coordination problem for different values of weight factors W 3 showcasing the trade-off between service regularity and passengers waiting times at transfer stations the best-performing one was selected; thus, increasing the chances that the bestperforming local optimum is closer to a globally optimal solution (Fig. 7).
In addition, two other discrete optimization heuristics are tested.The first heuristic is a problem-tailored simulated annealing search.Starting from the penalty function of Eq. 11 we generate a randomly selected initial guess of departure time variations for every bus trip denoted as x {0} .This initial solution guess returns a penalty function score p(x {0} ) and we perform linear cooling by setting the initial tempera- ture to Temp = 1500 according to the following algorithm.

3
Transfer synchronization of regularity-based bus operations return the parent which has the lowest penalty function score from the list of generated parents 27: end function Those two heuristics were selected instead of differential evolution or particle swarm optimization because of the discrete nature of the optimization problem.After optimizing the problem of Eq. 11 with the above-mentioned methods, the results are summarized in Fig. 7.The left sub-figure of Fig. 7 presents the Fig. 7 Summary results comparing the penalty function optimization and the computational costs of heuristic search methods improvement of the initial penalty function score (which was equal to 30.13) after optimization with simulated annealing with linear cooling, sequential GA, B&B with multi-start SQP and sequential hill-climbing demonstrating the significantly better convergence of the last two solution methods to the global optimum.The right sub-figure of Fig. 7 presents the computational times of the above-mentioned solution methods in minutes where the B&B with multi-start SQP had clearly the slowest convergence because of the implementation of the SQP optimizer many times from different starting points.

Discussion
This work addressed the problem of bus operations re-scheduling for the coordination of regularity-based services.First, the regularity-based bus coordination problem was modeled considering the regularity of bus services expressed in passenger EWT at stations, the coordination level expressed in passenger waiting time at transfer stations and operational constraints in the form of frequency ranges.The resulting problem was a constrained, integer nonlinear program with non-convex objective function that hinders the implementation of exact optimization methods for realistic problem instances.For this, a sequential hill-climbing heuristic method was introduced for converging to an approximate global optimum after using an exterior point penalty function to transform the constrained problem to an unconstrained one with the use of penalty terms.
Experimental results in two bi-directional bus services in Stockholm showcased a 50-70% room for improvement in terms of improving the regularity of existing services.In addition, the coordination at the transfer stations was able to improve 15-30% the transfer waiting times of passengers.Finally, we optimized the two bidirectional services based solely on regularity for identifying which is the optimal regularity of services if they are not coordinated and what is the trade-off between regularity deterioration and coordination improvement.Results summarized in Fig. 6 showcased that with 2.8% of regularity deterioration the transfer waiting times can be improved by 12.77% while after that point the coordination gains outweigh the regularity losses.
In addition, the proposed sequential heuristic search for coordinating regularitybased bus services was tested against other heuristic search methods such as (1) simulated annealing with linear cooling, (2) sequential GA and (3) B&B with SQP and a multi-start strategy.The results of optimizing the operations of the same buses were summarized in Fig. 7.In particular, the simulated annealing with linear cooling search method had the worst performance because the generation mechanism of better performing elements was random and was not capable of exploring the solution space in an efficient way.The sequential GA search inherited was not a simple GA implementation, but inherited the advantageous characteristic of the sequential hill-climbing since the crossover and the mutation were conducted in a sequential approach via an element-by-element investigation.The computational cost of this method was higher compared to the sequential hill climbing due to the presence of more than one population generators and it did not return a good approximation to 1 3 Transfer synchronization of regularity-based bus operations the global optimum because it did not search all potential departure time change options at the sequential search but mainly relied on the mutation to explore more options.Finally, the best performing option of B&B with SQP and a multi-start strategy was significantly slower than all other methods.The reason was the implementation of the SQP optimizer many times from different starting points to compute more local optima in order to find the closest one to the global optimum due to the non-convexity of the objective function and the need to run the entire process repeatedly for computing different instances of the B&B method since the regularity-based bus coordination problem is a discrete one.Therefore, its practical implementation is hindered.

Concluding remarks
This work examined the problem of bus operations re-scheduling for the coordination of regularity-based bus services.Summarizing the results from the previous sections, in the concluding remarks we present the lessons learned from this study.
First, the excess waiting times of passengers of existing services can be improved by more than 50% by modifying the departure times of the scheduled trips without considering the synchronization of services.This is in line with the on-site results from the work of Strathman et al. (1999) that used limited dispatching time modifications to improve the service regularity in Portland (transit provider: Tri-met) demonstrating a practical benefit of up to 23% on the coefficient of variation of the bus headways even in the case where the recorded trip travel times from the daily operations differ significantly from the expected ones.
Second, the waiting times at the transfer stops were able to be improved by 15-30% at the examined two-line network of Stockholm demonstrating that modifying the scheduled timetables of trips can have a positive impact to the transfer waiting times.
Third, an optimal synchronization among lines that minimizes the transfer waiting times of passengers has a significantly negative impact to the service regularity because the service regularity should be reduced for achieving an improved synchronization.For instance, from the experiments we learned that the total transfer waiting times can be improved by up to ∼ 13% by deteriorating the service regularity by ∼ 3% .Nevertheless, any further transfer waiting time improvement impacts the service regularity in a disproportionate manner.
Finally, the proposed sequential hill-climbing method has almost similar performance with the method of B&B that employs multi-start SQP and requires only a limited number of computations for converging to a nearly optimal solution; something that enables its use for larger-scale applications in contrast to the B&B which is impractical.
In future research, the proposed coordination method for regularity-based services can be expanded also to other modes (i.e., trains) which operate under a regularity scheme allowing the cooperation between multiple modes.Finally, one can argue that the capacity of the vehicles might have an impact on synchronized transfers if a synchronized transfer leads to overcrowding.However, this analysis requires

Fig. 2
Fig. 2 Examined bus lines in Stockholm-bus line 1 (light blue color) and bus line 4 (purple color) (colour figure online)

Fig. 3
Fig. 3 Penalty function changes until convergence with hill-climbing.The penalty function value converged to the objective function value after 404 changes and convergence achieved after 542 changes

Fig.
Fig.EWT at every station and service-level EWT for every service before and after the heuristic search re-scheduling for the afternoon EWT phase (time period 2:10 p.m.-7:30 p.m.)

Table 2
Summary results showcasing the improvement of service-wide average EWTs and total daily transfer waiting times for services Q 1 , Q 4 after bus operations re-scheduling with the heuristic hill-climbing search