This work proposes a periodic bus holding control method where the bus holding times of all run- ning trips are computed simultaneously within each optimization time period; thus, increasing the coordination among running buses to avoid bus bunching. We consider the adverse effects of the bus holding control on the in-vehicle travel times of on-board passengers and formulate holistic bus holding decisions by modeling the bus holding problem as a discrete, nonlinear, constrained opti- mization problem. Given the computational complexity of the bus holding problem, an alternating minimization approach is introduced for computing the optimal holding times at each optimization instance. The performance of the periodic control method is evaluated against the performance of event-based control methods using 5-month automated vehicle location and automated passenger count data from bus line 1 in Stockholm demonstrating an improvement potential of 5% for the in-vehicle travel times and 11% for the service regularity. Keywords : Dynamic Bus Holding; Nonlinear Programming; Periodic Control; Discrete Opti- mization; Bus Bunching Gkiotsalitis and Cats 2 INTRODUCTION Transport authorities have different metrics to evaluate the performance of bus services ranging from safety Key Performance Indicators (KPIs) to planning and operational KPIs. Focusing on the operational KPIs, transport authorities use punctuality-based KPIs to evaluate the performance of low frequency bus services and regularity-based KPIs to evaluate the performance of high fre- quency services (Gkiotsalitis and Maslekar ( 1), Gkiotsalitis and Stathopoulos ( 2)). In the case of high frequency bus services (i.e., bus lines with more than 5 buses per hour), the most common objective of the bus operations is to maintain an appropriate headway among buses in order to ensure that the waiting times of passengers at stops are close to the planned ones. Not having a speciﬁc, punctuality-based arrival time target provides more ﬂexibility to bus operators in applying control measures that can reduce the waiting time variation of passengers. Given the spatio-temporal variations of trafﬁc and passenger demand, some buses are delayed at some speciﬁc segments of the line. These delays have two adverse effects: (i) they increase the headway of the delayed bus from its preceding one; (ii) the delayed bus arrives late at the next stop and, as a result, it needs to accommodate more passengers who are waiting for an extended period of time. This causes a "domino effect" that makes it difﬁcult for a delayed bus to recover and leads to bus bunching. Several metropolises with high-frequency bus services such as London, Singapore, Hong Kong, Barcelona (to name a few) have adopted speciﬁc KPIs to evaluate the regularity of the running services based on the waiting times of passengers (TfL ( 3), LTA ( 4)). For this purpose, real-time Automated Vehicle Location (A VL) information from the Fleet Telematic System (FTS) of each bus is collected and the recorded arrival times of buses at stops are used for evaluating the performance of the services. This constant monitoring of the bus operations has increased the pressure on bus operators to improve the regularity of their services and has led to the adop- tion of real-time control systems that can communicate actuation measures to bus drivers in near real-time. Nevertheless, most control options are currently based on ad-hoc rules (Cats et al. ( 5), Laskaris et al. ( 6)) which are triggered when a speciﬁc event occurs (i.e., a bus arrives late at a stop) and are not capable of ﬁnding an adequate trade-off between the in-vehicle travel times and the waiting times of passengers at stops. RELATED STUDIES In general, the daily operations of a bus line are formed on the basis of an extensive tactical plan- ning process that involves the planning of bus frequencies, crew scheduling and several other con- straints such as the minimum layover times between successive trips and the permitted dispatching headway ranges (Kepaptsoglou and Karlaftis ( 7), Farahani et al. ( 8), Gkiotsalitis and Cats ( 9), Gkiotsalitis and Cats ( 10)). During the tactical planning phase, the dispatching time of each bus trip is deﬁned and the daily timetable is generated. Nevertheless, tactical planning assumes that the trip travel times remain stable during the actual operations or variate very little from their expected values (Gkiotsalitis and Kumar ( 11), Gkiotsalitis and Maslekar ( 12)). This is not the case though in real-world operations and several studies have proposed the introduction of real-time corrective actions, such as bus holding control, that can reduce the waiting time variation of passengers at stops. The plurality of studies on the bus holding problem (Daganzo ( 13), Koehler et al. ( 14), Hick- man ( 15), Eberlein et al. ( 16), Berrebi et al. ( 17), Sánchez-Martínez et al. ( 18), Gkiotsalitis and Maslekar ( 19)) focus on optimizing the bus holding times at a pre-determined sub-set of bus stops, Gkiotsalitis and Cats 3 known as control points, for reducing the waiting times of passengers. Such studies have a pre- deﬁned target for the waiting time of passengers at each stop (i.e., scheduled waiting time) and employ heuristics for computing the optimal bus holding times at some bus stops for adhering to the target waiting times. The work of Bartholdi and Eisenstein ( 20) is an exception because its objective was to reduce the waiting time variance without adhering to a speciﬁc target waiting time. In most of these works, the implications of the bus holding measures to the future trips are not considered and it is assumed that the slack times are long enough to ensure that delays are not propagated. Fu and Yang ( 21) tested two of the most common bus holding strategies: (i) the one- headway-based control where a bus is held at a control point stop if it is decided upon its arrival there that it is too close to its preceding bus and (ii) the two-headway-based control that estimates also the arrival time of the following bus at the same stop and considers the headways between its preceding and its following bus when computing a holding time. In addition, several works have also considered travel time variations when computing the bus holding times (Chen et al. (22), Daganzo ( 13)). For instance, Hickman ( 15) introduced stochasticity at the single holding problem by modeling it as a convex quadratic program in a single variable. Finally, Daganzo and Pilachowski ( 23) and Zhao et al. ( 24) proposed distributed control models where buses act as agents that communicate in real-time to achieve dynamic coordination. To the best of the authors’ knowledge, most of the existing studies on bus holding are based on the assumption that the slack times are long enough so that the bus holding measures do not impact future trips and the respective crew and vehicle schedules (Fu and Yang ( 21), Zhao et al. (24)). In addition, most studies focus on myopic, event-based rules that decide about the holding time of a bus trip based on its headway from the preceding and following buses when it arrives at a control point. Such myopic approaches might fail to capture the repercussions of the bus holding measure to the entirety of the running trips and can be therefore counterproductive (Fu and Yang (21)). Our study is one of the ﬁrst to examine how a set of holistic holding time decisions can be determined while ensuring that the potential repercussions due to the holding time decisions are taken into consideration. For making such holistic holding time decisions, we consider the effects of different holding time combinations and this results in a decision making process where all holding time options for all running buses are evaluated simultaneously. The main contributions of this work to the state-of-the-art are (i) the introduction of a time-window-based model that makes holistic decisions regarding holding times instead of performing a myopic optimization considering the status of neighboring buses; and (ii) the explicit consideration of domino effects due to bus holdings, such as the increase of the in-vehicle travel times and the delay of future trip departures (also known as "schedule sliding"), with the use of a multi-constrained modeling approach. BUS KINEMATIC MODEL FOR A PERIODIC OPTIMIZATION INSTANCE During the daily bus operations, the running buses, N=f1;2;:::;n;:::;jNjg, are expected to serve a number of predeﬁned bus stops. Before proceeding further, it is important to present the main assumptions of this work: Assumption 1 : The number of trips of a bus line are deﬁned during the tactical planning phase and must be fulﬁlled during the actual operations (i.e., no trip cancellations) Assumption 2 : Bus operators can apply bus holding control measures at control stops and cannot Gkiotsalitis and Cats 4 modify the scheduled dispatching times of buses except in the case of emergency (i.e., due to heavy trafﬁc delays) Assumption 3 : Overtaking between buses of the same line is not allowed (a common assumption used in most related works Xuan et al. ( 25), Sun and Hickman ( 26), Chen et al. ( 22)) As it will be demonstrated later on, this work applies bus holdings measures that are not expected to postpone the departures of future trips operated by the same buses. This is a signiﬁcant step towards ensuring that the optimization of the holding times of the running buses will not result in schedule sliding. In addition, this allows the discretization of the daily time horizon into a set of times periods M=ft;t+d;t+2d;:::;t+jMjdgwithin which periodic optimizations of the daily operations are allowed to carry out. The periodic, time-window-based optimization differs from the event-based headway opti- mization proposed in most studies (e.g., Yu and Yang ( 27), Fu and Yang ( 21)). In a periodic optimization holding control, it is assumed that within a time period d, the link travel times of trips are available from a link travel time prediction model and the time period dis short enough so that these travel times are relatively stable. In addition, the periodic optimization control is time-based and an optimization that computes all the holding times of the running buses is applied only at the beginning of the optimization time period m2M. Under the following assumptions, the expected arrival time of a running bus at one future bus stop, s0, during the periodic optimization period m=ft+md;t+ (m+1)dgis: am;n;s0(x) = (t+md) +xbn;m;rn+s01 å i=rntm;n;i+s01 å i=rndm;n;i+s01 å i=rnxm;n;i (1) where am;n;s0is the expected arrival time of trip nat stop s0during the mthoptimization time period, rnis the next bus stop of trip ngiven its current position at the beginning of the mth optimization timeperiod and xbn;m;rnis the expected travel time from the position of trip nat the beginning of the optimization, which is denoted as bn;m, to its next bus stop rn. The expected arrival of the bus at stop s0is affected also by the expected dwell time of trip nat each stop iwhich is denoted as dm;n;iand the holding times at each stop xm;n;i. Obviously, the expected arrival time of a bus trip at a stop s0is equal to the accumulation of the link travel times, the dwell times and the bus holding times from its current position, bn;m, until stop s0plus the timestamp t+mdthat denotes the beginning of the optimization time period m2M. In Eq.1 the link travel times and the dwell times of the bus trip nare deterministic. However, in real-world operations these two parameters have some degree of randomness. Several works such as Shalaby and Farhan ( 28), Mazloumi et al. ( 29), Chien et al. ( 30) have focused on predicting the travel times of bus operations based on bus data with the use of several techniques spanning from time series analysis to artiﬁcial neural networks. The randomness of the link travel times is explicitly considered in our analysis which studies the effects of travel time disturbances to the waiting times of passengers. By using the predicted link travel time and dwell time values at the time period mth, the headway between a bus trip nand its preceding trip n1 at stop scan be estimated as: hm;n;s(x) =am;n;s(x)am;n1;s(x) (2) Gkiotsalitis and Cats 5 where its value depends on the applied bus holding measures x. Finally, assuming random passenger arrivals at stops that follow the uniform distribution (Osuna and Newell ( 31)), the waiting time of passengers who are willing to board bus trip nis: wm;n;s(x) =hm;n;s(x) 2=am;n;s(x)am;n1;s(x) 2(3) MODEL FORMULATION In this section we formulate the mathematical model in order to solve the periodic bus holding optimization problem. The following nomenclature is used for developing the model. Parameters S=f1;2;:::;s;:::;jSjg set of bus stops N=f1;:::;n;:::;jNjg set of daily trips J set of bus stops where bus holdings are allowed un scheduled dispatching time of each trip n2N un scheduled arrival time of each trip n2Nat the terminaljSj kn embedded slack time for each trip n2N rn required resting time for each trip n2Nbefore the dispatch- ing of the next trip operated by the same bus fn the next trip after trip nthat is operated by the same bus d time duration of each periodic optimization time period M=ft;:::;t+jMjdgset of daily periodic optimization time periods m identiﬁer of a periodic optimization time period bn;m location of bus trip nat the start of optimization period m bs;s+1 travel distance between bus stops sands+1 an;s the recorded arrival time of trip nat stop swhich is set equal to zero if such information is not available yet am;n;s expected arrival time of trip nat stop sduring the mthopti- mization time period without considering bus holdings am;n;s(x) expected arrival time of trip nat stop sduring the mthopti- mization time period as a function of the bus holding times hm;n;s(x) headway between bus trips n1 and nat stop sduring the mthoptimization time period dm;n;s expected dwell time of bus trip nat stop sduring the mth optimization time period tm;n;s expected link travel time of trip nbetween consecutive bus stops sands+1 during the mthoptimization time period zs desired waiting time of passengers at stop saccording to the planned schedule Decision Variables x=fxm;n;sg holding times for each trip n2Nat each stop s2Sduring themthoptimization time period Gkiotsalitis and Cats 6 Note that our aim is to ﬁnd the optimal bus holding times, x, for the running buses at each optimization time period. During the mthperiodic optimization period, the active and inactive decision variables can be deﬁned by the following inequalities: xm;n;s=( xm;n;s=0 ifam;n;s<t+mdoram;n;s>t+ (m+1)dors2J xm;n;s2qift+mdam;n;st+ (m+1)dands=2J(4) where Jis the set of control point stops at which a but trip can be held at. Eq.4 includes several inequalities that determine whether a bus holding decision variable xm;n;sis equal to zero (inactive) or can take values from a candidate set q. From Eq.4 it is evident that if the expected arrival time of a bus trip nat one stop sis before the start of the mthoptimization window { t+md} or after the end of the mthoptimization window { t+ (m+1)d}, then there is no need to make a decision for that holding time during this optimization time window. Many studies such as Lee et al. ( 32) have shown that a driver has a minimum response time of'1:5 second and if a driver is requested to hold a bus at one stop for a very small amount of time, he/she will not be able to implement the holding time request with high accuracy. Several works such as Cortés et al. ( 33), Gkiotsalitis and Maslekar ( 34) have proposed to discretize the holding time instructions to 30-second periods (i.e., 0, 30, 60, ...) seconds. The same 30-second discretization for the holding time decisions is used in commercial software as it is described in Cats et al. ( 5). In this study, we propose a 10-second discretization scheme for the holding time options. In such case, the minimum driver response time of '1:5sec. can only affect the holding time suggestion by 15% at most. Slack Time Constraints During the daily scheduling of the bus operations a slack time is allocated to each trip in order to ensure that the time delay of a bus trip due to trafﬁc will not result in a knock-down effect for the departure of the next trip operated by the same bus. If a bus trip, n, is expected to arrive at terminal jSjat time am;n;jSj(x), then the slack time kn2Z+is the extra time that the bus is allowed to remain at the terminal before starting its next trip. The expected arrival time of a running trip nat terminaljSj,am;n;jSj(x), can be directly linked to the slack time knif we consider the scheduled arrival time of trip nat the terminal, un. Satisfying the slack time constraints can be expressed by the following set of inequalities that consider the scheduled arrival time of a trip at the terminal and the associated slack time: am;n;jSj(x)un+kn8n2R (5) where Ris the set of running trips during the mthoptimization period. The inequalities of Eq.5 denote that a bus trip nshould be anticipated to arrive at the last stop before its scheduled arrival time at that stop plus the slack time. Note that if, due to heavy trafﬁc, any of these inequality constraints cannot be satisﬁed (i.e., the examined trip is delayed beyond the slack time buffer); then, Eq.5 does not allow to hold this trip at control point stops. Gkiotsalitis and Cats 7 Travel Time Limits for On-Board Passengers Excessive holding for regularizing the headways of the running buses is inconvenient for the on- board passengers. The obvious side effects are (i) the inﬂicted increases on the on-board passenger travel times due to holdings, and (ii) the inconvenience of on-board passengers that are held at several stops beyond the required dwell time for the boardings/alightings. Let assume that gndenotes the total bus holding time that is permitted for a trip n. This cumulative limit for the holding times, gn, can be used for controlling the maximum amount of holding times that are allowed for one trip and can ensure that the on-board passengers will not have a travel time increase of more than gndue to holdings. Depending on the passenger demand of each trip n2N, the gnvalues can vary (i.e., trips with higher demand might be allowed to operate without signiﬁcant holdings ( gn'0), whereas trips with lower demand might be more ﬂexible and allow more time for holdings). The cumulative holding time limits, gn, per each trip can be provided from the transport authority or derived from historical data analysis regarding the passenger demand levels imposing the following inequality constraints: jSj å s=1jMj å m=1xm;n;sgn8n2N (6) One should note that Eq.6 accumulates all the holding times that are calculated during different periodic optimization time periods because the state-space model of this work is implemented in a rolling-horizon. Bounding the Upper Limits of Bus Holdings As described in Eq.4, each holding decision is allowed to take a value from the discrete set q that contains all potential holding time options with a 10-second difference for improving the probability that a bus driver will be able to implement accurately this holding option. However, each holding time should have an upper limit because one bus cannot be held at one stop for a very long time since (i) this increases the on-board passengers’ and driver inconvenience; (ii) other buses from the same or different lines should use the same stop for boardings/alightings. Several studies such as Cortés et al. ( 33), Gkiotsalitis ( 35) and Cats et al. ( 5) have proposed to impose an upper limit to the bus holding time which is in the range of 60-90 seconds depending on the characteristics of the control point stop. Adopting such upper limit for each holding time, the admissible set of holding times is modiﬁed to: q=f0;10;20;30;40;50;60;70;80;90gseconds (7) Reducing the Variance of the Actual Passenger Waiting Times from the Planned Waiting Times Bus operators strive to improve the service reliability by minimizing the passenger waiting time variance from the planned waiting time, zs, at each stop s2S(Daganzo ( 13), LTA ( 4)). If the scheduled waiting time zsat one stop is different than the actual waiting time of passengers, then the service reliability deteriorates and the bus operators should take corrective actions. We illustrate this using an example where the trips R=fn;n+1;n+2;n+3gare the running trips at the start of the mthperiodic optimization instance as they are presented in Fig.1. In addition, letJ=fs=3;s=5;s=7;s=9gbe the set of the control point stops of the bus service in Fig.1. Gkiotsalitis and Cats 8 s=1s=2s=3s=4s=5s=6s=7s=8s=9s=10 Timen n+1 n+2 n+3Control Point Stops n n+1 n+2 n+3 Running BusesBus Stops current time instance FIGURE 1 : Illustration of the running buses at the beginning of the mthperiodic optimization instance If trips R=fn;n+1;n+2;n+3gare the running trips at the mthoptimization instance, then trips Rprior =f1;2;:::;n1ghave already completed their service prior to the start of the mth optimization time period. In this way, we can split the daily trips N=f1;:::;n;:::;jNjginto three sets during the mth optimization instance: (i) the trips that are already completed Rprior =f1;2;:::;n1g; (ii) the running trips that are expected to operate during the mthtime period R=fn;n+1;n+2;n+3g; and (iii) the remaining trips, Rf uture =fn+4;n+5;:::;jNjg, until the end of the day. Having deﬁned the set of running trips, we proceed one step further by identifying the set of trips that are expected to arrive at each bus stop s2Sduring the mthoptimization time period. For each bus stop s2Swe initialize a set Rm;sthat contains all the bus trips that are expected to arrive at that stop during the mthoptimization time period. As a result, at the start of the periodic optimization each trip n2Nis assigned to a set Rm;saccording to the following inequalities: n=( n2Rm;sifan;jSj=0 and t+mdam;n;st+ (m+1)dandn2R n=2Rm;sotherwise(8) where an;jSj=0 denotes that trip nhas not been completed at the start of the mthoptimization time period and t+mdam;n;st+ (m+1)ddenotes that the expected arrival time of trip nat stop slies within the range of the mthoptimization time period where: am;n;s= (t+md) +xbn;m;rn+s1 å i=rntm;n;i+s1 å i=rndm;n;i8n2R (9) In Eq.9, am;n;sis the expected arrival time of trip nat stop sduring the mthoptimization time period without considering bus holding measures, rnis the next bus stop of trip ngiven its current position at the beginning of the mthoptimization timeperiod and xbn;m;rnis the expected travel time from the position of trip nat the beginning of the optimization, which is denoted as bn;m, to its next bus stop rn. After deﬁning the set of trips Rm;sthat are expected to arrive at stops during the mthopti- mization time period, the actual waiting time variance from the scheduled waiting times at stops, zs, should be minimized. Assuming random passenger arrivals at stops, those arrivals follow the Gkiotsalitis and Cats 9 uniform distribution and the waiting times of passengers are half the headway between two con- secutive bus trips (Osuna and Newell ( 31)) resulting in the following objective function: minxf(x) =å s2Så n2Rm;sam;n;s(x)am;n1;s(x) 2zs2 (10) where am;n;s(x) = (t+md) +xbn;m;rn+s1 å i=rntm;n;i+s1 å i=rndm;n;i+s1 å i=rnxm;n;i8n2R (11) From the above descriptions, a mathematical model that (i) minimizes the passenger waiting time variation at stops, (ii) satisﬁes the slack time constraints and (iii) limits the travel times of on- board passengers according to the on-board passenger demand levels is formulated by combining equations 10, 5, 6 for the optimization of the mthrolling horizon: minxf(x) =å s2Så n2Rm;sam;n;s(x)am;n1;s(x) 2zs2 subject to am;n;jSj(x)un+kn;8n2R jSj å s=1jMj å m=1xm;n;sgn;n=1;:::;jNj: xm;n;s=( 0 ifam;n;s<t+mdoram;n;s>t+ (m+1)dors2J xm;n;s2qift+mdam;n;st+ (m+1)dands=2J q=f0;10;20;:::;90gseconds(12) Modeling the Dwell Times The dwell time dm;n;sof a bus trip nat stop sduring the mthoptimization time period can be associated with the number of boardings and alightings. Detailed studies such as Levinson ( 36), Guenthner and Sinha ( 37) and Bertini and El-Geneidy ( 38) showed that there is a linear relationship between the dwell time and the number of boardings where the dwell time can be equal to a minimum value p0that can range from 2-5 seconds in case of no boardings and it increases by p1=1.5-4.5 seconds for each extra passenger boarding depending on the fare payment structure. In accordance with the above studies, the underlying relationship between the observed dwell time and the number of boardings can be expressed by a linear equation: dm;n;s=p0+p1Bn;s (13) where p0is expressed in seconds, p1in seconds per passenger and Bn;sis the expected number of passengers that will board on trip nat stop s. If we assume a passenger arrival rate Bsthat expresses the number of passengers per second that are expected to arrive at stop s, the number of passenger boardings, Bn;s, on one trip nat Gkiotsalitis and Cats 10 stop sis associated with the headway between trip nand its preceding trip n1 at that stop. Consequently, the boarding time of passengers dm;n;sis: dm;n;s=p0+p1Bs(am;n;sam;n1;s) (14) Eq.14 denotes that the dwell time at one stop is proportional to the headway between two consecutive bus trips n1;nmultiplied by the passenger arrival rate at that stop. As a result, holding a running bus n1 at stop simpacts the dwell time of trip nat stop s: dm;n;s(x) =p0+p1Bs(am;n;s(x)am;n1;s(x)) (15) Plugging Eq.11 and Eq.15 into the model of Eq.12 yields: min xf(x) =å s2Så n2Rm;s1 2 (t+md) +xbn;m;rn+s1 å i=rntm;n;i+s1 å i=rndm;n;i(x) +s1 å i=rnxm;n;i 1 2 (t+md) +xbn1;m;rn1+s1 å i=rn1tm;n1;i+s1 å i=rn1dm;n1;i(x) +s1 å i=rn1xm;n1;i zs2 s.t. am;n;jSj(x)un+kn;8n2R jSj å s=1jMj å m=1xm;n;sgn;n=1;:::;jNj: xm;n;s=( 0 ifam;n;s<t+mdoram;n;s>t+ (m+1)dors2J xm;n;s2qift+mdam;n;st+ (m+1)dands=2J dm;n;s(x) =p0+p1Bs(am;n;s(x)am;n1;s(x)) q=f0;10;20;:::;90gseconds (16) Remark: From the model of the mthperiodic optimization described in Eq.16, one can note that the dwell time of a trip at a stop is affected by the holding times. This results in an open-loop relationship between the arrival times of trips at stops and the dwell times because the dwell time modiﬁcations impact the arrival times of trips at stops according to Eq.15 and the arrival time modiﬁcations are modifying again the dwell times expressed in Eq.11 leading to a recursive opti- mization problem with a non-convex objective function. SOLUTION METHOD Decoupling the Arrival and Dwell Times with Alternating Minimization The authors utilize the alternating minimization method for non-convex optimization by solving iteratively the model of Eq.18 where at each iteration the dwell time values are assumed to be deterministic resulting in a simpler sub-problem with a convex objective function. The alternating minimization method initializes the problem by assuming that the optimal bus holding times are equal to zero. For this initial solution, xk, where xk m;n;s=0;8xk m;n;s2xk, the resulting dwell times Gkiotsalitis and Cats 11 are computed: dm;n;s=p0+p1Bs(am;n;s(xk)am;n1;s(xk)) (17) Then, Eq.18 computes the new optimal solution xk+1under the assumption of ﬁxed dwell time values which were computed by Eq.17. After that, the dwell time values are updated for the bus holding solution xk+1and the model of Eq.18 is solved iteratively until a termination criterion is satisﬁed. min xk+1f(xk+1) =å s2Så n2Rm;s1 2 (t+md) +xbn;m;rn+s1 å i=rntm;n;i+s1 å i=rndm;n;i+s1 å i=rnxk+1 m;n;i 1 2 (t+md) +xbn1;m;rn1+s1 å i=rn1tm;n1;i+s1 å i=rn1dm;n1;i+s1 å i=rn1xk+1 m;n1;i zs2 s.t. am;n;jSj(xk+1)un+kn;8n2R jSj å s=1jMj å m=1xk+1 m;n;sgn;n=1;:::;jNj: xk+1 m;n;s=( 0 ifam;n;s<t+mdoram;n;s>t+ (m+1)dors2J xk+1 m;n;s2qift+mdam;n;st+ (m+1)dands=2J dm;n;s=p0+p1Bs(am;n;s(xk)am;n1;s(xk)) q=f0;10;20;:::;90gseconds (18) The alternating minimization approach is formed in such a way that, even if the main problem is non-convex, each alternating step in isolation is tractable. For solving this problem, one can deploy the branch-and-bound (B&B) method that is proposed by Land and Doig ( 39) and is one of the most common approaches for solving discrete, NP-Hard problems. The B&B method explores branches of a dynamically generated tree which represent subsets of the solution set of discrete options. In contrary to the brute-force algorithm that explores the solution space exhaustively, the B&B method uses upper and lower bounds for pruning the set of candidate solutions. In such way, if the estimation of the upper and lower bounds is effective, the exploration effort is signiﬁcantly reduced. In order to deﬁne the lower bound (LB) of our problem, we allow the holding times to take any real value from 0 to 90 seconds instead of values from the discrete set q=f0;10;20;:::;90g. Then, after expanding the problem of Eq.18 by plugging Eq.11 into Eq.18, it becomes a contin- uous Quadratic Programming (QP) problem that can be solved in polynomial time with an exact optimization algorithm such as the Active Set method (for more information regarding the Active Set method, one can refer to Murty and Yu ( 40) or Nocedal and Wright ( 41)). By solving the continuous QP, the lower bound of the problem is established because any dis- crete solution has inferior performance compared to the global optimum solution of the continuous Gkiotsalitis and Cats 12 Starting stopControl Points FIGURE 2 : Eastbound direction of Bus Line 1 in Stockholm that serves 33 bus stops from which 4 are control point stops problem and the steps of the B&B algorithm can be performed1. CASE STUDY: BUS LINE 1 IN STOCKHOLM Application Our application uses data from the eastbound direction of bus line 1 in Stockholm. Line 1 connects the main eastern harbor to a residential area in the western part of the city through the commercial center. In this work we processed detailed 5-month datasets of Automatic Vehicle Location (A VL) and Automatic Passenger Counting (APC) data for this line that serves 33 bus stops and has 4 control point stops as presented in Fig.2. The 5-month A VL data is used for estimating the link travel times of the 118 daily trips that operate in the eastbound direction (direction 1) from 07:00 until 19:00. Note that the time-period 07:00-19:00 is selected as the main focus area of the daily operations since trips prior to 07:00 and after 19:00 do not experience signiﬁcant spatio-temporal link travel time variations due to road trafﬁc. In this work, we perform a time-window-based bus holding control on the daily operations where the entire day is split into 10-minute periods within which the optimization is performed. For applying the computed bus holding measures, we build the Stockholm city network in the open source Simulation of Urban MObility (SUMO) using information from OpenStreetMaps. The daily A VL data from the actual bus operations was used for replicating the daily operations in SUMO. For performing such action, we updated automatically the parameters of the trafﬁc scenarios with the use of the Constrained Optimization BY Linear Approximation (COBYLA) algorithm from the SciPy library in Python that computes the optimal parameter values of trafﬁc 1since the B&B algorithm is a well-known solution method, we will not provide here its detailed steps. Instead, we refer the interested reader to the work of Land and Doig ( 39)) Gkiotsalitis and Cats 13 ﬂows and trafﬁc signal cycles which minimize the root-mean-square error (RMSE) between the actual and simulated arrival times of buses at stops. The vehicle parameters of the running buses (i.e., acceleration, max. speed, deceleration, vehicle length) are speciﬁed in SUMO using an .xml ﬁle that follows the guidelines of SUMO ( 42). After the calibration stage, the acceleration, deceleration and maximum speed values are set as 2.61 m/s2, 4.68 m/s2and 100.0 km/h, respectively. In addition, the daily boarding rates, expressed as number of boarding passengers per hour, are derived from the APC data. For demonstration purposes, the hourly boardings from 08:00 until 12:00 are presented in Fig.3. 15101520253033 Bus Stop050100150200Number of Passengers 08:00-09:00 09:00-10:00 15101520253033 Bus Stop050100 10:00-11:00 11:00-12:00 FIGURE 3 : Hourly passenger boardings rates at each bus stop (from 08:00 until 12:00) The results of the bus holding control for the eastbound direction of line 1 in Stockholm are presented in Fig.4 where the computed bus holding measures at every 10-minute optimization horizon are implemented in SUMO (from 07:00 until 19:00 are performed jMj=72 periodic op- timizations and the respective bus holding times are applied in SUMO). In the same ﬁgure, the results of the one-headway-based control method (Fu et al. ( 43)) are presented where each bus trip that arrives at one control point stop is held there according to its headway from its preceding trip and the variance of the passenger waiting times from their target values is evaluated every 10 minutes. In this way, the proposed periodic bus holding optimization approach is tested against the one-headway-based control method which is one of the most used event-based approaches. Gkiotsalitis and Cats 14 07:00 07:30 08:00 08:30 09:00 09:30 10:00 10:30 11:00 11:30 12:00 12:30 13:00 13:30 14:00 14:30 15:00 15:30 16:00 16:30 17:00 17:30 18:00 18:30 18:5001020304050St. Deviation of the Passenger Waiting Times per trip in secondsOne-Headway-based Control Time-Window-based Control One- Headway- basedTime- Window- based30003100320033003400Average In-Vehicle Travel Time in seconds3303 3140 FIGURE 4 : Results of the Bus Holding control in the eastbound direction of line 1 in Stockholm Sensitivity to Travel Time Variation Finally, Fig.5 presents different instances of the results of the one-headway-based bus holding measures and the periodic optimization bus holding measures when the actual link travel times exhibit noise levels of 10%, 20%, 50% and 80%. An important observation from Fig.5 is that the bus holding measures for noise levels of 10% maintain a passenger waiting time variance close to the desired passenger waiting times whereas when the noise levels are in the range of 80% the waiting times variate signiﬁcantly from their target values due to bus bunching. 10% T-W-b10% O-H-b20% T-W-b20% O-H-b50% T-W-b50% O-H-b80% T-W-b80% O-H-b St. Dev. as a percentage of the Estimeted Link Travel Times (%)0100000200000300000400000500000600000Passenger Waiting Time VarianceT-W-b: Time-Window-based optimization O-H-b: One-Headway-based optimization 123456789 Running Trips050100150200250300350400450Passenger Waiting Times62th Simulation σ=10% of the est. Link Travel Times Desired Waiting Time One-Headway-based Optimization Periodic Optimization 123456789 Running Trips050100150200250300350400450350th Simulation σ=80% of the est. Link Travel Times Desired Waiting Time One-Headway-based Optimization Periodic Optimization FIGURE 5 : Comparative Analysis of the periodic optimization and the one-headway based ap- proach in the presence of travel time noise. The overall daily variance of the actual passenger waiting times from the expected ones is measured in seconds Gkiotsalitis and Cats 15 Discussion Our work replicates the daily operations of the eastbound direction of bus line 1 in Stockholm using the A VL and APC data from the daily operations for building a realistic simulation scenario in SUMO by minimizing the root-mean-square error (RMSE) between the actual and simulated arrival times of buses at stops. This calibration enabled us to apply the periodically optimized bus holding times to a realistic simulation scenario that is close to the real operations. In this simulation scenario, this work evaluated the performance of the bus holding control measures computed by the proposed time-window-based optimization method. The one-headway- based bus holding method was also applied in the same scenario resulting in higher waiting time variance from the desired passenger waiting times compared to the time-window-based opti- mization method. The main reason behind the improved performance of the time-window-based method compared to the event-based methods is that the proposed time-window-based optimiza- tion method evaluates the bus holding options holistically within each time window; thus, enabling the coordination of the running buses. The improvement with the use of the time-window-based optimization method is in the range of'017% as it is presented in Fig.4 (an improvement of 11% on average). The average waiting time difference of passengers at stops from the desired waiting time values was only 3 seconds in some cases whereas in others increased up to 58 seconds (i.e., during the time period 10:00-10:10). The main ﬁndings of the analysis are: 1. the consideration of the slack time and in-vehicle travel time constraints in the optimiza- tion process of the time-window-based method can improve the waiting time variance of passengers (by up to '17%) while improving also the average in-vehicle travel times (by'5%); 2. the performances of the time-window-based control and the one-headway-based control deteriorate following a similar trend when the external disturbance is signiﬁcant and the estimated link travel time values are not close to the actual ones (travel time noise levels of 80%). CONCLUSION This study proposed a time-window-based optimization method for improving the waiting time variance of passengers from the planned waiting time values with the use of bus holding control. Past research works Cats et al. ( 5), Fu and Yang ( 21) have shown that event-based bus holding control methods and bus holding methods that hold the bus trips at more than two control point stops can be effective in reducing the passenger waiting time variance only up to a certain level. For this reason, this study investigated whether a time-window-based approach that computes the bus holding times of all running trips within a time period in a holistic manner can reduce further the passenger waiting time variance from the desired waiting time values. The proposed time-window-based bus holding method was evaluated against the one-headway- based method and exhibited an improved performance with improvements of the passengers’ wait- ing time variance of up to 17% in some speciﬁc time periods for bus line 1 in Stockholm where the control measures computed within each 10-minute optimization period were applied to a simula- tion scenario in SUMO. This simulation scenario replicated the actual operations by using the A VL and APC data from the daily operations for its calibration and the application of the bus holding Gkiotsalitis and Cats 16 measures computed by the time-window-based optimization demonstrated an improvement poten- tial of'5% for the in-vehicle travel times and 11% (on average) for the waiting time variance of passengers at stops. Although the effect of the travel time noise was considered in the validation phase, the de- veloped time-window-based model does not consider the travel time stochasticity as a problem variable when computing the bus holding measures for improving the robustness of the solution. Works in event-based bus holding control such as Hickman ( 15) have considered the stochasticity of travel times in simple holding models that hold a bus at a single control point stop (the departure stop) because of the computational complexity. This challenge however is a very interesting topic for future investigation in order to determine the improvement potential of robust time-window- based control. Author Contribution Statement The authors conﬁrm contribution to the paper as follows: study conception and design: K. Gkiot- salitis, O. Cats; data collection: O. Cats; analysis and interpretation of results: K. Gkiotsalitis, O. Cats; draft manuscript preparation: K. Gkiotsalitis, O. Cats. All authors reviewed the results and approved the ﬁnal version of the manuscript.