Applied Mathematical Modelling A model for the periodic optimization of bus dispatching times

We model the problem of dispatching time control in rolling horizons following a peri- odic optimization approach reactionary to travel time and passenger demand disturbances. This model provides more ﬂexibility to transport planners allowing them to adjust the bus schedules during the daily operations. We prove that our periodic optimization model is a convex quadratic program, guaranteeing the global optimality of its solution. To reduce the computational burden, we introduce an iterative algorithm that uses gradient approximations to obtain an approximate dispatching solution. The proposed solution method is found to be signiﬁcantly faster than exact optimization approaches for quadratic programming and maintains an (almost) negligible optimality gap in realistic bus operation scenar- ios. Finally, we show that our periodic optimization method outperforms myopic methods that adjust the dispatching time of each bus trip in isolation using operational data from bus line 302 in Singapore.


Introduction
Scheduling the dispatching times of all trips operating in a bus line is a multivariable optimization problem where the dispatching time of every trip is a decision variable (Strathman et al. [1] , Gkiotsalitis and Alesiani [2] ). Typically, its objective function is scalar, f ( x ) , where x is a vector and x i ∈ x is the dispatching time of the i th trip. The computational cost of solving this problem increases with the number of decision variables and the number of bus stops served by the bus line. Therefore, the entire daily schedule cannot be re-adjusted every time we encounter travel time disruptions or unexpected passenger demand fluctuations. This motivates the use of periodic optimization that reschedules only a small fraction of the daily bus trips based on the currently available information and short-term travel time predictions.
There exists a wide range of theoretical models for modifying the scheduled dispatching times of buses or holding the buses at a control point stop to reduce the waiting time variation of passengers (Hickman [3] , Zhao et al. [4] , Delgado et al. [5] , Gkiotsalitis and Cats [6] ). Although there are several lines of research on dynamic bus control (e.g., stop-skipping, speed control), we hereby discuss the works on periodic optimization of dispatching control because they are closer to our study. A distinct line of research acknowledges that the periodic dispatching time control problem cannot be solved in quasi-realtime and embeds proactively slack times in every daily trip to allow more flexibility during the operations (see Zhao et al. [4] , Adamski and Turnau [7] , Daganzo [8] ). Introducing long slack times though requires more buses to maintain the same service frequency, leading to the under-utilization of resources. Wu et al. [15] scheduling (offline) integer program Yes genetic algorithm heuristic Adamski and Turnau [7] rescheduling multiple programs No sub-optimal, simulation-based control Eberlein et al. [18] periodic bus holding non-convex quadratic program No heuristic Hickman [3] single-vehicle holding control convex quadratic program in a single variable

Yes line search
Gkiotsalitis and van Berkum [21] periodic rescheduling non-convex quadratic program No CPLEX when capacity is not considered This study periodic rescheduling convex quadratic program No iterative gradient approximations Other models on timetable optimization strive to ensure that the dispatching times of trips are evenly spaced throughout the day. For instance, Ceder [9] and Ceder et al. [10] strive to achieve a desired even-load level for all buses at their maximum loading point by determining trip dispatching times that do not deviate significantly from the desired even headways. Similarly, Shafahi and Khani [11] and Gkiotsalitis et al. [12] generated timetables with evenly spaced dispatching times incorporating the additional objective of synchronizing passenger transfers. However, such models do not consider the travel time variations during the actual operations and cannot react to real-time changes in the travel conditions (Fayyaz et al. [13] ). Another distinct line of works models the scheduling problem as a stochastic optimization one incorporating the travel time and passenger demand variability when determining the dispatching times of the daily trips at the tactical planning stage (Xuan et al. [14] , Wu et al. [15] ). Despite the above, solving such stochastic models during the actual operations is not computationally feasible. Thus, control methods that approximate the optimal dispatching policy are typically applied (Berrebi et al. [16] ).
Closer to our work are studies that reschedule the original timetable at the operational level with periodic optimization. Adamski and Turnau [7] proposed a simulation support tool for real-time dispatching control. Similar to our work, their dispatching control objective is the elimination of the deviations between the actual operations and the planned schedule. However, the work of Adamski and Turnau [7] requires the performance of simulations inside the optimization process to produce a sub-optimal full state feedback gain controller. Liu et al. [17] used also a modeling and microsimulation approach to improve the efficiency and reliability of bus rapid transit systems. In Eberlein et al. [18] , Koehler et al. [19] and Gkiotsalitis [20] the dispatching control and holding problem is treated as a periodic optimization problem and is modeled as a nonconvex quadratic mathematical program that cannot be solved to global optimality. Those approaches differ from our work, which models the dispatching time control as a convex quadratic program and does not consider bus holding decisions.
In light of the above, our closest prior art is the work of Hickman [3] which also formulates the dispatching control problem as a quadratic convex program that can be solved to global optimality. In Hickman [3] though, the decisions about the holding or dispatching time modifications are made for single vehicles that are about to be dispatched. That is, every bus trip is treated in isolation resulting in single-variable optimization problems that do not consider the interactions among multiple vehicles. The work of Gkiotsalitis and van Berkum [21] is also close to our work and provides a periodic rescheduling formulation for the bus dispatching problem. In Gkiotsalitis and van Berkum [21] , the periodic rescheduling formulation is proved to be non-convex when considering the capacity of buses. Gkiotsalitis and van Berkum [21] propose a model that estimates the dwell times at stops by explicitly considering the passenger alighting times and employ a commercial optimization solver to solve simplified problem instances that do not consider the vehicle capacity limits. In contrast, our work introduces a convex mathematical program for the periodic rescheduling problem by extending the vehicle movement model of Daganzo [8] , and introduces an iterative gradient approximation solution method based on analytic solutions. A summary of the closest prior arts to our work is presented in Table 1 .
From the above studies, there are no periodic optimization works that model the dispatching control problem as a convex optimization problem. Hence, past periodic optimization models cannot be solved to global optimality. An exception is the work of Hickman [3] , which does not solve a periodic optimization problem. Instead, Hickman [3] addresses a simplified, event-based problem where the decision about the dispatching/holding time of each vehicle is made in isolation without considering the trajectories of following trips. This results in myopic, event-based control instead of holistic periodic optimization for all vehicles operating in a time horizon. Since the existing periodic optimization models for dispatching control are not convex (e.g., Eberlein et al. [18] ), their resulting solutions are sub-optimal and one might exhibit high computation costs when trying to improve their optimality gaps. The latter is a main technical challenge when dispatching control needs to be performed in quasi-real-time. This technical challenge and its associated research gap motivate our work: we propose a periodic optimization model that is proved to be convex and can determine the globally optimal dispatching times of trips in rolling horizons. Additionally, we introduce a solution method that can reduce the computation costs by using gradient approximations and improving iteratively the incumbent solution. Extensive numerical experiments demonstrate that the optimality gap of our solution method is negligible even at large-scale scenarios and can result in significant computational cost reductions. The main contributions of our work to the state-of-the-art are summarized below: • the introduction of an iterative gradient approximation solution method for reducing the computational burden of the periodic dispatching control; • the investigation of the efficiency of our solution method against exact solution methods for quadratic programming; • the performance evaluation of our periodic dispatching control approach in terms of passenger waiting times and service regularity using operational data from a high-frequency bus line.
The contributions of this work lie both in the modeling and the solution method phase of the periodic dispatching control problem. Their practical application can improve the regularity of unstable bus operations by re-adjusting the dispatching times of trips at the operational level.

Problem description
Periodic dispatching control can be applied to a bus line with a set of ordered stops S = 1 , 2 , . . . , s, . . . . The periodic dispatching control can be triggered every time a new trip is about to be dispatched. That is, the daily operations can be split into M = 1 , 2 , . . . , m, . . . rolling horizons equal to (at most) the number of daily trips. If j is a trip that is about to be dispatched at the m th rolling horizon, we can determine its dispatching time and the dispatching times of its immediate following trips, N m . N m is a sub-set of the ordered set of daily trips, N , and includes all trips for which the dispatching times should be determined when trip j is about to be dispatched. Note that in extreme cases N m can contain only one trip (namely, trip j ), or all remaining trips of the day, j, j + 1 , . . . , |N | . One should note that the periodic dispatching control might yield different results based on the length of the sub-set N m . Therefore, the length choice of N m is an important decision of the periodic dispatching control (see Eberlein et al. [18] , Sáez et al. [22] ).
To demonstrate the implementation of our periodic control, we start from the first trip of the day. If the pre-determined length of every set N m is n , we determine the dispatching times of all 1 , 2 , . . . , n trips when the first trip of the day is about to be dispatched. A new rolling horizon can start when a new trip is about to be dispatched. For instance, if the second trip of the day is about to be dispatched, we need to determine the dispatching times of trips 2 , 3 , . . . , n + 1 given the realized dispatching time of trip 1. For simplification, when a trip j ∈ N is about to be dispatched, it is re-indexed as trip 1. Similarly, trips j + 1 , . . . , j + n − 1 are re-indexed as 2 , 3 , . . . , n . If trip j is not the first trip of the day, it has a previously dispatched trip j − 1 which is re-indexed as 0. In addition, trip j + n, which is outside of this rolling horizon, is re-indexed as n + 1 . Following this re-indexing, trips 0 and n + 1 do not belong to N m and are the "boundaries" of rolling horizon m because their dispatching times cannot be modified within this rolling horizon.
To determine the optimal dispatching times within each rolling horizon, we need to model the bus operations (e.g., vehicle trajectories, expected arrival times at stops). To this end, we list the main assumptions of our vehicle movement model: (1) The arrivals of passengers at stops are random because passengers do not coordinate their arrival times at stops in high-frequency services (Bartholdi and Eisenstein [23] ); (2) Service supply, which is determined at the frequency settings stage, suffices for satisfying the passenger demand. Thus, all passengers can be served by the first arriving vehicle. This is the most common assumption on the operational control of high-frequency services (Daganzo [8] , Eberlein et al. [18] ); (3) Bus delays due to the dwell times at stops depend predominantly on passenger boardings since alightings are considerably faster. This assumption is used in Daganzo [8] and is supported by several works on dwell time estimation (see the early work of Kraft and Bergen [24] ); (4) The incremental increase of dwell times arising from headway increases is constant. This assumption is used in Daganzo [8] and is a direct result of the works of Hickman [3] , Eberlein et al. [18] that consider constant incremental boarding times; (5) Dwell times are negligibly small so that we do not have additional passenger arrivals during the limited time a bus is waiting at a stop (Hickman [3] ).
Based on assumptions 1-5, we can model the arrival times of buses at stops and the expected headways. The notation of our model is presented in Table 2 .
The variables a j,s , k j,s , h j,s in our nomenclature vary based on our dispatching control decisions, x j , and are defined as follows. The expected arrival time a j,s of a bus trip j ∈ { 1 , 2 , . . . , n } at stop s ∈ { 2 , 3 , . . . , | S|} is where | S | is the cardinality of set S denoting the set's size, δ j + x j is its determined dispatching time, s −1 φ=1 τ j,φ the sum of inter-station travel times from the first bus stop until bus stop s , and s −1 φ=2 k j,φ the sum of dwell times at stops 2 , 3 , . . . , s − 1 .
If overtaking is not allowed, the inter-arrival headway is always non-negative. If an overtaking occurs, one can re-index the impacted trips j − 1 , j at the stop level as in Hickman [3] . This local re-indexing maintains the non-negative headway value.
The dwell times of buses can be affected by the time headway changes because a longer time headway will require from the trailing bus to board a proportionally higher number of passengers. That is, we adopt the dwell time modeling approach of Daganzo [8] . In Daganzo [8] , the dwell time of each trip j at stop s is defined as

Inter-arrival headways considering boundary conditions
In each rolling horizon, the inter-arrival headways between bus trips 1 and 0 are: The inter-arrival headways between any other trip j ∈ N m \ { 1 } and its preceding trip j − 1 are:

Objective function
Our study focuses on high-frequency services that strive to maintain the service regularity. The service regularity can be assessed with several key performance indicators (see Trompet et al. [25] ). For instance, the minimization of the total passenger waiting times has been used as a proxy of the service regularity in Eberlein et al. [18] . Bartholdi and Eisenstein [23] tried to equalize all inter-arrival headways to improve regularity. Daganzo [8] tried to minimize the variance between the actual and the planned (target) headways.
Similar to Daganzo [8] , in this study we try to meet the target headways by minimizing the squared deviation of the inter-arrival headways, h j,s ( x ) , from their target values, h * j,s . The target headway h * j,s between each pair of trips j − 1 , j at stop s is determined at the tactical planning stage. If this target headway is met, the operational regularity is maintained. To reduce the squared deviation of the inter-arrival headways from their target values, one needs to minimize the following objective function: ≥0 is a vector of weight factors w 2 , w 3 , . . . , w |S| . The weight factors allow the bus operator to prioritize the adherence to the target headway at some significant bus stops in the expense of others.

Schedule sliding constraint
In each rolling horizon m ∈ M , we determine the dispatching times of a number of trips 1 , 2 , . . . , n by minimizing the objective function of Eq. (6) . If, however, the solution of this minimization problem results in a major dispatching time delay of the last trip in the rolling horizon, n , this will result in a "domino effect" where all trailing trips in next rolling horizons will have to postpone their dispatching resulting in changes to the vehicle and crew rosters (Gkiotsalitis and Maslekar [26] ).
If the rolling horizon has a slack time ζ ≥ 0 (see Zhao et al. [4] ), we can enforce the dispatching time of its last trip to not exceed this slack time: With the schedule sliding constraint of Eq. (7) we ensure that the dispatching delays do not propagate to the next rolling horizon because they do not exceed the pre-determined slack time. Based on our objective function in Eq. (6) and the associated constraints, our main mathematical program can be succinctly written as: (5) and (7) } (8) where F ( x ) is the feasible set. The equality constraints of Eqs. (4) and (5) simply set the values of headways for a given x and can always be satisfied for any x ∈ R n since the headways are not bounded from any other constraint. The same is not always true for the inequality constraint of Eq. (7) . I.e., n } such that x 0 n > ζ , and thus x 0 ∈ F ( x ) . As shown in Theorem 2.1 , the objective function f ( x ) is quadratic and convex (with respect to x ). Therefore, ( Q ) can be solved to global optimality if the feasible set F ( x ) is non-empty.

Theorem 2.1. ( Q ) is a convex optimization problem, which has a unique global minimizer (if any) with respect to x.
Proof. Note that the feasible set F ( x ) is defined by linear (in)equalities, and thus is a closed polyhedron. Consequently, it is sufficient to prove that to h j,s . Indeed, it is a sum of strictly convex functions.

Analytic solution of the periodic dispatching control problem
To obtain an analytic solution of the multi-variate quadratic program ( Q ), we should derive the first-order conditions. Our solution is based on gradient approximations and we prove that it is the global optimum of ( Then, the deviation of the time headway h 1 , 2 ( x ) from its target headway is expressed as We now define Similarly, we define Finally, we define Then, the deviation of the time headway h j,s from its target headway value h * The core of our analytic solution is that the values of γ s are significantly small ( γ s ≈ 0 , ∀ s ∈ { 2 , 3 , . . . , |S| − 1 } ) so that the product γ s h j,s ( x ) remains constant for a very small change of the dispatching times. In Section 5.3 we will show that even if the values of γ s are not close to zero, a gradient approximation that assumes the product γ s h j,s ( x ) as constant can provide a good initial solution guess and lead to convergence after a few iterations.
or, if the above solution violates the schedule sliding constraint, Note that the inequality constraint of the schedule sliding (x n ≤ ζ ) is not binding when λ = 0 . The complementary slackness condition is λ(x n − ζ ) = 0 which has two solutions: (i) when the schedule sliding constraint is not binding (inactive), its solution is λ = 0 , and (ii) when it is binding, x n = ζ . To find critical point(s), one should derive the first-order conditions of L ( x ) . Expanding L ( x ) using Eqs. (10) , (12) , (14) and (16) yields The critical points of L ( x ) can be found by setting ∇L ( x ) = 0 and satisfying the complementary slackness condition. This yields the KKT system of equations: By solving the above system of equations we obtain the local extrema which are globally optimal solutions if L ( x ) is convex.
To investigate convexity, we compute the Hessian matrix H ∈ R n ×n with elements Noticing that β s ∈S\{ 1 } w s = 1 n and setting ν = n −1 , the Hessian is: which is non-negative for all possible values of x . Thus, L ( x ) is convex. Given the convexity of L ( x ) , the local extrema derived by solving the system of equations in (18) are globally optimal solutions. The system of equations in (18) has two solutions (one when the schedule sliding-related constraint is binding and one when it is not). By setting λ = 0 in the system of equations in (18) we get: Considering the above, the globally optimal solution when the schedule sliding constraint is inactive can be written as: which completes the first part of our proof. Now, if the above solution violates the schedule sliding constraint, we should consider the schedule sliding constraint as binding.
That is, λ > 0 which requires that x n = ζ to meet the complementary slackness condition. For λ > 0, the second solution is derived. In more detail, by setting x n = ζ in the KKT conditions in Eq. (18) we get the following globally optimal solution when schedule sliding constraint is binding: Proof. This can be easily proved by setting ζ = 0 in Theorem 3.1 .
Another interesting case arises when we want to maintain the target headway at only one important bus stop (i.e., a central bus stop with many interconnections). In this case, from Theorem 3.1 follows the corollary: or, if the above solution violates the schedule sliding constraint, Proof. This can be easily proved from Theorem 3.1 by setting all w s values to zero, except w s * > 0 .

Gradient approximation-based iterative algorithm for sensitive dwell times
In Theorem 3.1 , we showed that our analytic solution can provide the global optimum when γ s ≈ 0 , ∀ s ∈ { 2 , 3 , . . . , |S| − 1 } . Daganzo [8] reported empirical γ s values in the range of 0.01 to 0.1 where 0.01 is mostly observed in off-peaks. Values of γ s 0 can have a significant impact on the solution quality because they result in a marginal effect to the dwell times at stops when x changes. As a result, we devise an iterative solution method which uses the output from the analytic solution as a solution guess for the next iteration.
For instance, let us consider the following analytic solution of ( Q ) when γ s ≈ 0 , ∀ s ∈ { 2 , 3 , . . . , |S| − 1 } and the schedule sliding constraint is active: and, as shown in Section 5.3 , exhibits optimality gaps of less than 1% even at large-scale scenarios.
Step 0: Re-index the trips that can change their dispatching times within the rolling horizon as 1 , 2 , 3 , . . . , n following their dispatching order; Step 1: Assume that the schedule sliding constraint is not binding by setting λ = 0 . Set κ ← 0 and initialize solution x κ = { 0 , 0 , . . . , 0 } ; Step 2: Compute the values of c j,s , ∀ j ∈ N m , ∀ s ∈ S \ { 1 } using x κ and calculate the analytic solution x κ+1 : in a sequential order starting from x κ+1 1 while updating the values of c j,s at each sequence; Step 3: If L ( x κ+1 ) ≤ L ( x κ ) , set κ ← κ + 1 and return to (Step 2). Else, proceed to (Step 4) to check the feasibility of x κ ; Step 4: If x κ n ≤ ζ , then STOP because our assumption of λ = 0 is correct. Else, discard the infeasible solution x κ and continue to (Step 5); Step 5: Provided that x κ n > ζ , set κ ← 0 and re-initialize solution x κ = { 0 , 0 , . . . , ζ } ; Step 6: Compute the values of c j,s , ∀ j ∈ N m , ∀ s ∈ S \ { 1 } using x κ and use the analytic solution method to derive solution x κ+1 : . . , n − 1 } in a (reverse) sequential order starting from x κ+1 n −1 while updating the values of c j,s at each sequence.
Step 7: If L ( x κ+1 ) ≤ L ( x κ ) , set κ ← κ + 1 and return to (Step 6). Else, return solution x κ and STOP. It should be noted that in the special case where γ s ≈ 0 , ∀ s ∈ { 2 , 3 , . . . , |S| − 1 } we can still use Algorithm 1 which will return the global minimizer (and not an approximate solution). In that case, the initial solution guess x κ=0 = { 0 , 0 , . . . , 0 } is irrelevant. That is, even if the initial solution guess is changed, Algorithm 1 will return the same result as in Theorem 3.1 and it will always convergence in a single iteration.
An important final note is that even if this algorithm decides about the dispatching time adjustments of all trips at N m , those decisions can be updated every time a bus is about to be dispatched by applying again Algorithm 1 in a "rolled" horizon. For instance, the same algorithm can be applied for the set of trips 2 , 3 , . . . , n + 1 when trip 2 is about to be dispatched to fully exploit the updated information regarding the inter-station travel times.

Comparison periodic dispatching control and one-by-one (myopic) dispatching control
One-by-one (myopic) decision methods decide the dispatching time of a trip j when it is about to be dispatched with the objective of meeting the target headway with its preceding trip j − 1 (Hickman [3] , Fu and Yang [27] ). This decision is made by solving the following minimization problem every time a trip j is about to be dispatched: Therefore, in a rolling horizon with 1 , 2 , . . . , n bus trips the one-by-one decision methods will make n decisions every time the respective trip j is about to be dispatched. In contrast, our periodic optimization method considers the operations of all trips in the rolling horizon when making a decision. To investigate the performance of our holistic periodic optimization approach against the (myopic) one-by-one dispatching control method, both approaches are applied to the following idealized scenario.
We consider a rolling horizon within which three new bus trips (namely 1, 2 and 3) are expected to be dispatched. Trip 0 that precedes trip 1 has been dispatched when our periodic optimization starts. The dispatching time of trip 0 is δ 0 = 0 s. Its expected arrival times at stops s = 2 and s = 3 are: a 0 , 2 = 900 s and a 0 , 3 = 1600 s .
From Fig. 1 one can note that if we decide the dispatching time of each trip independently, the performance deteriorates. The reason is that each bus trip tries to maintain its own target headways and potential problems are transferred to the next trip(s) of the rolling horizon. The problem intensifies if the slack time is smaller (e.g., ζ = 0) because, in that case, the last trip of the rolling horizon has less flexibility to correct the accumulated delays from previous trips.

Time complexity tests
To investigate the time complexity of our approach, we report the computation costs of solving the dispatching time model ( Q ) with the use of (i) the interior-reflective Newton method for quadratic programming that returns a globally optimal solution (see algorithm 5 on p. 1048 in Coleman and Li [28] ), and (ii) our approximate gradient solution method presented in Algorithm 1 . Our experiments are executed in a general-purpose computer with Intel Core i7-7700HQ CPU @ 2.80GHz and 16 GB RAM. All solution methods and external libraries are in Python 3.6. The algorithm of Coleman and Li [28] that is used as a benchmark is implemented with the use of the 'quadprog' library.
To perform the computational tests, we consider rolling horizon optimization scenarios with 3, 4, 5 and 6 trips. In Fig. 2 , we present the computational costs of the interior-reflective Newton algorithm in Coleman and Li [28] and our approximate gradient solution method. Note that the number of recursive calculations that express the bus motion law increases quickly with the number of bus stops and, after a certain point, it becomes time consuming to evaluate the objective function. Therefore, the interior-reflective Newton algorithm in the 'quadprog' solver -which needs to evaluate the performance of the objective function multiple times -exhibits a rapid computation time increase. Hence, it cannot solve program ( Q ) within a reasonable time (i.e., 10 min) for lines with more than ≈ 16 bus stops.
The computational efficiency tests presented in Fig. 2 demonstrate a significant advantage of our solution method because its computation time remains in the range of 0 to 90 s. On the contrary, the self-imposed computational threshold of 10 min allowed the algorithm of Coleman and Li [28] to solve networks with up to 16 stops when we have 3 trips in a rolling horizon, 15 stops when we have 4, 13 stops when we have 5, and 13 stops when we have 6.

Solution quality tests
Using the randomly generated networks from Section 5.2 , we investigate the optimality gap of our approximate gradient solution method with respect to the solution of Coleman and Li [28] which is a global minimizer of ( Q ). We report the Table 4 Optimality gap of our gradient approximation solution method compared to the globally optimal solution of Coleman and Li [28] in multiple idealized networks.
Networks  Table 4 for each randomized network that can be solved to global optimality within the self-imposed computation time threshold of 10 min. In Table 4 , the first two columns indicate the number of trips optimized in the rolling horizon and the number of stops of the idealized scenario. When solving each idealized scenario, we initially set κ ← 0 and we start with a solution guess x κ=0 .
In the third column we present the performance of the solution of Algorithm 1 after completing its first iteration, x κ+1 = x 1 . Column 4 indicates the required number of iterations from the first iteration until the termination of the algorithm. Column 5 indicates the performance of the solution of Algorithm 1 when the algorithm terminates: f ( x * ) .
Column 6 indicates the performance improvement (%) of solution x * compared to solution x 1 from the first iteration of Algorithm 1 . This demonstrates how the use of our analytic solution at each iteration can provide a direction towards solutions that are closer to the global optimum. In column 7, we present the performance of the solution of Coleman and Li [28] , x , which is the global minimizer of ( Q ). Finally, column 8 presents the optimality gap of our solution method which is computed as From Table 4 one can note that the optimality gap remains below 1% at almost all cases (except the case of 5 trips and 8 stops) and it does not exhibit a substantial increase when the number of bus stops or the number of trips increases (e.g., when the scenarios increase in complexity). This demonstrates that one can apply our solution method to large-scale scenarios without a significant optimality loss.
In the numerical experiments in Table 4 , we observe that the number of required iterations of Algorithm 1 tends to increase with the number of stops when the trips in the rolling horizon remain the same. For instance, in the case of 4 trips the required iterations of Algorithm 1 increased from 0 to 4 when the stops increased from 3 to 15. This can be explained because the marginal effect of the neglected γ s values in our gradient approximation is larger when we have more stops. In such case, Algorithm 1 requires more iterations to find dispatching time adjustments that are closer to a globally optimal solution.
To provide an example, in Fig. 3 we present the performance of the solution of Algorithm 1 after its first iteration, f ( x 1 ) , and its final performance after 5 iterations when Algorithm 1 terminates.
Finally, to show the effect of the number of trips on the solution times and solution quality, in Table 5 we present the optimality gap and the computation time of our gradient approximation algorithm compared to the interior-reflective  Newton algorithm of Coleman and Li [28] . Note that in Table 5 we focus on a scenario with 13 bus stops and the exact optimization algorithm of Coleman and Li [28] cannot return a globally optimal solution for more than 6 trips in the rolling horizon due to the self-imposed computation time limit of 10 min.

Case study description
Our case study is the high-frequency, circular bus line 302 in Singapore. Bus line 302 has 22 stops departing from Choa Chu Kang Loop -Choa Chu Kang Int (44,009) and ending at the same stop. A subsidiary of the government of Singapore's Temasek Holdings (SMRT) operates this line and the Land Transport Authority (LTA) monitors its regularity. Normally it starts operating at 05:30 and ends at 00:55. Its route length is 8.1 km and its total travel time typically ranges from 35 to 40 min. Bus line 302 is selected because it is one of the seven high-frequency bus lines in Singapore that are monitored in terms of service regularity and are placed under the Bus Service Reliability Framework (BSRF) (see Leong et al. [29] ). The bus operator's objective is in line with our proposed dispatching control method that aims at improving the service regularity. Under the BSRF framework, bus lines that do not maintain their scheduled headways are penalized, whereas well-performing lines receive monetary incentives (up to 30 0 0$ for every 0.1 min improvement in regularity at the end of each month, as of May 2014).
Bus line 302 satisfies several other assumptions of our work. Its planned service supply is sufficient and the number of stranded passengers is negligibly small. Additionally, its dwell times depend predominantly on passenger boardings (the observed average time for an extra passenger boarding is 1.47 s). The topology of bus line 302 is presented in Fig. 4 .
Bus line 302 is a feeder service that serves residential blocks, schools, and public amenities, connecting them to Choa Chu Kang Town Center and Yew Tee Mass Rapid Transit (MRT) station. Its primary area of service is Choa Chu Kang neighborhoods 5 and 6. Typically, on this bus line operate 12-meter single-decker buses with a seated capacity of 42 passengers and a standing capacity of 33 passengers (75 passengers in total). In this line operate also high capacity, articulated buses because of the high demand from residents at specific times of the day. As presented in Table 6 , the total number of operating trips per day is 245, and the scheduled (target) headways differ among peak/off-peak hours. Note that from 05:30 until 6:30, the service is not headway-based due to the low passenger demand; thus, the scheduled headways in Table 6 Table 7 summarizes the values of the remaining parameters in our case study.

Evaluation of periodic dispatching control and one-by-one dispatching control
In this experimentation, we demonstrate the application of our periodic dispatching time control and the one-by-one dispatching control proposed in the works of Hickman [3] and Fu and Yang [27] . Since both Hickman [3] and Fu and Yang [27] considered also bus holding, our benchmark one-by-one control logic is modified to consider only dispatching time changes. Herein, we describe the implementation of the one-by-one dispatching control and our proposed periodic dispatching control in our case study. Note that both the one-by-one dispatching control and our periodic dispatching control solve convex quadratic programs when a decision needs to be made. This is why the one-by-one control logic is used as a benchmark. Their difference is that the one-by-one dispatching control solves single-variable decision problems where only the dispatching time of a single trip is considered, whereas our approach modifies the dispatching times of all trips N m in a rolling horizon. When implementing the one-by-one dispatching control logic, we decide the dispatching time adjustment of every trip j from the set of 31 trips operating from 6:30-8:30 when trip it is about to be dispatched. The decision of the dispatching time adjustment, x j , is based on the sole objective of meeting the target headways with its preceding trip j − 1 : Note that to calculate the headway h j,s ( x j ) we use the realized/expected arrival time of the already dispatched trip j − 1 at stop s ∈ S \ { 1 } , the expected arrival time of trip j at stop s based on the inter-station travel time expectations, E [ t j,s ] , and its determined dispatching time, x j + δ j . To evaluate the dispatching control strategies in a realistic setting, we use the expected inter-station travel times, E [ t j,s ] , when making a dispatching control decision because we do not know the realized travel times of trips that have not been dispatched when the decision is made (e.g., trip j ). When this decision is implemented, trip j starts to operate based on its realized inter-station travel times (e.g., the ones observed in practice), ˜ t j,s , and not the expected ones, E [ t j,s ] . This provides a realistic evaluation of the dispatching control strategies because our decisions, which were based on future travel time expectations, are evaluated in the actual operations. With this approach, we consider in our evaluation framework the fact that an optimal dispatching decision x j might not maintain its optimality when applied in the actual operations. When implementing our proposed periodic dispatching control logic, we modify the dispatching times of N m trips that are following trip j . Ergo, we do not determine the dispatching time of only the trip which is about to be dispatched, but of all other trips in N m = j, j + 1 , . . . , j + n − 1 . The number of trips in N m plays an important role in the optimization process because if we consider only 1 trip our approach is reduced to the one-by-one dispatching control described in Eq. (21) . To investigate the influence of the number of considered trips in each rolling horizon, n , we apply our approach when n varies from 2 to 7 trips. We implement our periodic optimization every time a trip is about to be dispatched, resulting in 31 control decisions from 6:30-8:30. Our approach is implemented 6 times in total (every time with a different number of considered trips, n , in the rolling horizons). Table 8 presents the results of the evaluation. In Table 8 we report the results of the following key performance indicators when different dispatching control methods are applied: (i) the average squared headway deviation from the target headway (indicating the service regularity); and (ii) the average passenger waiting time.
The operational performance when applying our periodic dispatching control logic results in the following improvements compared to the baseline (e.g., one-by-one dispatching control logic): • improvement up to 21% in terms of service regularity expressed in the mean squared headway deviation of the realized headways from their target values • improvement up to 15% in terms of average passenger waiting times Another interesting finding is the importance of the number of trips in a rolling horizon, n , in the operational performance of the periodic dispatching control method. As expected, if we only consider two trips, N m = 1 , 2 , the improvement Table 8 Performance of the bus operations from 6:30-8:30 when applying the oneby-one dispatching control logic and our proposed periodic optimization control logic for different numbers of trips, n , in the rolling horizon.

Mean squared headway
Average passenger deviation (min 2 ) waiting times (min) potential compared to the one-by-one decision control is limited (only 5% in terms of service regularity). The performance of the periodic dispatching control increases progressively when considering up to five trips in a rolling horizon. After that, the improvement is marginal because trips that will not be dispatched within a short time period do not affect considerably our current dispatching decisions.

Summary and future research direction
In this study, we introduced a mathematical model for the periodic bus dispatching control of high-frequency services and we proved its strict convexity. Additionally, we introduced a gradient approximation solution method that exploits the analytic solution of our main mathematical program in the case where the dwell times are not affected considerably by the headway deviations.
After testing our solution method against the interior-reflective Newton algorithm of Coleman and Li [28] , we showed that the optimality gap of our solution method remains stable in the range of 0-1% for different scenarios with bus lines of up to 16 bus stops and up to 6 trip dispatches in a rolling horizon. The small optimality gap of our approximate gradient solution method is compensated by its significantly lower computation time which remains below 1 min even in cases where the interior-reflective Newton algorithm requires more than 10 min to return a solution. The reason behind this improved computation cost is that, unlike quadratic optimization algorithms that need to evaluate the objective function several times, our solution method uses an easy-to-compute analytic solution at each iteration to converge within a very limited number of iterations. Thus, it can return a good approximation of the globally optimal solution within a limited time.
In addition to the computational experiments, we applied our periodic optimization model for dispatching control using real data from the high-frequency bus line 302 in Singapore. In this application, the dispatching times of a number of trips, N m , were re-optimized every time a new trip was about to be dispatched. The performance of our periodic dispatching control model when applied from 6:30-8:30 in one day of operations was compared against the performance of the oneby-one dispatching control model adapted from the works of Hickman [3] and Fu and Yang [27] . The application results demonstrated an up to 21% improved performance in terms of service regularity and a 15% reduction in passenger waiting times when considering 4 or more trips in each rolling horizon. For fewer trips, the performance of our periodic optimization model deteriorates and becomes exactly equal to the performance of the one-by-one dispatching control when considering only a single trip in each rolling horizon.
In future research, our approach can be expanded to account for the available vehicle capacity and the total travel times of passengers who are using multiple public transport lines for their origin-destination trips.

Limitations
To facilitate the reproduction and extension of our work, we hereby report its main limitations: • rescheduling the dispatching times in rolling horizons can account for mild disruptions that do not require the addition of more vehicles. If the disruptions are severe, there should be changes in the planned service provision (e.g., vehicle and crew schedules); • our rescheduled dispatching times aim to maintain the target headway, which is the main key performance indicator of high-frequency services. In the case of low frequencies, this objective should be replaced by the adherence to the planned arrival times at stops resulting in a different optimization problem; • our approach is suitable if the planned vehicle capacity suffices to accommodate the passenger demand. In contexts where the passenger demand exceeds the vehicle capacity, we will have stranded passengers with increased waiting times that should be considered in our objective function.